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