File indexing completed on 2025-01-31 09:21:53
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044 #ifdef CEXMC_USE_ROOT
0045
0046 #include <iostream>
0047 #include <iomanip>
0048 #include <TH1.h>
0049 #include <TH1F.h>
0050 #include <TH2F.h>
0051 #include <TH3F.h>
0052 #include <TFile.h>
0053 #include <TObject.h>
0054 #include <TCollection.h>
0055 #include <TDirectory.h>
0056 #include <TString.h>
0057 #ifdef CEXMC_USE_ROOTQT
0058 #include <TCanvas.h>
0059 #include <TList.h>
0060 #include <QApplication>
0061 #include <QFont>
0062 #include <QMenu>
0063 #include <G4UIQt.hh>
0064 #include "CexmcMessenger.hh"
0065 #endif
0066 #include <G4LogicalVolume.hh>
0067 #include <G4Box.hh>
0068 #include <G4Tubs.hh>
0069 #include <G4SystemOfUnits.hh>
0070 #include "CexmcHistoManager.hh"
0071 #include "CexmcHistoManagerMessenger.hh"
0072 #include "CexmcProductionModel.hh"
0073 #include "CexmcPhysicsManager.hh"
0074 #include "CexmcRunManager.hh"
0075 #include "CexmcSetup.hh"
0076 #include "CexmcException.hh"
0077 #include "CexmcHistoWidget.hh"
0078
0079
0080 namespace
0081 {
0082 const G4double CexmcHistoBeamMomentumMin( 0.0 * GeV );
0083 const G4double CexmcHistoBeamMomentumMax( 1.0 * GeV );
0084 const G4double CexmcHistoBeamMomentumResolution( 0.5 * MeV );
0085 const G4double CexmcHistoTPResolution( 0.1 * cm );
0086 const G4double CexmcHistoTPSafetyArea( 1.0 * cm );
0087 const G4double CexmcHistoMassResolution( 1.0 * MeV );
0088 const G4double CexmcHistoEnergyMax( 1.0 * GeV );
0089 const G4double CexmcHistoEnergyResolution( 1.0 * MeV );
0090 const G4double CexmcHistoMissEnergyMin( -0.1 * GeV );
0091 const G4double CexmcHistoMissEnergyMax( 0.2 * GeV );
0092 const G4double CexmcHistoMissEnergyResolution( 0.2 * MeV );
0093 const G4double CexmcHistoAngularResolution( 0.5 );
0094 const G4double CexmcHistoAngularCResolution( 0.001 );
0095 #ifdef CEXMC_USE_ROOTQT
0096 const G4int CexmcHistoCanvasWidth( 800 );
0097 const G4int CexmcHistoCanvasHeight( 600 );
0098 #endif
0099 const G4String CexmcHistoDirectoryHandle( "histograms" );
0100 const G4String CexmcHistoDirectoryTitle( "Histograms" );
0101 }
0102
0103
0104 CexmcHistoManager * CexmcHistoManager::instance( NULL );
0105
0106
0107 CexmcHistoManager * CexmcHistoManager::Instance( void )
0108 {
0109 if ( instance == NULL )
0110 instance = new CexmcHistoManager;
0111
0112 return instance;
0113 }
0114
0115
0116 void CexmcHistoManager::Destroy( void )
0117 {
0118 delete instance;
0119 instance = NULL;
0120 }
0121
0122
0123 CexmcHistoManager::CexmcHistoManager() : outFile( NULL ),
0124 isInitialized( false ), opName( "" ), nopName( "" ), opMass( 0. ),
0125 nopMass( 0. ), verboseLevel( 0 ),
0126 #ifdef CEXMC_USE_ROOTQT
0127 rootCanvas( NULL ), areLiveHistogramsEnabled( false ),
0128 isHistoMenuInitialized( false ), drawOptions1D( "" ), drawOptions2D( "" ),
0129 drawOptions3D( "" ), histoMenuHandle( "" ), histoMenuLabel( "" ),
0130 #endif
0131 messenger( NULL )
0132 {
0133 for ( int i( 0 ); i < CexmcHistoType_SIZE; ++i )
0134 {
0135 histos.insert( CexmcHistoPair( CexmcHistoType( i ),
0136 CexmcHistoVector() ) );
0137 }
0138
0139 messenger = new CexmcHistoManagerMessenger( this );
0140 }
0141
0142
0143 CexmcHistoManager::~CexmcHistoManager()
0144 {
0145 if ( outFile )
0146 {
0147 outFile->Write();
0148 outFile->Close();
0149 }
0150
0151
0152 delete outFile;
0153 #ifdef CEXMC_USE_ROOTQT
0154 delete rootCanvas;
0155 #endif
0156 delete messenger;
0157 }
0158
0159
0160 void CexmcHistoManager::AddHisto( const CexmcHistoData & data,
0161 const CexmcAngularRange & aRange )
0162 {
0163 G4String fullName( data.name );
0164 G4String fullTitle( data.title );
0165 G4String rangeTypeLabel;
0166 G4String triggerTypeLabel;
0167 G4String decorTriggerTypeLabel;
0168
0169 if ( data.isARHisto )
0170 {
0171 if ( data.isARRec )
0172 {
0173 fullName += "_arrec";
0174 rangeTypeLabel = "rec";
0175 }
0176 else
0177 {
0178 fullName += "_arreal";
0179 rangeTypeLabel = "real";
0180 }
0181 }
0182
0183 switch ( data.triggerType )
0184 {
0185 case CexmcTPT :
0186 fullName += "_tpt";
0187 decorTriggerTypeLabel = " --tpt--";
0188 fullTitle += decorTriggerTypeLabel;
0189 triggerTypeLabel = "tpt";
0190 break;
0191 case CexmcEDT :
0192 fullName += "_edt";
0193 decorTriggerTypeLabel = " --edt--";
0194 fullTitle += decorTriggerTypeLabel;
0195 triggerTypeLabel = "edt";
0196 break;
0197 case CexmcRT :
0198 fullName += "_rt";
0199 decorTriggerTypeLabel = " --rt--";
0200 fullTitle += decorTriggerTypeLabel;
0201 triggerTypeLabel = "rt";
0202 break;
0203 default :
0204 break;
0205 }
0206
0207 CexmcHistosMap::iterator found( histos.find( data.type ) );
0208
0209 if ( found == histos.end() )
0210 throw CexmcException( CexmcWeirdException );
0211
0212 CexmcHistoVector & histoVector( found->second );
0213
0214 if ( data.isARHisto )
0215 {
0216 G4bool dirOk( false );
0217
0218 dirOk = gDirectory->Get( fullName ) != NULL;
0219
0220 if ( ! dirOk )
0221 dirOk = ( gDirectory->mkdir( fullName, fullTitle ) != NULL );
0222
0223 if ( dirOk )
0224 gDirectory->cd( fullName );
0225
0226 std::ostringstream histoName;
0227 std::ostringstream histoTitle;
0228 histoName << data.name << "_r" << aRange.index + 1 << rangeTypeLabel <<
0229 "_" << triggerTypeLabel;
0230 histoTitle << data.title << " {range " << aRange.index + 1 <<
0231 rangeTypeLabel << " [" << std::fixed <<
0232 std::setprecision( 4 ) << aRange.top << ", " <<
0233 aRange.bottom << ")}" << decorTriggerTypeLabel;
0234 CreateHisto( histoVector, data.impl, histoName.str(), histoTitle.str(),
0235 data.axes );
0236
0237 if ( dirOk )
0238 gDirectory->cd( ".." );
0239 }
0240 else
0241 {
0242 CreateHisto( histoVector, data.impl, fullName, fullTitle, data.axes );
0243 }
0244 }
0245
0246
0247 void CexmcHistoManager::CreateHisto( CexmcHistoVector & histoVector,
0248 CexmcHistoImpl histoImpl, const G4String & name,
0249 const G4String & title, const CexmcHistoAxes & axes )
0250 {
0251 TH1 * histo( NULL );
0252
0253 switch ( histoImpl )
0254 {
0255 case Cexmc_TH1F :
0256 histo = new TH1F( name, title, axes.at( 0 ).nBins,
0257 axes.at( 0 ).nBinsMin, axes.at( 0 ).nBinsMax );
0258 break;
0259 case Cexmc_TH2F :
0260 histo = new TH2F( name, title, axes.at( 0 ).nBins,
0261 axes.at( 0 ).nBinsMin, axes.at( 0 ).nBinsMax,
0262 axes.at( 1 ).nBins, axes.at( 1 ).nBinsMin,
0263 axes.at( 1 ).nBinsMax );
0264 break;
0265 case Cexmc_TH3F :
0266 histo = new TH3F( name, title, axes.at( 0 ).nBins,
0267 axes.at( 0 ).nBinsMin, axes.at( 0 ).nBinsMax,
0268 axes.at( 1 ).nBins, axes.at( 1 ).nBinsMin,
0269 axes.at( 1 ).nBinsMax, axes.at( 2 ).nBins,
0270 axes.at( 2 ).nBinsMin, axes.at( 2 ).nBinsMax );
0271 break;
0272 default :
0273 break;
0274 }
0275
0276 if ( histo )
0277 histoVector.push_back( histo );
0278 }
0279
0280
0281 void CexmcHistoManager::Initialize( void )
0282 {
0283 if ( isInitialized )
0284 return;
0285
0286 CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
0287 G4RunManager::GetRunManager() ) );
0288 CexmcPhysicsManager * physicsManager( runManager->GetPhysicsManager() );
0289
0290 if ( ! physicsManager )
0291 throw CexmcException( CexmcWeirdException );
0292
0293 CexmcProductionModel * productionModel( physicsManager->
0294 GetProductionModel() );
0295
0296 if ( ! productionModel )
0297 throw CexmcException( CexmcWeirdException );
0298
0299 G4ParticleDefinition * outputParticle(
0300 productionModel->GetOutputParticle() );
0301 G4ParticleDefinition * nucleusOutputParticle(
0302 productionModel->GetNucleusOutputParticle() );
0303
0304 if ( ! outputParticle || ! nucleusOutputParticle )
0305 throw CexmcException( CexmcIncompleteProductionModel );
0306
0307 opName = outputParticle->GetParticleName();
0308 nopName = nucleusOutputParticle->GetParticleName();
0309 opMass = outputParticle->GetPDGMass();
0310 nopMass = nucleusOutputParticle->GetPDGMass();
0311
0312 G4String title;
0313 Int_t nBinsX;
0314 Int_t nBinsY;
0315 Double_t nBinsMinX;
0316 Double_t nBinsMaxX;
0317 Double_t nBinsMinY;
0318 Double_t nBinsMaxY;
0319 CexmcHistoAxes axes;
0320
0321 if ( runManager->ProjectIsSaved() )
0322 {
0323 G4String projectsDir( runManager->GetProjectsDir() );
0324 G4String resultsFile( projectsDir + "/" + runManager->GetProjectId() +
0325 ".root" );
0326 outFile = new TFile( resultsFile, "recreate" );
0327 }
0328 else
0329 {
0330 outFile = new TDirectoryFile( CexmcHistoDirectoryHandle,
0331 CexmcHistoDirectoryTitle );
0332 gDirectory->cd( CexmcHistoDirectoryHandle );
0333 }
0334
0335 if ( ! outFile )
0336 throw CexmcException( CexmcWeirdException );
0337
0338 const CexmcSetup * setup( static_cast< const CexmcSetup * >(
0339 runManager->GetUserDetectorConstruction() ) );
0340 if ( ! setup )
0341 throw CexmcException( CexmcWeirdException );
0342
0343 const G4LogicalVolume * lVolume( setup->GetVolume( CexmcSetup::Monitor ) );
0344
0345 if ( ! lVolume )
0346 throw CexmcException( CexmcIncompatibleGeometry );
0347
0348 nBinsMinX = CexmcHistoBeamMomentumMin;
0349 nBinsMaxX = CexmcHistoBeamMomentumMax;
0350 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) /
0351 CexmcHistoBeamMomentumResolution );
0352 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0353 AddHisto( CexmcHistoData( CexmcMomentumBP_TPT_Histo, Cexmc_TH1F, false,
0354 false, CexmcTPT, "mombp", "Beam momentum at the monitor", axes ) );
0355 AddHisto( CexmcHistoData( CexmcMomentumBP_RT_Histo, Cexmc_TH1F, false,
0356 false, CexmcRT, "mombp", "Beam momentum at the monitor", axes ) );
0357 if ( verboseLevel > 0 )
0358 {
0359 AddHisto( CexmcHistoData( CexmcMomentumIP_TPT_Histo, Cexmc_TH1F, false,
0360 false, CexmcTPT, "momip", "Momentum of the incident particle",
0361 axes ) );
0362 }
0363
0364 G4Box * box( dynamic_cast< G4Box * >( lVolume->GetSolid() ) );
0365
0366 if ( ! box )
0367 throw CexmcException( CexmcIncompatibleGeometry );
0368
0369 G4double width( box->GetXHalfLength() * 2 );
0370 G4double height( box->GetYHalfLength() * 2 );
0371 G4double halfWidth( width / 2 + CexmcHistoTPSafetyArea );
0372 G4double halfHeight( height / 2 + CexmcHistoTPSafetyArea );
0373
0374 nBinsX = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
0375 nBinsY = Int_t( halfHeight * 2 / CexmcHistoTPResolution );
0376 axes.clear();
0377 axes.push_back( CexmcHistoAxisData( nBinsX, -halfWidth, halfWidth ) );
0378 axes.push_back( CexmcHistoAxisData( nBinsY, -halfHeight, halfHeight ) );
0379 AddHisto( CexmcHistoData( CexmcTPInMonitor_TPT_Histo, Cexmc_TH2F, false,
0380 false, CexmcTPT, "tpmon", "Track points (mon)", axes ) );
0381
0382 lVolume = setup->GetVolume( CexmcSetup::Target );
0383 G4Tubs * tube( dynamic_cast< G4Tubs * >( lVolume->GetSolid() ) );
0384
0385 if ( ! tube )
0386 throw CexmcException( CexmcIncompatibleGeometry );
0387
0388 G4double radius( tube->GetOuterRadius() );
0389 height = tube->GetZHalfLength() * 2;
0390 halfWidth = radius + CexmcHistoTPSafetyArea;
0391 halfHeight = height / 2 + CexmcHistoTPSafetyArea;
0392
0393 nBinsX = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
0394 nBinsY = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
0395 Int_t nBinsZ( Int_t( halfHeight * 2 / CexmcHistoTPResolution ) );
0396 axes.clear();
0397 axes.push_back( CexmcHistoAxisData( nBinsX, -halfWidth, halfWidth ) );
0398 axes.push_back( CexmcHistoAxisData( nBinsY, -halfWidth, halfWidth ) );
0399 axes.push_back( CexmcHistoAxisData( nBinsZ, -halfHeight, halfHeight ) );
0400 AddHisto( CexmcHistoData( CexmcTPInTarget_TPT_Histo, Cexmc_TH3F, false,
0401 false, CexmcTPT, "tptar", "Track points (tar)", axes ) );
0402 AddHisto( CexmcHistoData( CexmcTPInTarget_RT_Histo, Cexmc_TH3F, false,
0403 false, CexmcRT, "tptar", "Track points (tar)", axes ) );
0404
0405 title = "Reconstructed masses (" + nopName + " vs. " + opName + ")";
0406 nBinsMinX = opMass / 2;
0407 nBinsMaxX = opMass + opMass / 2;
0408 nBinsMinY = nopMass / 2;
0409 nBinsMaxY = nopMass + nopMass / 2;
0410 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoMassResolution );
0411 nBinsY = Int_t( ( nBinsMaxY - nBinsMinY ) / CexmcHistoMassResolution );
0412 axes.clear();
0413 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0414 axes.push_back( CexmcHistoAxisData( nBinsY, nBinsMinY, nBinsMaxY ) );
0415 AddHisto( CexmcHistoData( CexmcRecMasses_EDT_Histo, Cexmc_TH2F, false,
0416 false, CexmcEDT, "recmasses", title, axes ) );
0417 AddHisto( CexmcHistoData( CexmcRecMasses_RT_Histo, Cexmc_TH2F, false,
0418 false, CexmcRT, "recmasses", title, axes ) );
0419
0420 nBinsMinX = 0.;
0421 nBinsMaxX = CexmcHistoEnergyMax;
0422 nBinsMinY = 0.;
0423 nBinsMaxY = CexmcHistoEnergyMax;
0424 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoEnergyResolution );
0425 nBinsY = Int_t( ( nBinsMaxY - nBinsMinY ) / CexmcHistoEnergyResolution );
0426 axes.clear();
0427 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0428 axes.push_back( CexmcHistoAxisData( nBinsY, nBinsMinY, nBinsMaxY ) );
0429 AddHisto( CexmcHistoData( CexmcAbsorbedEnergy_EDT_Histo, Cexmc_TH2F, false,
0430 false, CexmcEDT, "ae", "Absorbed energy (rc vs. lc)", axes ) );
0431 AddHisto( CexmcHistoData( CexmcAbsorbedEnergy_RT_Histo, Cexmc_TH2F, false,
0432 false, CexmcRT, "ae", "Absorbed energy (rc vs. lc)", axes ) );
0433
0434 SetupARHistos( runManager->GetPhysicsManager()->GetProductionModel()->
0435 GetAngularRanges() );
0436
0437 isInitialized = true;
0438 }
0439
0440
0441 void CexmcHistoManager::SetupARHistos( const CexmcAngularRangeList & aRanges )
0442 {
0443 TIter objs( gDirectory->GetList() );
0444 TObject * obj( NULL );
0445
0446 while ( ( obj = ( TObject * )objs() ) )
0447 {
0448 TString name( obj->GetName() );
0449
0450 if ( obj->IsFolder() && ( name.Contains( TString( "_arreal_" ) ) ||
0451 name.Contains( TString( "_arrec_" ) ) ) )
0452 {
0453 gDirectory->cd( name );
0454 gDirectory->DeleteAll();
0455 gDirectory->cd( ".." );
0456 }
0457 }
0458
0459 for ( CexmcHistosMap::iterator k( histos.begin() ); k != histos.end();
0460 ++k )
0461 {
0462 if ( k->second.empty() )
0463 continue;
0464
0465 if ( k->first >= CexmcHistoType_ARReal_START &&
0466 k->first <= CexmcHistoType_ARReal_END )
0467 {
0468 k->second.clear();
0469 }
0470 }
0471
0472 for ( CexmcAngularRangeList::const_iterator k( aRanges.begin() );
0473 k != aRanges.end(); ++k )
0474 {
0475 AddARHistos( *k );
0476 }
0477 }
0478
0479
0480 void CexmcHistoManager::AddARHistos( const CexmcAngularRange & aRange )
0481 {
0482 G4String title;
0483 Int_t nBinsX;
0484 Double_t nBinsMinX;
0485 Double_t nBinsMaxX;
0486 CexmcHistoAxes axes;
0487
0488 title = "Reconstructed mass of " + opName;
0489 nBinsMinX = opMass / 2;
0490 nBinsMaxX = opMass + opMass / 2;
0491 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoMassResolution );
0492 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0493 AddHisto( CexmcHistoData( CexmcRecMassOP_ARReal_RT_Histo, Cexmc_TH1F, true,
0494 false, CexmcRT, "recmassop", title, axes ), aRange );
0495
0496 title = "Reconstructed mass of " + nopName;
0497 nBinsMinX = nopMass / 2;
0498 nBinsMaxX = nopMass + nopMass / 2;
0499 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoMassResolution );
0500 axes.clear();
0501 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0502 AddHisto( CexmcHistoData( CexmcRecMassNOP_ARReal_RT_Histo, Cexmc_TH1F, true,
0503 false, CexmcRT, "recmassnop", title, axes ), aRange );
0504
0505 G4RunManager * runManager( G4RunManager::GetRunManager() );
0506 const CexmcSetup * setup( static_cast< const CexmcSetup * >(
0507 runManager->GetUserDetectorConstruction() ) );
0508 if ( ! setup )
0509 throw CexmcException( CexmcWeirdException );
0510
0511 const G4LogicalVolume * lVolume( setup->GetVolume(
0512 CexmcSetup::Calorimeter ) );
0513
0514 G4Box * box( dynamic_cast< G4Box * >( lVolume->GetSolid() ) );
0515
0516 if ( ! box )
0517 throw CexmcException( CexmcIncompatibleGeometry );
0518
0519 G4double width( box->GetXHalfLength() * 2 );
0520 G4double height( box->GetYHalfLength() * 2 );
0521 G4double halfWidth( width / 2 + CexmcHistoTPSafetyArea );
0522 G4double halfHeight( height / 2 + CexmcHistoTPSafetyArea );
0523
0524 nBinsX = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
0525 Int_t nBinsY( Int_t( halfHeight * 2 / CexmcHistoTPResolution ) );
0526 axes.clear();
0527 axes.push_back( CexmcHistoAxisData( nBinsX, -halfWidth, halfWidth ) );
0528 axes.push_back( CexmcHistoAxisData( nBinsY, -halfHeight, halfHeight ) );
0529
0530
0531
0532
0533 AddHisto( CexmcHistoData( CexmcOPDPAtLeftCalorimeter_ARReal_EDT_Histo,
0534 Cexmc_TH2F, true, false, CexmcEDT, "opdpcl",
0535 "Gamma position on the surface (lc)", axes ), aRange );
0536 AddHisto( CexmcHistoData( CexmcOPDPAtRightCalorimeter_ARReal_EDT_Histo,
0537 Cexmc_TH2F, true, false, CexmcEDT, "opdpcr",
0538 "Gamma position on the surface (rc)", axes ), aRange );
0539 AddHisto( CexmcHistoData( CexmcOPDPAtLeftCalorimeter_ARReal_RT_Histo,
0540 Cexmc_TH2F, true, false, CexmcRT, "opdpcl",
0541 "Gamma position on the surface (lc)", axes ), aRange );
0542 AddHisto( CexmcHistoData( CexmcOPDPAtRightCalorimeter_ARReal_RT_Histo,
0543 Cexmc_TH2F, true, false, CexmcRT, "opdpcr",
0544 "Gamma position on the surface (rc)", axes ), aRange );
0545 AddHisto( CexmcHistoData( CexmcRecOPDPAtLeftCalorimeter_ARReal_EDT_Histo,
0546 Cexmc_TH2F, true, false, CexmcEDT, "recopdpcl",
0547 "Reconstructed gamma position on the surface (lc)", axes ), aRange );
0548 AddHisto( CexmcHistoData( CexmcRecOPDPAtRightCalorimeter_ARReal_EDT_Histo,
0549 Cexmc_TH2F, true, false, CexmcEDT, "recopdpcr",
0550 "Reconstructed gamma position on the surface (rc)", axes ), aRange );
0551 AddHisto( CexmcHistoData( CexmcRecOPDPAtLeftCalorimeter_ARReal_RT_Histo,
0552 Cexmc_TH2F, true, false, CexmcRT, "recopdpcl",
0553 "Reconstructed gamma position on the surface (lc)", axes ), aRange );
0554 AddHisto( CexmcHistoData( CexmcRecOPDPAtRightCalorimeter_ARReal_RT_Histo,
0555 Cexmc_TH2F, true, false, CexmcRT, "recopdpcr",
0556 "Reconstructed gamma position on the surface (rc)", axes ), aRange );
0557
0558 nBinsMinX = 0.;
0559 nBinsMaxX = CexmcHistoEnergyMax;
0560 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoEnergyResolution );
0561 axes.clear();
0562 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0563 AddHisto( CexmcHistoData( CexmcKinEnAtLeftCalorimeter_ARReal_TPT_Histo,
0564 Cexmc_TH1F, true, false, CexmcTPT, "kecl",
0565 "Kinetic energy of gamma (lc)", axes ), aRange );
0566 AddHisto( CexmcHistoData( CexmcKinEnAtRightCalorimeter_ARReal_TPT_Histo,
0567 Cexmc_TH1F, true, false, CexmcTPT, "kecr",
0568 "Kinetic energy of gamma (rc)", axes ), aRange );
0569 AddHisto( CexmcHistoData( CexmcKinEnAtLeftCalorimeter_ARReal_RT_Histo,
0570 Cexmc_TH1F, true, false, CexmcRT, "kecl",
0571 "Kinetic energy of gamma (lc)", axes ), aRange );
0572 AddHisto( CexmcHistoData( CexmcKinEnAtRightCalorimeter_ARReal_RT_Histo,
0573 Cexmc_TH1F, true, false, CexmcRT, "kecr",
0574 "Kinetic energy of gamma (rc)", axes ), aRange );
0575 AddHisto( CexmcHistoData( CexmcAbsEnInLeftCalorimeter_ARReal_EDT_Histo,
0576 Cexmc_TH1F, true, false, CexmcEDT, "aecl",
0577 "Absorbed energy (lc)", axes ), aRange );
0578 AddHisto( CexmcHistoData( CexmcAbsEnInRightCalorimeter_ARReal_EDT_Histo,
0579 Cexmc_TH1F, true, false, CexmcEDT, "aecr",
0580 "Absorbed energy (rc)", axes ), aRange );
0581 AddHisto( CexmcHistoData( CexmcAbsEnInLeftCalorimeter_ARReal_RT_Histo,
0582 Cexmc_TH1F, true, false, CexmcRT, "aecl",
0583 "Absorbed energy (lc)", axes ), aRange );
0584 AddHisto( CexmcHistoData( CexmcAbsEnInRightCalorimeter_ARReal_RT_Histo,
0585 Cexmc_TH1F, true, false, CexmcRT, "aecr",
0586 "Absorbed energy (rc)", axes ), aRange );
0587
0588 nBinsMinX = CexmcHistoMissEnergyMin;
0589 nBinsMaxX = CexmcHistoMissEnergyMax;
0590 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) /
0591 CexmcHistoMissEnergyResolution );
0592 axes.clear();
0593 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0594 AddHisto( CexmcHistoData( CexmcMissEnFromLeftCalorimeter_ARReal_RT_Histo,
0595 Cexmc_TH1F, true, false, CexmcRT, "mecl",
0596 "Missing energy (lc)", axes ), aRange );
0597 AddHisto( CexmcHistoData( CexmcMissEnFromRightCalorimeter_ARReal_RT_Histo,
0598 Cexmc_TH1F, true, false, CexmcRT, "mecr",
0599 "Missing energy (rc)", axes ), aRange );
0600
0601 title = "Kinetic energy of newborn " + opName + " (lab)";
0602 nBinsMinX = 0.;
0603 nBinsMaxX = CexmcHistoEnergyMax;
0604 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoEnergyResolution );
0605 axes.clear();
0606 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0607 AddHisto( CexmcHistoData( CexmcKinEnOP_LAB_ARReal_TPT_Histo,
0608 Cexmc_TH1F, true, false, CexmcTPT, "keop_lab", title, axes ), aRange );
0609 AddHisto( CexmcHistoData( CexmcKinEnOP_LAB_ARReal_RT_Histo,
0610 Cexmc_TH1F, true, false, CexmcRT, "keop_lab", title, axes ), aRange );
0611
0612 title = "Angle of newborn " + opName + " (scm)";
0613 nBinsMinX = -1.0;
0614 nBinsMaxX = 1.0;
0615 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoAngularCResolution );
0616 axes.clear();
0617 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0618 AddHisto( CexmcHistoData( CexmcAngleOP_SCM_ARReal_TPT_Histo,
0619 Cexmc_TH1F, true, false, CexmcTPT, "aop_scm", title, axes ), aRange );
0620 AddHisto( CexmcHistoData( CexmcAngleOP_SCM_ARReal_RT_Histo,
0621 Cexmc_TH1F, true, false, CexmcRT, "aop_scm", title, axes ), aRange );
0622
0623 title = "Reconstruced angle of newborn " + opName + " (scm)";
0624 AddHisto( CexmcHistoData( CexmcRecAngleOP_SCM_ARReal_RT_Histo,
0625 Cexmc_TH1F, true, false, CexmcRT, "recaop_scm", title, axes ), aRange );
0626
0627 title = "Real - reconstruced angle of newborn " + opName + " (scm)";
0628 AddHisto( CexmcHistoData( CexmcDiffAngleOP_SCM_ARReal_RT_Histo,
0629 Cexmc_TH1F, true, false, CexmcRT, "diffaop_scm", title, axes ),
0630 aRange );
0631
0632 nBinsMinX = 0.;
0633 nBinsMaxX = 360.;
0634 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoAngularResolution );
0635 axes.clear();
0636 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0637 AddHisto( CexmcHistoData( CexmcOpenAngle_ARReal_TPT_Histo,
0638 Cexmc_TH1F, true, false, CexmcTPT, "oa",
0639 "Open angle between the gammas", axes ), aRange );
0640 AddHisto( CexmcHistoData( CexmcOpenAngle_ARReal_RT_Histo,
0641 Cexmc_TH1F, true, false, CexmcRT, "oa",
0642 "Open angle between the gammas", axes ), aRange );
0643 AddHisto( CexmcHistoData( CexmcRecOpenAngle_ARReal_RT_Histo,
0644 Cexmc_TH1F, true, false, CexmcRT, "recoa",
0645 "Reconstructed open angle between the gammas", axes ), aRange );
0646
0647 nBinsMinX = -180.;
0648 nBinsMaxX = 180.;
0649 nBinsX = Int_t( ( nBinsMaxX - nBinsMinX ) / CexmcHistoAngularResolution );
0650 axes.clear();
0651 axes.push_back( CexmcHistoAxisData( nBinsX, nBinsMinX, nBinsMaxX ) );
0652 AddHisto( CexmcHistoData( CexmcDiffOpenAngle_ARReal_RT_Histo,
0653 Cexmc_TH1F, true, false, CexmcRT, "diffoa",
0654 "Real - reconstructed open angle between the gammas", axes ), aRange );
0655
0656 lVolume = setup->GetVolume( CexmcSetup::Target );
0657 G4Tubs * tube( dynamic_cast< G4Tubs * >( lVolume->GetSolid() ) );
0658
0659 if ( ! tube )
0660 throw CexmcException( CexmcIncompatibleGeometry );
0661
0662 G4double radius( tube->GetOuterRadius() );
0663 height = tube->GetZHalfLength() * 2;
0664 halfWidth = radius + CexmcHistoTPSafetyArea;
0665 halfHeight = height / 2 + CexmcHistoTPSafetyArea;
0666
0667 nBinsX = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
0668 nBinsY = Int_t( halfWidth * 2 / CexmcHistoTPResolution );
0669 Int_t nBinsZ( Int_t( halfHeight * 2 / CexmcHistoTPResolution ) );
0670 axes.clear();
0671 axes.push_back( CexmcHistoAxisData( nBinsX, -halfWidth, halfWidth ) );
0672 axes.push_back( CexmcHistoAxisData( nBinsY, -halfWidth, halfWidth ) );
0673 axes.push_back( CexmcHistoAxisData( nBinsZ, -halfHeight, halfHeight ) );
0674 AddHisto( CexmcHistoData( CexmcTPInTarget_ARReal_TPT_Histo, Cexmc_TH3F,
0675 true, false, CexmcTPT, "tptar", "Track points (tar)", axes ), aRange );
0676 AddHisto( CexmcHistoData( CexmcTPInTarget_ARReal_RT_Histo, Cexmc_TH3F,
0677 true, false, CexmcRT, "tptar", "Track points (tar)", axes ), aRange );
0678 }
0679
0680
0681 void CexmcHistoManager::Add( CexmcHistoType histoType, unsigned int index,
0682 G4double x )
0683 {
0684 CexmcHistosMap::iterator found( histos.find( histoType ) );
0685 if ( found == histos.end() || histos[ histoType ].size() <= index )
0686 throw CexmcException( CexmcWeirdException );
0687
0688 histos[ histoType ][ index ]->Fill( x );
0689 }
0690
0691
0692 void CexmcHistoManager::Add( CexmcHistoType histoType, unsigned int index,
0693 G4double x, G4double y )
0694 {
0695 CexmcHistosMap::iterator found( histos.find( histoType ) );
0696 if ( found == histos.end() || histos[ histoType ].size() <= index )
0697 throw CexmcException( CexmcWeirdException );
0698
0699
0700
0701 histos[ histoType ][ index ]->Fill( x, y );
0702 }
0703
0704
0705 void CexmcHistoManager::Add( CexmcHistoType histoType, unsigned int index,
0706 G4double x, G4double y, G4double z )
0707 {
0708 CexmcHistosMap::iterator found( histos.find( histoType ) );
0709 if ( found == histos.end() || histos[ histoType ].size() <= index )
0710 throw CexmcException( CexmcWeirdException );
0711
0712
0713
0714 TH3 * histo( static_cast< TH3 * >( histos[ histoType ][ index ] ) );
0715
0716 histo->Fill( x, y, z );
0717 }
0718
0719
0720 void CexmcHistoManager::Add( CexmcHistoType histoType, unsigned int index,
0721 G4int binX, G4int binY, G4double value )
0722 {
0723 CexmcHistosMap::iterator found( histos.find( histoType ) );
0724 if ( found == histos.end() || histos[ histoType ].size() <= index )
0725 throw CexmcException( CexmcWeirdException );
0726
0727 ++binX;
0728 ++binY;
0729 Double_t curValue( histos[ histoType ][ index ]->GetBinContent(
0730 binX, binY ) );
0731 histos[ histoType ][ index ]->SetBinContent( binX, binY,
0732 curValue + value / GeV );
0733 }
0734
0735
0736 void CexmcHistoManager::List( void ) const
0737 {
0738
0739 gDirectory->ls();
0740 }
0741
0742
0743 void CexmcHistoManager::Print( const G4String & value )
0744 {
0745 TObject * histo( gDirectory->FindObjectAny( value.c_str() ) );
0746
0747 if ( ! histo )
0748 {
0749 G4cout << "Histogram '" << value << "' was not found" << G4endl;
0750 return;
0751 }
0752
0753
0754 histo->Print( "range" );
0755 }
0756
0757
0758 #ifdef CEXMC_USE_ROOTQT
0759
0760 void CexmcHistoManager::Draw( const G4String & histoName,
0761 const G4String & histoDrawOptions )
0762 {
0763 if ( ! areLiveHistogramsEnabled )
0764 {
0765 G4cout << "Live histograms option is disabled" << G4endl;
0766 return;
0767 }
0768
0769 TObject * histo( gDirectory->FindObjectAny( histoName ) );
0770
0771 if ( ! histo )
0772 {
0773 G4cout << "Histogram '" << histoName << "' was not found" << G4endl;
0774 return;
0775 }
0776
0777 if ( ! rootCanvas )
0778 {
0779
0780 QFont defaultAppFont( QApplication::font() );
0781 rootCanvas = new CexmcHistoWidget;
0782 QApplication::setFont( defaultAppFont );
0783 rootCanvas->resize( CexmcHistoCanvasWidth, CexmcHistoCanvasHeight );
0784 rootCanvas->GetCanvas()->cd();
0785 }
0786
0787 histo->Draw( histoDrawOptions );
0788 rootCanvas->show();
0789 rootCanvas->GetCanvas()->Update();
0790 }
0791
0792
0793 void CexmcHistoManager::EnableLiveHistograms( G4UIsession * session,
0794 G4bool on )
0795 {
0796 areLiveHistogramsEnabled = on;
0797
0798 if ( ! on || ! isInitialized )
0799 return;
0800
0801 G4UIQt * qtSession( dynamic_cast< G4UIQt * >( session ) );
0802
0803 if ( ! qtSession )
0804 return;
0805
0806 if ( ! histoMenuHandle.empty() && ! isHistoMenuInitialized )
0807 {
0808 qtSession->AddMenu( histoMenuHandle, histoMenuLabel );
0809 BuildMenuTree( qtSession, histoMenuHandle, gDirectory->GetList() );
0810 isHistoMenuInitialized = true;
0811 }
0812 }
0813
0814
0815 void CexmcHistoManager::BuildMenuTree( G4UIQt * session,
0816 const G4String & menu, TList * ls )
0817 {
0818 TIter objs( ls );
0819 TObject * obj( NULL );
0820
0821 while ( ( obj = ( TObject * )objs() ) )
0822 {
0823 G4String name( obj->GetName() );
0824 G4String title( obj->GetTitle() );
0825
0826 if ( obj->IsFolder() )
0827 {
0828 AddSubmenu( session, menu, name, title );
0829 BuildMenuTree( session, name, ( ( TDirectory * )obj )->GetList() );
0830 }
0831 else
0832 {
0833 G4String options( name );
0834
0835 do
0836 {
0837 if ( obj->InheritsFrom( TH3::Class() ) &&
0838 ! drawOptions3D.empty() )
0839 {
0840 title = G4String( "3: " ) + title;
0841 options += G4String( " " ) + drawOptions3D;
0842 break;
0843 }
0844 if ( obj->InheritsFrom( TH2::Class() ) &&
0845 ! drawOptions2D.empty() )
0846 {
0847 title = G4String( "2: " ) + title;
0848 options += G4String( " " ) + drawOptions2D;
0849 break;
0850 }
0851 if ( obj->InheritsFrom( TH1::Class() ) &&
0852 ! drawOptions1D.empty() )
0853 {
0854 options += G4String( " " ) + drawOptions1D;
0855 break;
0856 }
0857 } while ( false );
0858
0859 G4String cmd( CexmcMessenger::histoDirName + "draw " + options );
0860 session->AddButton( menu, title.c_str(), cmd );
0861 }
0862 }
0863 }
0864
0865
0866 void CexmcHistoManager::AddSubmenu( G4UIQt * session,
0867 const G4String & parent,
0868 const G4String & name,
0869 const G4String & label )
0870 {
0871 QMenu * menu( new QMenu( label.c_str() ) );
0872 QMenu * parentMenu( ( QMenu * )session->GetInteractor( parent ) );
0873
0874 parentMenu->addMenu( menu );
0875 session->AddInteractor( name, ( G4Interactor )menu );
0876 }
0877
0878 #endif
0879
0880 #endif
0881
