| Geant4 Cross Reference |
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
882