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