Geant4_10
CexmcHistoManager.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /*
27  * ============================================================================
28  *
29  * Filename: CexmcHistoManager.cc
30  *
31  * Description: histograming manager (singleton)
32  *
33  * Version: 1.0
34  * Created: 26.11.2009 21:00:03
35  * Revision: none
36  * Compiler: gcc
37  *
38  * Author: Alexey Radkov (),
39  * Company: PNPI
40  *
41  * ============================================================================
42  */
43 
44 #ifdef CEXMC_USE_ROOT
45 
46 #include <iostream>
47 #include <iomanip>
48 #include <TH1.h>
49 #include <TH1F.h>
50 #include <TH2F.h>
51 #include <TH3F.h>
52 #include <TFile.h>
53 #include <TObject.h>
54 #include <TCollection.h>
55 #include <TDirectory.h>
56 #include <TString.h>
57 #include <TRegexp.h>
58 #ifdef CEXMC_USE_ROOTQT
59 #include <TCanvas.h>
60 #include <QApplication>
61 #include <QFont>
62 #endif
63 #include <G4LogicalVolume.hh>
64 #include <G4Box.hh>
65 #include <G4Tubs.hh>
66 #include "CexmcHistoManager.hh"
67 #include "CexmcHistoManagerMessenger.hh"
68 #include "CexmcProductionModel.hh"
69 #include "CexmcPhysicsManager.hh"
70 #include "CexmcRunManager.hh"
71 #include "CexmcSetup.hh"
72 #include "CexmcException.hh"
73 #include "CexmcHistoWidget.hh"
74 
75 
76 namespace
77 {
78  const G4double CexmcHistoBeamMomentumMin( 0.0 * GeV );
79  const G4double CexmcHistoBeamMomentumMax( 1.0 * GeV );
80  const G4double CexmcHistoBeamMomentumResolution( 0.5 * MeV );
81  const G4double CexmcHistoTPResolution( 0.1 * cm );
82  const G4double CexmcHistoTPSafetyArea( 1.0 * cm );
83  const G4double CexmcHistoMassResolution( 1.0 * MeV );
84  const G4double CexmcHistoEnergyMax( 1.0 * GeV );
85  const G4double CexmcHistoEnergyResolution( 1.0 * MeV );
86  const G4double CexmcHistoMissEnergyMin( -0.1 * GeV );
87  const G4double CexmcHistoMissEnergyMax( 0.2 * GeV );
88  const G4double CexmcHistoMissEnergyResolution( 0.2 * MeV );
89  const G4double CexmcHistoAngularResolution( 0.5 );
90  const G4double CexmcHistoAngularCResolution( 0.001 );
91  const G4int CexmcHistoCanvasWidth( 800 );
92  const G4int CexmcHistoCanvasHeight( 600 );
93 }
94 
95 
96 CexmcHistoManager * CexmcHistoManager::instance( NULL );
97 
98 
99 CexmcHistoManager * CexmcHistoManager::Instance( void )
100 {
101  if ( instance == NULL )
102  instance = new CexmcHistoManager;
103 
104  return instance;
105 }
106 
107 
108 void CexmcHistoManager::Destroy( void )
109 {
110  delete instance;
111  instance = NULL;
112 }
113 
114 
115 CexmcHistoManager::CexmcHistoManager() : outFile( NULL ),
116  isInitialized( false ), opName( "" ), nopName( "" ), opMass( 0. ),
117  nopMass( 0. ), verboseLevel( 0 ),
118 #ifdef CEXMC_USE_ROOTQT
119  rootCanvas( NULL ),
120 #endif
121  messenger( NULL )
122 {
123  for ( int i( 0 ); i < CexmcHistoType_SIZE; ++i )
124  {
125  histos.insert( CexmcHistoPair( CexmcHistoType( i ),
126  CexmcHistoVector() ) );
127  }
128 
129  messenger = new CexmcHistoManagerMessenger( this );
130 }
131 
132 
133 CexmcHistoManager::~CexmcHistoManager()
134 {
135  if ( outFile )
136  {
137  outFile->Write();
138  outFile->Close();
139  }
140 
141  /* all histograms will be deleted by outFile destructor! */
142  delete outFile;
143 #ifdef CEXMC_USE_ROOTQT
144  delete rootCanvas;
145 #endif
146  delete messenger;
147 }
148 
149 
150 void CexmcHistoManager::AddHisto( const CexmcHistoData & data,
151  const CexmcAngularRange & aRange )
152 {
153  G4String fullName( data.name );
154  G4String fullTitle( data.title );
155  G4String rangeTypeLabel;
156  G4String triggerTypeLabel;
157  G4String decorTriggerTypeLabel;
158 
159  if ( data.isARHisto )
160  {
161  if ( data.isARRec )
162  {
163  fullName += "_arrec";
164  rangeTypeLabel = "rec";
165  }
166  else
167  {
168  fullName += "_arreal";
169  rangeTypeLabel = "real";
170  }
171  }
172 
173  switch ( data.triggerType )
174  {
175  case CexmcTPT :
176  fullName += "_tpt";
177  decorTriggerTypeLabel = " --tpt--";
178  fullTitle += decorTriggerTypeLabel;
179  triggerTypeLabel = "tpt";
180  break;
181  case CexmcEDT :
182  fullName += "_edt";
183  decorTriggerTypeLabel = " --edt--";
184  fullTitle += decorTriggerTypeLabel;
185  triggerTypeLabel = "edt";
186  break;
187  case CexmcRT :
188  fullName += "_rt";
189  decorTriggerTypeLabel = " --rt--";
190  fullTitle += decorTriggerTypeLabel;
191  triggerTypeLabel = "rt";
192  break;
193  default :
194  break;
195  }
196 
197  CexmcHistosMap::iterator found( histos.find( data.type ) );
198 
199  if ( found == histos.end() )
200  throw CexmcException( CexmcWeirdException );
201 
202  CexmcHistoVector & histoVector( found->second );
203 
204  if ( data.isARHisto )
205  {
206  G4bool dirOk( false );
207 
208  if ( outFile )
209  {
210  dirOk = gDirectory->Get( fullName ) != NULL;
211 
212  if ( ! dirOk )
213  dirOk = ( gDirectory->mkdir( fullName, fullTitle ) != NULL );
214 
215  if ( dirOk )
216  gDirectory->cd( fullName );
217  }
218 
219  std::ostringstream histoName;
220  std::ostringstream histoTitle;
221  histoName << data.name << "_r" << aRange.index + 1 << rangeTypeLabel <<
222  "_" << triggerTypeLabel;
223  histoTitle << data.title << " {range " << aRange.index + 1 <<
224  rangeTypeLabel << " [" << std::fixed <<
225  std::setprecision( 4 ) << aRange.top << ", " <<
226  aRange.bottom << ")}" << decorTriggerTypeLabel;
227  CreateHisto( histoVector, data.impl, histoName.str(), histoTitle.str(),
228  data.axes );
229 
230  if ( outFile )
231  {
232  if ( dirOk )
233  gDirectory->cd( ".." );
234  }
235  }
236  else
237  {
238  CreateHisto( histoVector, data.impl, fullName, fullTitle, data.axes );
239  }
240 }
241 
242 
243 void CexmcHistoManager::CreateHisto( CexmcHistoVector & histoVector,
244  CexmcHistoImpl histoImpl, const G4String & name,
245  const G4String & title, const CexmcHistoAxes & axes )
246 {
247  TH1 * histo( NULL );
248 
249  switch ( histoImpl )
250  {
251  case Cexmc_TH1F :
252  histo = new TH1F( name, title, axes.at( 0 ).nBins,
253  axes.at( 0 ).nBinsMin, axes.at( 0 ).nBinsMax );
254  break;
255  case Cexmc_TH2F :
256  histo = new TH2F( name, title, axes.at( 0 ).nBins,
257  axes.at( 0 ).nBinsMin, axes.at( 0 ).nBinsMax,
258  axes.at( 1 ).nBins, axes.at( 1 ).nBinsMin,
259  axes.at( 1 ).nBinsMax );
260  break;
261  case Cexmc_TH3F :
262  histo = new TH3F( name, title, axes.at( 0 ).nBins,
263  axes.at( 0 ).nBinsMin, axes.at( 0 ).nBinsMax,
264  axes.at( 1 ).nBins, axes.at( 1 ).nBinsMin,
265  axes.at( 1 ).nBinsMax, axes.at( 2 ).nBins,
266  axes.at( 2 ).nBinsMin, axes.at( 2 ).nBinsMax );
267  break;
268  default :
269  break;
270  }
271 
272  if ( histo )
273  histoVector.push_back( histo );
274 }
275 
276 
277 void CexmcHistoManager::Initialize( void )
278 {
279  if ( isInitialized )
280  return;
281 
282  CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
283  G4RunManager::GetRunManager() ) );
284  CexmcPhysicsManager * physicsManager( runManager->GetPhysicsManager() );
285 
286  if ( ! physicsManager )
287  throw CexmcException( CexmcWeirdException );
288 
289  CexmcProductionModel * productionModel( physicsManager->
290  GetProductionModel() );
291 
292  if ( ! productionModel )
293  throw CexmcException( CexmcWeirdException );
294 
295  G4ParticleDefinition * outputParticle(
296  productionModel->GetOutputParticle() );
297  G4ParticleDefinition * nucleusOutputParticle(
298  productionModel->GetNucleusOutputParticle() );
299 
300  if ( ! outputParticle || ! nucleusOutputParticle )
301  throw CexmcException( CexmcIncompleteProductionModel );
302 
303  opName = outputParticle->GetParticleName();
304  nopName = nucleusOutputParticle->GetParticleName();
305  opMass = outputParticle->GetPDGMass();
306  nopMass = nucleusOutputParticle->GetPDGMass();
307 
308  G4String title;
309  Int_t nBinsX;
310  Int_t nBinsY;
311  Double_t nBinsMinX;
312  Double_t nBinsMaxX;
313  Double_t nBinsMinY;
314  Double_t nBinsMaxY;
315  CexmcHistoAxes axes;
316 
317  if ( runManager->ProjectIsSaved() )
318  {
319  G4String projectsDir( runManager->GetProjectsDir() );
320  G4String resultsFile( projectsDir + "/" + runManager->GetProjectId() +
321  ".root" );
322  outFile = new TFile( resultsFile, "recreate" );
323  }
324 
325  const CexmcSetup * setup( static_cast< const CexmcSetup * >(
326  runManager->GetUserDetectorConstruction() ) );
327  if ( ! setup )
328  throw CexmcException( CexmcWeirdException );
329 
330  const G4LogicalVolume * lVolume( setup->GetVolume( CexmcSetup::Monitor ) );
331 
332  if ( ! lVolume )
333  throw CexmcException( CexmcIncompatibleGeometry );
334 
335  nBinsMinX = CexmcHistoBeamMomentumMin;
336  nBinsMaxX = CexmcHistoBeamMomentumMax;
337  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) /
338  CexmcHistoBeamMomentumResolution );
339  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
340  AddHisto( CexmcHistoData( CexmcMomentumBP_TPT_Histo, Cexmc_TH1F, false,
341  false, CexmcTPT, "mombp", "Beam momentum at the monitor", axes ) );
342  AddHisto( CexmcHistoData( CexmcMomentumBP_RT_Histo, Cexmc_TH1F, false,
343  false, CexmcRT, "mombp", "Beam momentum at the monitor", axes ) );
344  if ( verboseLevel > 0 )
345  {
346  AddHisto( CexmcHistoData( CexmcMomentumIP_TPT_Histo, Cexmc_TH1F, false,
347  false, CexmcTPT, "momip", "Momentum of the incident particle",
348  axes ) );
349  }
350 
351  G4Box * box( dynamic_cast< G4Box * >( lVolume->GetSolid() ) );
352 
353  if ( ! box )
354  throw CexmcException( CexmcIncompatibleGeometry );
355 
356  G4double width( box->GetXHalfLength() * 2 );
357  G4double height( box->GetYHalfLength() * 2 );
358  G4double halfWidth( width / 2 + CexmcHistoTPSafetyArea );
359  G4double halfHeight( height / 2 + CexmcHistoTPSafetyArea );
360 
361  nBinsX = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
362  nBinsY = Int_t( halfHeight * 2 / CexmcHistoTPResolution );
363  axes.clear();
364  axes.push_back( CexmcHistoAxisData( nBinsX, -halfWidth, halfWidth ) );
365  axes.push_back( CexmcHistoAxisData( nBinsY, -halfHeight, halfHeight ) );
366  AddHisto( CexmcHistoData( CexmcTPInMonitor_TPT_Histo, Cexmc_TH2F, false,
367  false, CexmcTPT, "tpmon", "Track points (mon)", axes ) );
368 
369  lVolume = setup->GetVolume( CexmcSetup::Target );
370  G4Tubs * tube( dynamic_cast< G4Tubs * >( lVolume->GetSolid() ) );
371 
372  if ( ! tube )
373  throw CexmcException( CexmcIncompatibleGeometry );
374 
375  G4double radius( tube->GetOuterRadius() );
376  height = tube->GetZHalfLength() * 2;
377  halfWidth = radius + CexmcHistoTPSafetyArea;
378  halfHeight = height / 2 + CexmcHistoTPSafetyArea;
379 
380  nBinsX = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
381  nBinsY = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
382  Int_t nBinsZ( Int_t( halfHeight * 2 / CexmcHistoTPResolution ) );
383  axes.clear();
384  axes.push_back( CexmcHistoAxisData( nBinsX, -halfWidth, halfWidth ) );
385  axes.push_back( CexmcHistoAxisData( nBinsY, -halfWidth, halfWidth ) );
386  axes.push_back( CexmcHistoAxisData( nBinsZ, -halfHeight, halfHeight ) );
387  AddHisto( CexmcHistoData( CexmcTPInTarget_TPT_Histo, Cexmc_TH3F, false,
388  false, CexmcTPT, "tptar", "Track points (tar)", axes ) );
389  AddHisto( CexmcHistoData( CexmcTPInTarget_RT_Histo, Cexmc_TH3F, false,
390  false, CexmcRT, "tptar", "Track points (tar)", axes ) );
391 
392  title = "Reconstructed masses (" + nopName + " vs. " + opName + ")";
393  nBinsMinX = opMass / 2;
394  nBinsMaxX = opMass + opMass / 2;
395  nBinsMinY = nopMass / 2;
396  nBinsMaxY = nopMass + nopMass / 2;
397  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoMassResolution );
398  nBinsY = Int_t( ( nBinsMaxY - nBinsMinY ) / CexmcHistoMassResolution );
399  axes.clear();
400  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
401  axes.push_back( CexmcHistoAxisData( nBinsY, nBinsMinY, nBinsMaxY ) );
402  AddHisto( CexmcHistoData( CexmcRecMasses_EDT_Histo, Cexmc_TH2F, false,
403  false, CexmcEDT, "recmasses", title, axes ) );
404  AddHisto( CexmcHistoData( CexmcRecMasses_RT_Histo, Cexmc_TH2F, false,
405  false, CexmcRT, "recmasses", title, axes ) );
406 
407  nBinsMinX = 0.;
408  nBinsMaxX = CexmcHistoEnergyMax;
409  nBinsMinY = 0.;
410  nBinsMaxY = CexmcHistoEnergyMax;
411  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoEnergyResolution );
412  nBinsY = Int_t( ( nBinsMaxY - nBinsMinY ) / CexmcHistoEnergyResolution );
413  axes.clear();
414  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
415  axes.push_back( CexmcHistoAxisData( nBinsY, nBinsMinY, nBinsMaxY ) );
416  AddHisto( CexmcHistoData( CexmcAbsorbedEnergy_EDT_Histo, Cexmc_TH2F, false,
417  false, CexmcEDT, "ae", "Absorbed energy (rc vs. lc)", axes ) );
418  AddHisto( CexmcHistoData( CexmcAbsorbedEnergy_RT_Histo, Cexmc_TH2F, false,
419  false, CexmcRT, "ae", "Absorbed energy (rc vs. lc)", axes ) );
420 
421  SetupARHistos( runManager->GetPhysicsManager()->GetProductionModel()->
422  GetAngularRanges() );
423 
424  isInitialized = true;
425 }
426 
427 
428 void CexmcHistoManager::SetupARHistos( const CexmcAngularRangeList & aRanges )
429 {
430  TIter objs( gDirectory->GetList() );
431  TObject * obj( NULL );
432 
433  while ( ( obj = ( TObject * )objs() ) )
434  {
435  TString name( obj->GetName() );
436 
437  if ( outFile )
438  {
439  if ( obj->IsFolder() && ( name.Contains( TString( "_arreal_" ) ) ||
440  name.Contains( TString( "_arrec_" ) ) ) )
441  {
442  gDirectory->cd( name );
443  gDirectory->DeleteAll();
444  gDirectory->cd( ".." );
445  }
446  }
447  else
448  {
449  if ( name.Contains( TRegexp( "_r[0-9]+" ) ) )
450  {
451  gDirectory->Remove( obj );
452  }
453  }
454  }
455 
456  for ( CexmcHistosMap::iterator k( histos.begin() ); k != histos.end();
457  ++k )
458  {
459  if ( k->second.empty() )
460  continue;
461 
462  if ( k->first >= CexmcHistoType_ARReal_START &&
463  k->first <= CexmcHistoType_ARReal_END )
464  {
465  k->second.clear();
466  }
467  }
468 
469  for ( CexmcAngularRangeList::const_iterator k( aRanges.begin() );
470  k != aRanges.end(); ++k )
471  {
472  AddARHistos( *k );
473  }
474 }
475 
476 
477 void CexmcHistoManager::AddARHistos( const CexmcAngularRange & aRange )
478 {
479  G4String title;
480  Int_t nBinsX;
481  Double_t nBinsMinX;
482  Double_t nBinsMaxX;
483  CexmcHistoAxes axes;
484 
485  title = "Reconstructed mass of " + opName;
486  nBinsMinX = opMass / 2;
487  nBinsMaxX = opMass + opMass / 2;
488  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoMassResolution );
489  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
490  AddHisto( CexmcHistoData( CexmcRecMassOP_ARReal_RT_Histo, Cexmc_TH1F, true,
491  false, CexmcRT, "recmassop", title, axes ), aRange );
492 
493  title = "Reconstructed mass of " + nopName;
494  nBinsMinX = nopMass / 2;
495  nBinsMaxX = nopMass + nopMass / 2;
496  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoMassResolution );
497  axes.clear();
498  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
499  AddHisto( CexmcHistoData( CexmcRecMassNOP_ARReal_RT_Histo, Cexmc_TH1F, true,
500  false, CexmcRT, "recmassnop", title, axes ), aRange );
501 
502  G4RunManager * runManager( G4RunManager::GetRunManager() );
503  const CexmcSetup * setup( static_cast< const CexmcSetup * >(
504  runManager->GetUserDetectorConstruction() ) );
505  if ( ! setup )
506  throw CexmcException( CexmcWeirdException );
507 
508  const G4LogicalVolume * lVolume( setup->GetVolume(
509  CexmcSetup::Calorimeter ) );
510 
511  G4Box * box( dynamic_cast< G4Box * >( lVolume->GetSolid() ) );
512 
513  if ( ! box )
514  throw CexmcException( CexmcIncompatibleGeometry );
515 
516  G4double width( box->GetXHalfLength() * 2 );
517  G4double height( box->GetYHalfLength() * 2 );
518  G4double halfWidth( width / 2 + CexmcHistoTPSafetyArea );
519  G4double halfHeight( height / 2 + CexmcHistoTPSafetyArea );
520 
521  nBinsX = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
522  Int_t nBinsY( Int_t( halfHeight * 2 / CexmcHistoTPResolution ) );
523  axes.clear();
524  axes.push_back( CexmcHistoAxisData( nBinsX, -halfWidth, halfWidth ) );
525  axes.push_back( CexmcHistoAxisData( nBinsY, -halfHeight, halfHeight ) );
526 
527  /* looks like there is no possibility to draw descending xaxis in ROOT,
528  * so imagine that you look at calorimeters from behind, i.e. your face to
529  * the beam */
530  AddHisto( CexmcHistoData( CexmcOPDPAtLeftCalorimeter_ARReal_EDT_Histo,
531  Cexmc_TH2F, true, false, CexmcEDT, "opdpcl",
532  "Gamma position on the surface (lc)", axes ), aRange );
533  AddHisto( CexmcHistoData( CexmcOPDPAtRightCalorimeter_ARReal_EDT_Histo,
534  Cexmc_TH2F, true, false, CexmcEDT, "opdpcr",
535  "Gamma position on the surface (rc)", axes ), aRange );
536  AddHisto( CexmcHistoData( CexmcOPDPAtLeftCalorimeter_ARReal_RT_Histo,
537  Cexmc_TH2F, true, false, CexmcRT, "opdpcl",
538  "Gamma position on the surface (lc)", axes ), aRange );
539  AddHisto( CexmcHistoData( CexmcOPDPAtRightCalorimeter_ARReal_RT_Histo,
540  Cexmc_TH2F, true, false, CexmcRT, "opdpcr",
541  "Gamma position on the surface (rc)", axes ), aRange );
542  AddHisto( CexmcHistoData( CexmcRecOPDPAtLeftCalorimeter_ARReal_EDT_Histo,
543  Cexmc_TH2F, true, false, CexmcEDT, "recopdpcl",
544  "Reconstructed gamma position on the surface (lc)", axes ), aRange );
545  AddHisto( CexmcHistoData( CexmcRecOPDPAtRightCalorimeter_ARReal_EDT_Histo,
546  Cexmc_TH2F, true, false, CexmcEDT, "recopdpcr",
547  "Reconstructed gamma position on the surface (rc)", axes ), aRange );
548  AddHisto( CexmcHistoData( CexmcRecOPDPAtLeftCalorimeter_ARReal_RT_Histo,
549  Cexmc_TH2F, true, false, CexmcRT, "recopdpcl",
550  "Reconstructed gamma position on the surface (lc)", axes ), aRange );
551  AddHisto( CexmcHistoData( CexmcRecOPDPAtRightCalorimeter_ARReal_RT_Histo,
552  Cexmc_TH2F, true, false, CexmcRT, "recopdpcr",
553  "Reconstructed gamma position on the surface (rc)", axes ), aRange );
554 
555  nBinsMinX = 0.;
556  nBinsMaxX = CexmcHistoEnergyMax;
557  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoEnergyResolution );
558  axes.clear();
559  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
560  AddHisto( CexmcHistoData( CexmcKinEnAtLeftCalorimeter_ARReal_TPT_Histo,
561  Cexmc_TH1F, true, false, CexmcTPT, "kecl",
562  "Kinetic energy of gamma (lc)", axes ), aRange );
563  AddHisto( CexmcHistoData( CexmcKinEnAtRightCalorimeter_ARReal_TPT_Histo,
564  Cexmc_TH1F, true, false, CexmcTPT, "kecr",
565  "Kinetic energy of gamma (rc)", axes ), aRange );
566  AddHisto( CexmcHistoData( CexmcKinEnAtLeftCalorimeter_ARReal_RT_Histo,
567  Cexmc_TH1F, true, false, CexmcRT, "kecl",
568  "Kinetic energy of gamma (lc)", axes ), aRange );
569  AddHisto( CexmcHistoData( CexmcKinEnAtRightCalorimeter_ARReal_RT_Histo,
570  Cexmc_TH1F, true, false, CexmcRT, "kecr",
571  "Kinetic energy of gamma (rc)", axes ), aRange );
572  AddHisto( CexmcHistoData( CexmcAbsEnInLeftCalorimeter_ARReal_EDT_Histo,
573  Cexmc_TH1F, true, false, CexmcEDT, "aecl",
574  "Absorbed energy (lc)", axes ), aRange );
575  AddHisto( CexmcHistoData( CexmcAbsEnInRightCalorimeter_ARReal_EDT_Histo,
576  Cexmc_TH1F, true, false, CexmcEDT, "aecr",
577  "Absorbed energy (rc)", axes ), aRange );
578  AddHisto( CexmcHistoData( CexmcAbsEnInLeftCalorimeter_ARReal_RT_Histo,
579  Cexmc_TH1F, true, false, CexmcRT, "aecl",
580  "Absorbed energy (lc)", axes ), aRange );
581  AddHisto( CexmcHistoData( CexmcAbsEnInRightCalorimeter_ARReal_RT_Histo,
582  Cexmc_TH1F, true, false, CexmcRT, "aecr",
583  "Absorbed energy (rc)", axes ), aRange );
584 
585  nBinsMinX = CexmcHistoMissEnergyMin;
586  nBinsMaxX = CexmcHistoMissEnergyMax;
587  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) /
588  CexmcHistoMissEnergyResolution );
589  axes.clear();
590  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
591  AddHisto( CexmcHistoData( CexmcMissEnFromLeftCalorimeter_ARReal_RT_Histo,
592  Cexmc_TH1F, true, false, CexmcRT, "mecl",
593  "Missing energy (lc)", axes ), aRange );
594  AddHisto( CexmcHistoData( CexmcMissEnFromRightCalorimeter_ARReal_RT_Histo,
595  Cexmc_TH1F, true, false, CexmcRT, "mecr",
596  "Missing energy (rc)", axes ), aRange );
597 
598  title = "Kinetic energy of newborn " + opName + " (lab)";
599  nBinsMinX = 0.;
600  nBinsMaxX = CexmcHistoEnergyMax;
601  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoEnergyResolution );
602  axes.clear();
603  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
604  AddHisto( CexmcHistoData( CexmcKinEnOP_LAB_ARReal_TPT_Histo,
605  Cexmc_TH1F, true, false, CexmcTPT, "keop_lab", title, axes ), aRange );
606  AddHisto( CexmcHistoData( CexmcKinEnOP_LAB_ARReal_RT_Histo,
607  Cexmc_TH1F, true, false, CexmcRT, "keop_lab", title, axes ), aRange );
608 
609  title = "Angle of newborn " + opName + " (scm)";
610  nBinsMinX = -1.0;
611  nBinsMaxX = 1.0;
612  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoAngularCResolution );
613  axes.clear();
614  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
615  AddHisto( CexmcHistoData( CexmcAngleOP_SCM_ARReal_TPT_Histo,
616  Cexmc_TH1F, true, false, CexmcTPT, "aop_scm", title, axes ), aRange );
617  AddHisto( CexmcHistoData( CexmcAngleOP_SCM_ARReal_RT_Histo,
618  Cexmc_TH1F, true, false, CexmcRT, "aop_scm", title, axes ), aRange );
619 
620  title = "Reconstruced angle of newborn " + opName + " (scm)";
621  AddHisto( CexmcHistoData( CexmcRecAngleOP_SCM_ARReal_RT_Histo,
622  Cexmc_TH1F, true, false, CexmcRT, "recaop_scm", title, axes ), aRange );
623 
624  title = "Real - reconstruced angle of newborn " + opName + " (scm)";
625  AddHisto( CexmcHistoData( CexmcDiffAngleOP_SCM_ARReal_RT_Histo,
626  Cexmc_TH1F, true, false, CexmcRT, "diffaop_scm", title, axes ),
627  aRange );
628 
629  nBinsMinX = 0.;
630  nBinsMaxX = 360.;
631  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoAngularResolution );
632  axes.clear();
633  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
634  AddHisto( CexmcHistoData( CexmcOpenAngle_ARReal_TPT_Histo,
635  Cexmc_TH1F, true, false, CexmcTPT, "oa",
636  "Open angle between the gammas", axes ), aRange );
637  AddHisto( CexmcHistoData( CexmcOpenAngle_ARReal_RT_Histo,
638  Cexmc_TH1F, true, false, CexmcRT, "oa",
639  "Open angle between the gammas", axes ), aRange );
640  AddHisto( CexmcHistoData( CexmcRecOpenAngle_ARReal_RT_Histo,
641  Cexmc_TH1F, true, false, CexmcRT, "recoa",
642  "Reconstructed open angle between the gammas", axes ), aRange );
643 
644  nBinsMinX = -180.;
645  nBinsMaxX = 180.;
646  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoAngularResolution );
647  axes.clear();
648  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
649  AddHisto( CexmcHistoData( CexmcDiffOpenAngle_ARReal_RT_Histo,
650  Cexmc_TH1F, true, false, CexmcRT, "diffoa",
651  "Real - reconstructed open angle between the gammas", axes ), aRange );
652 
653  lVolume = setup->GetVolume( CexmcSetup::Target );
654  G4Tubs * tube( dynamic_cast< G4Tubs * >( lVolume->GetSolid() ) );
655 
656  if ( ! tube )
657  throw CexmcException( CexmcIncompatibleGeometry );
658 
659  G4double radius( tube->GetOuterRadius() );
660  height = tube->GetZHalfLength() * 2;
661  halfWidth = radius + CexmcHistoTPSafetyArea;
662  halfHeight = height / 2 + CexmcHistoTPSafetyArea;
663 
664  nBinsX = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
665  nBinsY = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
666  Int_t nBinsZ( Int_t( halfHeight * 2 / CexmcHistoTPResolution ) );
667  axes.clear();
668  axes.push_back( CexmcHistoAxisData( nBinsX, -halfWidth, halfWidth ) );
669  axes.push_back( CexmcHistoAxisData( nBinsY, -halfWidth, halfWidth ) );
670  axes.push_back( CexmcHistoAxisData( nBinsZ, -halfHeight, halfHeight ) );
671  AddHisto( CexmcHistoData( CexmcTPInTarget_ARReal_TPT_Histo, Cexmc_TH3F,
672  true, false, CexmcTPT, "tptar", "Track points (tar)", axes ), aRange );
673  AddHisto( CexmcHistoData( CexmcTPInTarget_ARReal_RT_Histo, Cexmc_TH3F,
674  true, false, CexmcRT, "tptar", "Track points (tar)", axes ), aRange );
675 }
676 
677 
678 void CexmcHistoManager::Add( CexmcHistoType histoType, unsigned int index,
679  G4double x )
680 {
681  CexmcHistosMap::iterator found( histos.find( histoType ) );
682  if ( found == histos.end() || histos[ histoType ].size() <= index )
683  throw CexmcException( CexmcWeirdException );
684 
685  histos[ histoType ][ index ]->Fill( x );
686 }
687 
688 
689 void CexmcHistoManager::Add( CexmcHistoType histoType, unsigned int index,
690  G4double x, G4double y )
691 {
692  CexmcHistosMap::iterator found( histos.find( histoType ) );
693  if ( found == histos.end() || histos[ histoType ].size() <= index )
694  throw CexmcException( CexmcWeirdException );
695 
696  /* no cast needed because TH1 has virtual method
697  * Fill( Double_t, Double_t ) */
698  histos[ histoType ][ index ]->Fill( x, y );
699 }
700 
701 
702 void CexmcHistoManager::Add( CexmcHistoType histoType, unsigned int index,
703  G4double x, G4double y, G4double z )
704 {
705  CexmcHistosMap::iterator found( histos.find( histoType ) );
706  if ( found == histos.end() || histos[ histoType ].size() <= index )
707  throw CexmcException( CexmcWeirdException );
708 
709  /* cast needed because TH1 does not have virtual method
710  * Fill( Double_t, Double_t, Double_t ) */
711  TH3 * histo( static_cast< TH3 * >( histos[ histoType ][ index ] ) );
712 
713  histo->Fill( x, y, z );
714 }
715 
716 
717 void CexmcHistoManager::Add( CexmcHistoType histoType, unsigned int index,
718  G4int binX, G4int binY, G4double value )
719 {
720  CexmcHistosMap::iterator found( histos.find( histoType ) );
721  if ( found == histos.end() || histos[ histoType ].size() <= index )
722  throw CexmcException( CexmcWeirdException );
723 
724  ++binX;
725  ++binY;
726  Double_t curValue( histos[ histoType ][ index ]->GetBinContent(
727  binX, binY ) );
728  histos[ histoType ][ index ]->SetBinContent( binX, binY,
729  curValue + value / GeV );
730 }
731 
732 
733 void CexmcHistoManager::List( void ) const
734 {
735  gDirectory->ls();
736 }
737 
738 
739 void CexmcHistoManager::Print( const G4String & value )
740 {
741  TObject * histo( gDirectory->FindObjectAny( value.c_str() ) );
742 
743  if ( ! histo )
744  {
745  G4cout << "Histogram '" << value << "' was not found" << G4endl;
746  return;
747  }
748 
749  histo->Print( "range" );
750 }
751 
752 
753 #ifdef CEXMC_USE_ROOTQT
754 
755 void CexmcHistoManager::Draw( const G4String & histoName,
756  const G4String & histoDrawOptions )
757 {
758  CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
759  G4RunManager::GetRunManager() ) );
760  if ( ! runManager->AreLiveHistogramsEnabled() )
761  {
762  G4cout << "Live histograms option is disabled" << G4endl;
763  return;
764  }
765 
766  TObject * histo( gDirectory->FindObjectAny( histoName ) );
767 
768  if ( ! histo )
769  {
770  G4cout << "Histogram '" << histoName << "' was not found" << G4endl;
771  return;
772  }
773 
774  if ( ! rootCanvas )
775  {
776  /* save default application font because rootCanvas will break it */
777  QFont defaultAppFont( QApplication::font() );
778  rootCanvas = new CexmcHistoWidget;
779  QApplication::setFont( defaultAppFont );
780  rootCanvas->resize( CexmcHistoCanvasWidth, CexmcHistoCanvasHeight );
781  rootCanvas->GetCanvas()->cd();
782  }
783 
784  histo->Draw( histoDrawOptions );
785  rootCanvas->show();
786  rootCanvas->GetCanvas()->Update();
787 }
788 
789 #endif
790 
791 #endif
792 
CexmcHistoManager.hh
UBox::GetXHalfLength
double GetXHalfLength() const
Definition: UBox.hh:116
index
Int_t index
Definition: macro.C:9
outFile
std::ofstream outFile
Definition: GammaRayTel.cc:68
Draw
hist1 Draw("HIST")
G4Box
Definition: G4Box.hh:63
Initialize
void Initialize()
Definition: errprop.cc:96
UBox::GetYHalfLength
double GetYHalfLength() const
Definition: UBox.hh:120
CexmcEDT
Definition: CexmcCommon.hh:101
G4Tubs
Definition: G4Tubs.hh:84
CexmcIncompatibleGeometry
Definition: CexmcException.hh:72
width
#define width
G4Box.hh
name
const XML_Char * name
Definition: expat.h:151
test.x
tuple x
Definition: test.py:50
CexmcProductionModel.hh
CexmcIncompleteProductionModel
Definition: CexmcException.hh:73
G4int
int G4int
Definition: G4Types.hh:78
y
Double_t y
Definition: plot.C:279
G4cout
G4GLOB_DLL std::ostream G4cout
G4INCL::Int_t
G4int Int_t
Definition: G4INCLEventInfo.hh:59
G4Tubs.hh
G4bool
bool G4bool
Definition: G4Types.hh:79
CexmcPhysicsManager.hh
G4RunManager::GetRunManager
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:74
G4LogicalVolume.hh
title
subroutine title(NA, NB, NCA, NCB)
Definition: dpm25nuc7.f:1744
CexmcException.hh
test.z
tuple z
Definition: test.py:28
CexmcAngularRangeList
std::vector< CexmcAngularRange > CexmcAngularRangeList
Definition: CexmcAngularRange.hh:72
value
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4endl
#define G4endl
Definition: G4ios.hh:61
CexmcRunManager.hh
CexmcRT
Definition: CexmcCommon.hh:102
G4INCL::Double_t
G4double Double_t
Definition: G4INCLEventInfo.hh:62
G4double
double G4double
Definition: G4Types.hh:76
box
UBox box
Definition: UPolycone.cc:32
CexmcHistoManagerMessenger.hh
CexmcHistoWidget.hh
CexmcTPT
Definition: CexmcCommon.hh:100
CexmcWeirdException
Definition: CexmcException.hh:88
G4INCL::Random::isInitialized
G4bool isInitialized()
Definition: G4INCLRandom.cc:158
data
const XML_Char const XML_Char * data
Definition: expat.h:268
G4ArrayOps::Add
void Add(G4int Elements, T *To, T *A1, T *A2=NULL)
Definition: G4ArrayOps.hh:77
CexmcSetup.hh
pyG4Element::Print
void Print(G4Element &ele)
Definition: pyG4Element.cc:56
Generated on Sat Dec 14 2013 14:32:55 for Geant4_10 by   doxygen 1.8.5