Geant4  10.02.p02
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 #ifdef CEXMC_USE_ROOTQT
58 #include <TCanvas.h>
59 #include <TList.h>
60 #include <QApplication>
61 #include <QFont>
62 #include <QMenu>
63 #include <G4UIQt.hh>
64 #include "CexmcMessenger.hh"
65 #endif
66 #include <G4LogicalVolume.hh>
67 #include <G4Box.hh>
68 #include <G4Tubs.hh>
69 #include <G4SystemOfUnits.hh>
70 #include "CexmcHistoManager.hh"
71 #include "CexmcHistoManagerMessenger.hh"
72 #include "CexmcProductionModel.hh"
73 #include "CexmcPhysicsManager.hh"
74 #include "CexmcRunManager.hh"
75 #include "CexmcSetup.hh"
76 #include "CexmcException.hh"
77 #include "CexmcHistoWidget.hh"
78 
79 
80 namespace
81 {
82  const G4double CexmcHistoBeamMomentumMin( 0.0 * GeV );
83  const G4double CexmcHistoBeamMomentumMax( 1.0 * GeV );
84  const G4double CexmcHistoBeamMomentumResolution( 0.5 * MeV );
85  const G4double CexmcHistoTPResolution( 0.1 * cm );
86  const G4double CexmcHistoTPSafetyArea( 1.0 * cm );
87  const G4double CexmcHistoMassResolution( 1.0 * MeV );
88  const G4double CexmcHistoEnergyMax( 1.0 * GeV );
89  const G4double CexmcHistoEnergyResolution( 1.0 * MeV );
90  const G4double CexmcHistoMissEnergyMin( -0.1 * GeV );
91  const G4double CexmcHistoMissEnergyMax( 0.2 * GeV );
92  const G4double CexmcHistoMissEnergyResolution( 0.2 * MeV );
93  const G4double CexmcHistoAngularResolution( 0.5 );
94  const G4double CexmcHistoAngularCResolution( 0.001 );
95  const G4int CexmcHistoCanvasWidth( 800 );
96  const G4int CexmcHistoCanvasHeight( 600 );
97  const G4String CexmcHistoDirectoryHandle( "histograms" );
98  const G4String CexmcHistoDirectoryTitle( "Histograms" );
99 }
100 
101 
102 CexmcHistoManager * CexmcHistoManager::instance( NULL );
103 
104 
105 CexmcHistoManager * CexmcHistoManager::Instance( void )
106 {
107  if ( instance == NULL )
108  instance = new CexmcHistoManager;
109 
110  return instance;
111 }
112 
113 
114 void CexmcHistoManager::Destroy( void )
115 {
116  delete instance;
117  instance = NULL;
118 }
119 
120 
121 CexmcHistoManager::CexmcHistoManager() : outFile( NULL ),
122  isInitialized( false ), opName( "" ), nopName( "" ), opMass( 0. ),
123  nopMass( 0. ), verboseLevel( 0 ),
124 #ifdef CEXMC_USE_ROOTQT
125  rootCanvas( NULL ), areLiveHistogramsEnabled( false ),
126  isHistoMenuInitialized( false ), drawOptions1D( "" ), drawOptions2D( "" ),
127  drawOptions3D( "" ), histoMenuHandle( "" ), histoMenuLabel( "" ),
128 #endif
129  messenger( NULL )
130 {
131  for ( int i( 0 ); i < CexmcHistoType_SIZE; ++i )
132  {
133  histos.insert( CexmcHistoPair( CexmcHistoType( i ),
134  CexmcHistoVector() ) );
135  }
136 
137  messenger = new CexmcHistoManagerMessenger( this );
138 }
139 
140 
141 CexmcHistoManager::~CexmcHistoManager()
142 {
143  if ( outFile )
144  {
145  outFile->Write();
146  outFile->Close();
147  }
148 
149  /* all histograms will be deleted by outFile destructor! */
150  delete outFile;
151 #ifdef CEXMC_USE_ROOTQT
152  delete rootCanvas;
153 #endif
154  delete messenger;
155 }
156 
157 
158 void CexmcHistoManager::AddHisto( const CexmcHistoData & data,
159  const CexmcAngularRange & aRange )
160 {
161  G4String fullName( data.name );
162  G4String fullTitle( data.title );
163  G4String rangeTypeLabel;
164  G4String triggerTypeLabel;
165  G4String decorTriggerTypeLabel;
166 
167  if ( data.isARHisto )
168  {
169  if ( data.isARRec )
170  {
171  fullName += "_arrec";
172  rangeTypeLabel = "rec";
173  }
174  else
175  {
176  fullName += "_arreal";
177  rangeTypeLabel = "real";
178  }
179  }
180 
181  switch ( data.triggerType )
182  {
183  case CexmcTPT :
184  fullName += "_tpt";
185  decorTriggerTypeLabel = " --tpt--";
186  fullTitle += decorTriggerTypeLabel;
187  triggerTypeLabel = "tpt";
188  break;
189  case CexmcEDT :
190  fullName += "_edt";
191  decorTriggerTypeLabel = " --edt--";
192  fullTitle += decorTriggerTypeLabel;
193  triggerTypeLabel = "edt";
194  break;
195  case CexmcRT :
196  fullName += "_rt";
197  decorTriggerTypeLabel = " --rt--";
198  fullTitle += decorTriggerTypeLabel;
199  triggerTypeLabel = "rt";
200  break;
201  default :
202  break;
203  }
204 
205  CexmcHistosMap::iterator found( histos.find( data.type ) );
206 
207  if ( found == histos.end() )
208  throw CexmcException( CexmcWeirdException );
209 
210  CexmcHistoVector & histoVector( found->second );
211 
212  if ( data.isARHisto )
213  {
214  G4bool dirOk( false );
215 
216  dirOk = gDirectory->Get( fullName ) != NULL;
217 
218  if ( ! dirOk )
219  dirOk = ( gDirectory->mkdir( fullName, fullTitle ) != NULL );
220 
221  if ( dirOk )
222  gDirectory->cd( fullName );
223 
224  std::ostringstream histoName;
225  std::ostringstream histoTitle;
226  histoName << data.name << "_r" << aRange.index + 1 << rangeTypeLabel <<
227  "_" << triggerTypeLabel;
228  histoTitle << data.title << " {range " << aRange.index + 1 <<
229  rangeTypeLabel << " [" << std::fixed <<
230  std::setprecision( 4 ) << aRange.top << ", " <<
231  aRange.bottom << ")}" << decorTriggerTypeLabel;
232  CreateHisto( histoVector, data.impl, histoName.str(), histoTitle.str(),
233  data.axes );
234 
235  if ( dirOk )
236  gDirectory->cd( ".." );
237  }
238  else
239  {
240  CreateHisto( histoVector, data.impl, fullName, fullTitle, data.axes );
241  }
242 }
243 
244 
245 void CexmcHistoManager::CreateHisto( CexmcHistoVector & histoVector,
246  CexmcHistoImpl histoImpl, const G4String & name,
247  const G4String & title, const CexmcHistoAxes & axes )
248 {
249  TH1 * histo( NULL );
250 
251  switch ( histoImpl )
252  {
253  case Cexmc_TH1F :
254  histo = new TH1F( name, title, axes.at( 0 ).nBins,
255  axes.at( 0 ).nBinsMin, axes.at( 0 ).nBinsMax );
256  break;
257  case Cexmc_TH2F :
258  histo = new TH2F( name, title, axes.at( 0 ).nBins,
259  axes.at( 0 ).nBinsMin, axes.at( 0 ).nBinsMax,
260  axes.at( 1 ).nBins, axes.at( 1 ).nBinsMin,
261  axes.at( 1 ).nBinsMax );
262  break;
263  case Cexmc_TH3F :
264  histo = new TH3F( name, title, axes.at( 0 ).nBins,
265  axes.at( 0 ).nBinsMin, axes.at( 0 ).nBinsMax,
266  axes.at( 1 ).nBins, axes.at( 1 ).nBinsMin,
267  axes.at( 1 ).nBinsMax, axes.at( 2 ).nBins,
268  axes.at( 2 ).nBinsMin, axes.at( 2 ).nBinsMax );
269  break;
270  default :
271  break;
272  }
273 
274  if ( histo )
275  histoVector.push_back( histo );
276 }
277 
278 
279 void CexmcHistoManager::Initialize( void )
280 {
281  if ( isInitialized )
282  return;
283 
284  CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
285  G4RunManager::GetRunManager() ) );
286  CexmcPhysicsManager * physicsManager( runManager->GetPhysicsManager() );
287 
288  if ( ! physicsManager )
289  throw CexmcException( CexmcWeirdException );
290 
291  CexmcProductionModel * productionModel( physicsManager->
292  GetProductionModel() );
293 
294  if ( ! productionModel )
295  throw CexmcException( CexmcWeirdException );
296 
297  G4ParticleDefinition * outputParticle(
298  productionModel->GetOutputParticle() );
299  G4ParticleDefinition * nucleusOutputParticle(
300  productionModel->GetNucleusOutputParticle() );
301 
302  if ( ! outputParticle || ! nucleusOutputParticle )
303  throw CexmcException( CexmcIncompleteProductionModel );
304 
305  opName = outputParticle->GetParticleName();
306  nopName = nucleusOutputParticle->GetParticleName();
307  opMass = outputParticle->GetPDGMass();
308  nopMass = nucleusOutputParticle->GetPDGMass();
309 
310  G4String title;
311  Int_t nBinsX;
312  Int_t nBinsY;
313  Double_t nBinsMinX;
314  Double_t nBinsMaxX;
315  Double_t nBinsMinY;
316  Double_t nBinsMaxY;
317  CexmcHistoAxes axes;
318 
319  if ( runManager->ProjectIsSaved() )
320  {
321  G4String projectsDir( runManager->GetProjectsDir() );
322  G4String resultsFile( projectsDir + "/" + runManager->GetProjectId() +
323  ".root" );
324  outFile = new TFile( resultsFile, "recreate" );
325  }
326  else
327  {
328  outFile = new TDirectoryFile( CexmcHistoDirectoryHandle,
329  CexmcHistoDirectoryTitle );
330  gDirectory->cd( CexmcHistoDirectoryHandle );
331  }
332 
333  if ( ! outFile )
334  throw CexmcException( CexmcWeirdException );
335 
336  const CexmcSetup * setup( static_cast< const CexmcSetup * >(
337  runManager->GetUserDetectorConstruction() ) );
338  if ( ! setup )
339  throw CexmcException( CexmcWeirdException );
340 
341  const G4LogicalVolume * lVolume( setup->GetVolume( CexmcSetup::Monitor ) );
342 
343  if ( ! lVolume )
344  throw CexmcException( CexmcIncompatibleGeometry );
345 
346  nBinsMinX = CexmcHistoBeamMomentumMin;
347  nBinsMaxX = CexmcHistoBeamMomentumMax;
348  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) /
349  CexmcHistoBeamMomentumResolution );
350  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
351  AddHisto( CexmcHistoData( CexmcMomentumBP_TPT_Histo, Cexmc_TH1F, false,
352  false, CexmcTPT, "mombp", "Beam momentum at the monitor", axes ) );
353  AddHisto( CexmcHistoData( CexmcMomentumBP_RT_Histo, Cexmc_TH1F, false,
354  false, CexmcRT, "mombp", "Beam momentum at the monitor", axes ) );
355  if ( verboseLevel > 0 )
356  {
357  AddHisto( CexmcHistoData( CexmcMomentumIP_TPT_Histo, Cexmc_TH1F, false,
358  false, CexmcTPT, "momip", "Momentum of the incident particle",
359  axes ) );
360  }
361 
362  G4Box * box( dynamic_cast< G4Box * >( lVolume->GetSolid() ) );
363 
364  if ( ! box )
365  throw CexmcException( CexmcIncompatibleGeometry );
366 
367  G4double width( box->GetXHalfLength() * 2 );
368  G4double height( box->GetYHalfLength() * 2 );
369  G4double halfWidth( width / 2 + CexmcHistoTPSafetyArea );
370  G4double halfHeight( height / 2 + CexmcHistoTPSafetyArea );
371 
372  nBinsX = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
373  nBinsY = Int_t( halfHeight * 2 / CexmcHistoTPResolution );
374  axes.clear();
375  axes.push_back( CexmcHistoAxisData( nBinsX, -halfWidth, halfWidth ) );
376  axes.push_back( CexmcHistoAxisData( nBinsY, -halfHeight, halfHeight ) );
377  AddHisto( CexmcHistoData( CexmcTPInMonitor_TPT_Histo, Cexmc_TH2F, false,
378  false, CexmcTPT, "tpmon", "Track points (mon)", axes ) );
379 
380  lVolume = setup->GetVolume( CexmcSetup::Target );
381  G4Tubs * tube( dynamic_cast< G4Tubs * >( lVolume->GetSolid() ) );
382 
383  if ( ! tube )
384  throw CexmcException( CexmcIncompatibleGeometry );
385 
386  G4double radius( tube->GetOuterRadius() );
387  height = tube->GetZHalfLength() * 2;
388  halfWidth = radius + CexmcHistoTPSafetyArea;
389  halfHeight = height / 2 + CexmcHistoTPSafetyArea;
390 
391  nBinsX = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
392  nBinsY = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
393  Int_t nBinsZ( Int_t( halfHeight * 2 / CexmcHistoTPResolution ) );
394  axes.clear();
395  axes.push_back( CexmcHistoAxisData( nBinsX, -halfWidth, halfWidth ) );
396  axes.push_back( CexmcHistoAxisData( nBinsY, -halfWidth, halfWidth ) );
397  axes.push_back( CexmcHistoAxisData( nBinsZ, -halfHeight, halfHeight ) );
398  AddHisto( CexmcHistoData( CexmcTPInTarget_TPT_Histo, Cexmc_TH3F, false,
399  false, CexmcTPT, "tptar", "Track points (tar)", axes ) );
400  AddHisto( CexmcHistoData( CexmcTPInTarget_RT_Histo, Cexmc_TH3F, false,
401  false, CexmcRT, "tptar", "Track points (tar)", axes ) );
402 
403  title = "Reconstructed masses (" + nopName + " vs. " + opName + ")";
404  nBinsMinX = opMass / 2;
405  nBinsMaxX = opMass + opMass / 2;
406  nBinsMinY = nopMass / 2;
407  nBinsMaxY = nopMass + nopMass / 2;
408  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoMassResolution );
409  nBinsY = Int_t( ( nBinsMaxY - nBinsMinY ) / CexmcHistoMassResolution );
410  axes.clear();
411  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
412  axes.push_back( CexmcHistoAxisData( nBinsY, nBinsMinY, nBinsMaxY ) );
413  AddHisto( CexmcHistoData( CexmcRecMasses_EDT_Histo, Cexmc_TH2F, false,
414  false, CexmcEDT, "recmasses", title, axes ) );
415  AddHisto( CexmcHistoData( CexmcRecMasses_RT_Histo, Cexmc_TH2F, false,
416  false, CexmcRT, "recmasses", title, axes ) );
417 
418  nBinsMinX = 0.;
419  nBinsMaxX = CexmcHistoEnergyMax;
420  nBinsMinY = 0.;
421  nBinsMaxY = CexmcHistoEnergyMax;
422  nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoEnergyResolution );
423  nBinsY = Int_t( ( nBinsMaxY - nBinsMinY ) / CexmcHistoEnergyResolution );
424  axes.clear();
425  axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
426  axes.push_back( CexmcHistoAxisData( nBinsY, nBinsMinY, nBinsMaxY ) );
427  AddHisto( CexmcHistoData( CexmcAbsorbedEnergy_EDT_Histo, Cexmc_TH2F, false,
428  false, CexmcEDT, "ae", "Absorbed energy (rc vs. lc)", axes ) );
429  AddHisto( CexmcHistoData( CexmcAbsorbedEnergy_RT_Histo, Cexmc_TH2F, false,
430  false, CexmcRT, "ae", "Absorbed energy (rc vs. lc)", axes ) );
431 
432  SetupARHistos( runManager->GetPhysicsManager()->GetProductionModel()->
433  GetAngularRanges() );
434 
435  isInitialized = true;
436 }
437 
438 
439 void CexmcHistoManager::SetupARHistos( const CexmcAngularRangeList & aRanges )
440 {
441  TIter objs( gDirectory->GetList() );
442  TObject * obj( NULL );
443 
444  while ( ( obj = ( TObject * )objs() ) )
445  {
446  TString name( obj->GetName() );
447 
448  if ( obj->IsFolder() && ( name.Contains( TString( "_arreal_" ) ) ||
449  name.Contains( TString( "_arrec_" ) ) ) )
450  {
451  gDirectory->cd( name );
452  gDirectory->DeleteAll();
453  gDirectory->cd( ".." );
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  /* BEWARE: list will be printed on stdout */
737  gDirectory->ls();
738 }
739 
740 
741 void CexmcHistoManager::Print( const G4String & value )
742 {
743  TObject * histo( gDirectory->FindObjectAny( value.c_str() ) );
744 
745  if ( ! histo )
746  {
747  G4cout << "Histogram '" << value << "' was not found" << G4endl;
748  return;
749  }
750 
751  /* BEWARE: histo will be printed on stdout */
752  histo->Print( "range" );
753 }
754 
755 
756 #ifdef CEXMC_USE_ROOTQT
757 
758 void CexmcHistoManager::Draw( const G4String & histoName,
759  const G4String & histoDrawOptions )
760 {
761  if ( ! 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 
791 void CexmcHistoManager::EnableLiveHistograms( G4UIsession * session,
792  G4bool on )
793 {
794  areLiveHistogramsEnabled = on;
795 
796  if ( ! on || ! isInitialized )
797  return;
798 
799  G4UIQt * qtSession( dynamic_cast< G4UIQt * >( session ) );
800 
801  if ( ! qtSession )
802  return;
803 
804  if ( ! histoMenuHandle.empty() && ! isHistoMenuInitialized )
805  {
806  qtSession->AddMenu( histoMenuHandle, histoMenuLabel );
807  BuildMenuTree( qtSession, histoMenuHandle, gDirectory->GetList() );
808  isHistoMenuInitialized = true;
809  }
810 }
811 
812 
813 void CexmcHistoManager::BuildMenuTree( G4UIQt * session,
814  const G4String & menu, TList * ls )
815 {
816  TIter objs( ls );
817  TObject * obj( NULL );
818 
819  while ( ( obj = ( TObject * )objs() ) )
820  {
821  G4String name( obj->GetName() );
822  G4String title( obj->GetTitle() );
823 
824  if ( obj->IsFolder() )
825  {
826  AddSubmenu( session, menu, name, title );
827  BuildMenuTree( session, name, ( ( TDirectory * )obj )->GetList() );
828  }
829  else
830  {
831  G4String options( name );
832 
833  do
834  {
835  if ( obj->InheritsFrom( TH3::Class() ) &&
836  ! drawOptions3D.empty() )
837  {
838  title = G4String( "3: " ) + title;
839  options += G4String( " " ) + drawOptions3D;
840  break;
841  }
842  if ( obj->InheritsFrom( TH2::Class() ) &&
843  ! drawOptions2D.empty() )
844  {
845  title = G4String( "2: " ) + title;
846  options += G4String( " " ) + drawOptions2D;
847  break;
848  }
849  if ( obj->InheritsFrom( TH1::Class() ) &&
850  ! drawOptions1D.empty() )
851  {
852  options += G4String( " " ) + drawOptions1D;
853  break;
854  }
855  } while ( false );
856 
857  G4String cmd( CexmcMessenger::histoDirName + "draw " + options );
858  session->AddButton( menu, title.c_str(), cmd );
859  }
860  }
861 }
862 
863 
864 void CexmcHistoManager::AddSubmenu( G4UIQt * session,
865  const G4String & parent,
866  const G4String & name,
867  const G4String & label )
868 {
869  QMenu * menu( new QMenu( label.c_str() ) );
870  QMenu * parentMenu( ( QMenu * )session->GetInteractor( parent ) );
871 
872  parentMenu->addMenu( menu );
873  session->AddInteractor( name, ( G4Interactor )menu );
874 }
875 
876 #endif
877 
878 #endif
879 
cm
static const double cm
Definition: G4SIunits.hh:118
CexmcHistoManager.hh
MeV
static const double MeV
Definition: G4SIunits.hh:211
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
G4UIQt.hh
CexmcIncompatibleGeometry
Definition: CexmcException.hh:72
width
#define width
G4Box.hh
CexmcProductionModel.hh
CexmcIncompleteProductionModel
Definition: CexmcException.hh:73
G4int
int G4int
Definition: G4Types.hh:78
G4Interactor
void * G4Interactor
Definition: G4VInteractorManager.hh:37
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:214
Print
void Print(const std::vector< T > &data)
Definition: DicomRun.hh:109
session
static G4UIterminal * session
Definition: pyG4UIterminal.cc:38
CexmcMessenger.hh
G4RunManager::GetRunManager
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
G4LogicalVolume.hh
CexmcException.hh
x
const G4double x[NPOINTSGL]
Definition: G4PreCompoundFragment.cc:56
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