g4simulation/g4main/PHG4Reco.cc

0001 #include "PHG4Reco.h"
0002
0003 #include "Fun4AllMessenger.h"
0004 #include "G4TBMagneticFieldSetup.hh"
0005 #include "PHG4DisplayAction.h"
0006 #include "PHG4InEvent.h"
0007 #include "PHG4PhenixDetector.h"
0008 #include "PHG4PhenixDisplayAction.h"
0009 #include "PHG4PhenixEventAction.h"
0010 #include "PHG4PhenixStackingAction.h"
0011 #include "PHG4PhenixSteppingAction.h"
0012 #include "PHG4PhenixTrackingAction.h"
0013 #include "PHG4PrimaryGeneratorAction.h"
0014 #include "PHG4Subsystem.h"
0015 #include "PHG4TrackingAction.h"
0016 #include "PHG4UIsession.h"
0017 #include "PHG4Utils.h"
0018
0019 #include <g4decayer/EDecayType.hh>
0020 #include <g4decayer/P6DExtDecayerPhysics.hh>
0021
0022 #include <g4decayer/EvtGenExtDecayerPhysics.hh>
0023
0024 #include <phgeom/PHGeomUtility.h>
0025
0026 #include <g4gdml/PHG4GDMLUtility.hh>
0027
0028 #include <phfield/PHFieldConfig.h>  // for PHFieldConfig
0029 #include <phfield/PHFieldConfigv1.h>
0030 #include <phfield/PHFieldConfigv2.h>
0031 #include <phfield/PHFieldUtility.h>
0032
0033 #include <ffamodules/CDBInterface.h>
0034
0035 #include <fun4all/Fun4AllReturnCodes.h>
0036 #include <fun4all/Fun4AllServer.h>
0037 #include <fun4all/SubsysReco.h>  // for SubsysReco
0038
0039 #include <phool/PHCompositeNode.h>
0040 #include <phool/PHDataNode.h>      // for PHDataNode
0041 #include <phool/PHNode.h>          // for PHNode
0042 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0043 #include <phool/PHObject.h>        // for PHObject
0044 #include <phool/PHRandomSeed.h>
0045 #include <phool/getClass.h>
0046 #include <phool/phool.h>  // for PHWHERE
0047 #include <phool/recoConsts.h>
0048
0049 #include <TSystem.h>  // for TSystem, gSystem
0050
0051 #include <CLHEP/Random/Random.h>
0052
0053 #include <G4HadronicParameters.hh>  // for G4HadronicParameters
0054 #include <Geant4/G4Cerenkov.hh>
0055 #include <Geant4/G4Element.hh>       // for G4Element
0056 #include <Geant4/G4EventManager.hh>  // for G4EventManager
0057 #include <Geant4/G4HadronicProcessStore.hh>
0058 #include <Geant4/G4IonisParamMat.hh>  // for G4IonisParamMat
0059 #include <Geant4/G4LossTableManager.hh>
0060 #include <Geant4/G4Material.hh>
0061 #include <Geant4/G4NistManager.hh>
0062 #include <Geant4/G4OpAbsorption.hh>
0063 #include <Geant4/G4OpBoundaryProcess.hh>
0064 #include <Geant4/G4OpMieHG.hh>
0065 #include <Geant4/G4OpRayleigh.hh>
0066 #include <Geant4/G4OpWLS.hh>
0067 #include <Geant4/G4OpticalPhoton.hh>
0068 #include <Geant4/G4ParticleDefinition.hh>
0069 #include <Geant4/G4ParticleTable.hh>
0070 #include <Geant4/G4PhotoElectricEffect.hh>  // for G4PhotoElectricEffect
0071 #include <Geant4/G4ProcessManager.hh>
0072 #include <Geant4/G4RunManager.hh>
0073 #include <Geant4/G4Scintillation.hh>
0074 #include <Geant4/G4StepLimiterPhysics.hh>
0075 #include <Geant4/G4String.hh>  // for G4String
0076 #include <Geant4/G4SystemOfUnits.hh>
0077 #include <Geant4/G4Types.hh>  // for G4double, G4int
0078 #include <Geant4/G4UIExecutive.hh>
0079 #include <Geant4/G4UImanager.hh>
0080 #include <Geant4/G4UImessenger.hh>          // for G4UImessenger
0081 #include <Geant4/G4VModularPhysicsList.hh>  // for G4VModularPhysicsList
0082 #include <Geant4/G4Version.hh>
0083 #include <Geant4/G4VisExecutive.hh>
0084 #include <Geant4/G4VisManager.hh>  // for G4VisManager
0085 #include <Geant4/Randomize.hh>     // for G4Random
0086
0087 // physics lists
0088 #include <Geant4/FTFP_BERT.hh>
0089 #include <Geant4/FTFP_BERT_HP.hh>
0090 #include <Geant4/FTFP_INCLXX.hh>
0091 #include <Geant4/FTFP_INCLXX_HP.hh>
0092 #include <Geant4/QGSP_BERT.hh>
0093 #include <Geant4/QGSP_BERT_HP.hh>
0094 #include <Geant4/QGSP_BIC.hh>
0095 #include <Geant4/QGSP_BIC_HP.hh>
0096 #include <Geant4/QGSP_INCLXX.hh>
0097 #include <Geant4/QGSP_INCLXX_HP.hh>
0098
0099 #include <cassert>
0100 #include <cstdlib>
0101 #include <exception>  // for exception
0102 #include <filesystem>
0103 #include <iostream>  // for operator<<, endl
0104 #include <memory>
0105
0106 class G4EmSaturation;
0107 class G4TrackingManager;
0108 class G4VPhysicalVolume;
0109 class PHField;
0110 class PHG4EventAction;
0111 class PHG4StackingAction;
0112 class PHG4SteppingAction;
0113
0114 //_________________________________________________________________
0115 PHG4Reco::PHG4Reco(const std::string &name)
0116   : SubsysReco(name)
0117   , m_Fun4AllMessenger(new Fun4AllMessenger(Fun4AllServer::instance()))
0118 {
0119   return;
0120 }
0121
0122 //_________________________________________________________________
0123 PHG4Reco::~PHG4Reco()
0124 {
0125   // one can delete null pointer (it results in a nop), so checking if
0126   // they are non zero is not needed
0127   delete m_Field;
0128   delete m_RunManager;
0129   delete m_UISession;
0130   delete m_VisManager;
0131   delete m_Fun4AllMessenger;
0132   while (m_SubsystemList.begin() != m_SubsystemList.end())
0133   {
0134     delete m_SubsystemList.back();
0135     m_SubsystemList.pop_back();
0136   }
0137   delete m_DisplayAction;
0138 }
0139
0140 //_________________________________________________________________
0141 int PHG4Reco::Init(PHCompositeNode *topNode)
0142 {
0143   if (Verbosity() > 0)
0144   {
0145     std::cout << "========================= PHG4Reco::Init() ================================" << std::endl;
0146   }
0147   unsigned int iseed = PHRandomSeed();
0148   G4Seed(iseed);  // fixed seed handled in PHRandomSeed()
0149
0150   // create GEANT run manager
0151   if (Verbosity() > 1)
0152   {
0153     std::cout << "PHG4Reco::Init - create run manager" << std::endl;
0154   }
0155
0156   // redirect GEANT verbosity to nowhere
0157   //  if (Verbosity() < 1)
0158   if (false)
0159   {
0160     G4UImanager *uimanager = G4UImanager::GetUIpointer();
0161     m_UISession = new PHG4UIsession();
0162     uimanager->SetSession(m_UISession);
0163     uimanager->SetCoutDestination(m_UISession);
0164   }
0165
0166   m_RunManager = new G4RunManager();
0167
0168   DefineMaterials();
0169   // create physics processes
0170   G4VModularPhysicsList *myphysicslist = nullptr;
0171   if (m_PhysicsList == "FTFP_BERT")
0172   {
0173     myphysicslist = new FTFP_BERT(Verbosity());
0174   }
0175   else if (m_PhysicsList == "FTFP_BERT_HP")
0176   {
0177     setenv("AllowForHeavyElements", "1", 1);
0178     myphysicslist = new FTFP_BERT_HP(Verbosity());
0179   }
0180   else if (m_PhysicsList == "FTFP_INCLXX")
0181   {
0182     myphysicslist = new FTFP_INCLXX(Verbosity());
0183   }
0184   else if (m_PhysicsList == "FTFP_INCLXX_HP")
0185   {
0186     setenv("AllowForHeavyElements", "1", 1);
0187     myphysicslist = new FTFP_INCLXX_HP(Verbosity());
0188   }
0189   else if (m_PhysicsList == "QGSP_BERT")
0190   {
0191     myphysicslist = new QGSP_BERT(Verbosity());
0192   }
0193   else if (m_PhysicsList == "QGSP_BERT_HP")
0194   {
0195     setenv("AllowForHeavyElements", "1", 1);
0196     myphysicslist = new QGSP_BERT_HP(Verbosity());
0197   }
0198   else if (m_PhysicsList == "QGSP_BIC")
0199   {
0200     myphysicslist = new QGSP_BIC(Verbosity());
0201   }
0202   else if (m_PhysicsList == "QGSP_BIC_HP")
0203   {
0204     setenv("AllowForHeavyElements", "1", 1);
0205     myphysicslist = new QGSP_BIC_HP(Verbosity());
0206   }
0207   else if (m_PhysicsList == "QGSP_INCLXX")
0208   {
0209     myphysicslist = new QGSP_INCLXX(Verbosity());
0210   }
0211   else if (m_PhysicsList == "QGSP_INCLXX_HP")
0212   {
0213     setenv("AllowForHeavyElements", "1", 1);
0214     myphysicslist = new QGSP_INCLXX_HP(Verbosity());
0215   }
0216   else
0217   {
0218     std::cout << "Physics List " << m_PhysicsList << " not implemented" << std::endl;
0219     gSystem->Exit(1);
0220     exit(1);
0221   }
0222
0223   if (m_Decayer == kPYTHIA6Decayer)
0224   {
0225     std::cout << "Use PYTHIA Decayer" << std::endl;
0226     G4HadronicParameters::Instance()->SetEnableBCParticles(false);  // Disable the Geant4 built in HF Decay and use external decayers for them
0227     P6DExtDecayerPhysics *decayer = new P6DExtDecayerPhysics();
0228     if (m_ActiveForceDecayFlag)
0229     {
0230       decayer->SetForceDecay(m_ForceDecayType);
0231     }
0232     myphysicslist->RegisterPhysics(decayer);
0233   }
0234
0235   if (m_Decayer == kEvtGenDecayer)
0236   {
0237     std::cout << "Use EvtGen Decayer" << std::endl;
0238     G4HadronicParameters::Instance()->SetEnableBCParticles(false);  // Disable the Geant4 built in HF Decay and use external decayers for them
0239     EvtGenExtDecayerPhysics *decayer = new EvtGenExtDecayerPhysics();
0240     if (CustomizeDecay)
0241     {
0242       decayer->CustomizeEvtGenDecay(EvtGenDecayFile);
0243     }
0244
0245     myphysicslist->RegisterPhysics(decayer);
0246   }
0247
0248   if (m_Decayer == kGEANTInternalDecayer)
0249   {
0250     std::cout << "Use GEANT Internal Decayer" << std::endl;
0251   }
0252
0253   myphysicslist->RegisterPhysics(new G4StepLimiterPhysics());
0254   // initialize cuts so we can ask the world region for it's default
0255   // cuts to propagate them to other regions in DefineRegions()
0256   myphysicslist->SetCutsWithDefault();
0257   m_RunManager->SetUserInitialization(myphysicslist);
0258
0259   DefineRegions();
0260   // initialize registered subsystems
0261   for (SubsysReco *reco : m_SubsystemList)
0262   {
0263     reco->Init(topNode);
0264   }
0265
0266   if (Verbosity() > 0)
0267   {
0268     std::cout << "===========================================================================" << std::endl;
0269   }
0270   return 0;
0271 }
0272
0273 int PHG4Reco::InitField(PHCompositeNode *topNode)
0274 {
0275   if (Verbosity() > 1)
0276   {
0277     std::cout << "PHG4Reco::InitField - create magnetic field setup" << std::endl;
0278     if (std::isfinite(m_MagneticField))
0279     {
0280       std::cout << "using constant file with " << m_MagneticField << " Tesla" << std::endl;
0281     }
0282     else
0283     {
0284       std::cout << "Using fieldmap: " << m_FieldMapFile << std::endl;
0285     }
0286   }
0287   std::unique_ptr<PHFieldConfig> default_field_cfg(nullptr);
0288   if (std::isfinite(m_MagneticField))
0289   {
0290     default_field_cfg.reset(new PHFieldConfigv2(0, 0, m_MagneticField * m_MagneticFieldRescale));
0291   }
0292   else
0293   {
0294     if (std::filesystem::path(m_FieldMapFile).extension() != ".root")
0295     {
0296       // loading from database
0297       m_FieldMapFile = CDBInterface::instance()->getUrl(m_FieldMapFile);
0298     }
0299     if (std::filesystem::exists(m_FieldMapFile))
0300     {
0301       default_field_cfg.reset(new PHFieldConfigv1(m_FieldConfigType, m_FieldMapFile, m_MagneticFieldRescale));
0302     }
0303     else
0304     {
0305       std::cout << PHWHERE << " Fieldmap " << m_FieldMapFile << " not found " << std::endl;
0306     }
0307   }
0308
0309   if (Verbosity() > 1)
0310   {
0311     std::cout << "PHG4Reco::InitField - create magnetic field setup" << std::endl;
0312   }
0313
0314   PHField *phfield = PHFieldUtility::GetFieldMapNode(default_field_cfg.get(), topNode, Verbosity() + 1);
0315   assert(phfield);
0316
0317   m_Field = new G4TBMagneticFieldSetup(phfield);
0318
0319   return Fun4AllReturnCodes::EVENT_OK;
0320 }
0321
0322 int PHG4Reco::InitRun(PHCompositeNode *topNode)
0323 {
0324   // this is a dumb protection against executing this twice.
0325   // we have cases (currently detector display or material scan) where
0326   // we need the detector but have not run any events (who wants to wait
0327   // for processing an event if you just want a detector display?).
0328   // Then the InitRun is executed from a macro. But if you decide to run events
0329   // afterwards, the InitRun is executed by the framework together with all
0330   // other modules. This prevents the PHG4Reco::InitRun() to be executed
0331   // again in this case
0332   static int InitRunDone = 0;
0333   if (InitRunDone)
0334   {
0335     return Fun4AllReturnCodes::EVENT_OK;
0336   }
0337   InitRunDone = 1;
0338   if (Verbosity() > 0)
0339   {
0340     std::cout << "========================= PHG4Reco::InitRun() ================================" << std::endl;
0341   }
0342
0343   recoConsts *rc = recoConsts::instance();
0344
0345   rc->set_StringFlag("WorldMaterial", m_WorldMaterial);
0346   // build world material - so in subsequent code we can call
0347   //  G4Material::GetMaterial(rc->get_StringFlag("WorldMaterial"))
0348   // if the world material is not in the nist DB, we need to implement it here
0349   G4NistManager::Instance()->FindOrBuildMaterial(m_WorldMaterial);
0350   // G4NistManager::Instance()->FindOrBuildMaterial("G4_Galactic");
0351   // G4NistManager::Instance()->FindOrBuildMaterial("G4_Be");
0352   // G4NistManager::Instance()->FindOrBuildMaterial("G4_Al");
0353   // G4NistManager::Instance()->FindOrBuildMaterial("G4_STAINLESS-STEEL");
0354   rc->set_FloatFlag("WorldSizex", m_WorldSize[0]);
0355   rc->set_FloatFlag("WorldSizey", m_WorldSize[1]);
0356   rc->set_FloatFlag("WorldSizez", m_WorldSize[2]);
0357
0358   // setup the global field
0359   const int field_ret = InitField(topNode);
0360   if (field_ret != Fun4AllReturnCodes::EVENT_OK)
0361   {
0362     std::cout << "PHG4Reco::InitRun- Error - Failed field init with status = " << field_ret << std::endl;
0363     return field_ret;
0364   }
0365
0366   // initialize registered subsystems
0367   for (SubsysReco *reco : m_SubsystemList)
0368   {
0369     if (Verbosity() >= 1)
0370     {
0371       std::cout << "PHG4Reco::InitRun - " << reco->Name() << "->InitRun" << std::endl;
0372     }
0373     reco->InitRun(topNode);
0374   }
0375
0376   // create phenix detector, add subsystems, and register to GEANT
0377   // create display settings before detector
0378   m_DisplayAction = new PHG4PhenixDisplayAction(Name());
0379   if (Verbosity() > 1)
0380   {
0381     std::cout << "PHG4Reco::Init - create detector" << std::endl;
0382   }
0383   m_Detector = new PHG4PhenixDetector(this);
0384   m_Detector->Verbosity(Verbosity());
0385   m_Detector->SetWorldSizeX(m_WorldSize[0] * cm);
0386   m_Detector->SetWorldSizeY(m_WorldSize[1] * cm);
0387   m_Detector->SetWorldSizeZ(m_WorldSize[2] * cm);
0388   m_Detector->SetWorldShape(m_WorldShape);
0389   m_Detector->SetWorldMaterial(m_WorldMaterial);
0390
0391   for (PHG4Subsystem *g4sub : m_SubsystemList)
0392   {
0393     if (g4sub->GetDetector())
0394     {
0395       m_Detector->AddDetector(g4sub->GetDetector());
0396     }
0397   }
0398   m_RunManager->SetUserInitialization(m_Detector);
0399
0400   if (m_disableUserActions)
0401   {
0402     std::cout << "PHG4Reco::InitRun - WARNING - event/track/stepping action disabled! "
0403               << "This is aimed to reduce resource consumption for G4 running only. E.g. dose analysis. "
0404               << "Meanwhile, it will disable all Geant4 based analysis. Toggle this feature on/off with PHG4Reco::setDisableUserActions()" << std::endl;
0405   }
0406
0407   setupInputEventNodeReader(topNode);
0408   // create main event action, add subsystemts and register to GEANT
0409   m_EventAction = new PHG4PhenixEventAction();
0410
0411   for (PHG4Subsystem *g4sub : m_SubsystemList)
0412   {
0413     PHG4EventAction *evtact = g4sub->GetEventAction();
0414     if (evtact)
0415     {
0416       m_EventAction->AddAction(evtact);
0417     }
0418   }
0419
0420   if (not m_disableUserActions)
0421   {
0422     m_RunManager->SetUserAction(m_EventAction);
0423   }
0424
0425   // create main stepping action, add subsystems and register to GEANT
0426   m_StackingAction = new PHG4PhenixStackingAction();
0427   for (PHG4Subsystem *g4sub : m_SubsystemList)
0428   {
0429     PHG4StackingAction *action = g4sub->GetStackingAction();
0430     if (action)
0431     {
0432       if (Verbosity() > 1)
0433       {
0434         std::cout << "Adding steppingaction for " << g4sub->Name() << std::endl;
0435       }
0436       m_StackingAction->AddAction(g4sub->GetStackingAction());
0437     }
0438   }
0439
0440   if (not m_disableUserActions)
0441   {
0442     m_RunManager->SetUserAction(m_StackingAction);
0443   }
0444
0445   // create main stepping action, add subsystems and register to GEANT
0446   m_SteppingAction = new PHG4PhenixSteppingAction();
0447   for (PHG4Subsystem *g4sub : m_SubsystemList)
0448   {
0449     PHG4SteppingAction *action = g4sub->GetSteppingAction();
0450     if (action)
0451     {
0452       if (Verbosity() > 1)
0453       {
0454         std::cout << "Adding steppingaction for " << g4sub->Name() << std::endl;
0455       }
0456
0457       m_SteppingAction->AddAction(g4sub->GetSteppingAction());
0458     }
0459   }
0460
0461   if (not m_disableUserActions)
0462   {
0463     m_RunManager->SetUserAction(m_SteppingAction);
0464   }
0465
0466   // create main tracking action, add subsystems and register to GEANT
0467   m_TrackingAction = new PHG4PhenixTrackingAction();
0468   for (PHG4Subsystem *g4sub : m_SubsystemList)
0469   {
0470     m_TrackingAction->AddAction(g4sub->GetTrackingAction());
0471
0472     // not all subsystems define a user tracking action
0473     if (g4sub->GetTrackingAction())
0474     {
0475       // make tracking manager accessible within user tracking action if defined
0476       if (G4TrackingManager *trackingManager = G4EventManager::GetEventManager()->GetTrackingManager())
0477       {
0478         g4sub->GetTrackingAction()->SetTrackingManagerPointer(trackingManager);
0479       }
0480     }
0481   }
0482
0483   if (not m_disableUserActions)
0484   {
0485     m_RunManager->SetUserAction(m_TrackingAction);
0486   }
0487
0488   // initialize
0489   m_RunManager->Initialize();
0490
0491 #if G4VERSION_NUMBER >= 1033
0492   G4EmSaturation *emSaturation = G4LossTableManager::Instance()->EmSaturation();
0493   if (!emSaturation)
0494   {
0495     std::cout << PHWHERE << "Could not initialize EmSaturation, Birks constants will fail" << std::endl;
0496   }
0497 #endif
0498
0499   // add cerenkov and optical photon processes
0500   // std::cout << std::endl << "Ignore the next message - we implemented this correctly" << std::endl;
0501   G4Cerenkov *theCerenkovProcess = new G4Cerenkov("Cerenkov");
0502   // std::cout << "End of bogus warning message" << std::endl << std::endl;
0503   G4Scintillation *theScintillationProcess = new G4Scintillation("Scintillation");
0504
0505   /*
0506     if (Verbosity() > 0)
0507     {
0508     // This segfaults
0509     theCerenkovProcess->DumpPhysicsTable();
0510     }
0511   */
0512   theCerenkovProcess->SetMaxNumPhotonsPerStep(300);
0513   theCerenkovProcess->SetMaxBetaChangePerStep(10.0);
0514   theCerenkovProcess->SetTrackSecondariesFirst(false);  // current PHG4TruthTrackingAction does not support suspect active track and track secondary first
0515 #if G4VERSION_NUMBER < 1100
0516   theScintillationProcess->SetScintillationYieldFactor(1.0);
0517 #endif
0518   theScintillationProcess->SetTrackSecondariesFirst(false);
0519   // theScintillationProcess->SetScintillationExcitationRatio(1.0);
0520
0521   // Use Birks Correction in the Scintillation process
0522
0523   // G4EmSaturation* emSaturation = G4LossTableManager::Instance()->EmSaturation();
0524   // theScintillationProcess->AddSaturation(emSaturation);
0525
0526   G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
0527   G4ParticleTable::G4PTblDicIterator *_theParticleIterator;
0528   _theParticleIterator = theParticleTable->GetIterator();
0529   _theParticleIterator->reset();
0530   while ((*_theParticleIterator)())
0531   {
0532     G4ParticleDefinition *particle = _theParticleIterator->value();
0533     G4String particleName = particle->GetParticleName();
0534     G4ProcessManager *pmanager = particle->GetProcessManager();
0535     if (theCerenkovProcess->IsApplicable(*particle))
0536     {
0537       pmanager->AddProcess(theCerenkovProcess);
0538       pmanager->SetProcessOrdering(theCerenkovProcess, idxPostStep);
0539     }
0540     if (theScintillationProcess->IsApplicable(*particle))
0541     {
0542       pmanager->AddProcess(theScintillationProcess);
0543       pmanager->SetProcessOrderingToLast(theScintillationProcess, idxAtRest);
0544       pmanager->SetProcessOrderingToLast(theScintillationProcess, idxPostStep);
0545     }
0546     for (PHG4Subsystem *g4sub : m_SubsystemList)
0547     {
0548       g4sub->AddProcesses(particle);
0549     }
0550   }
0551   G4ProcessManager *pmanager = G4OpticalPhoton::OpticalPhoton()->GetProcessManager();
0552   // std::cout << " AddDiscreteProcess to OpticalPhoton " << std::endl;
0553   pmanager->AddDiscreteProcess(new G4OpAbsorption());
0554   pmanager->AddDiscreteProcess(new G4OpRayleigh());
0555   pmanager->AddDiscreteProcess(new G4OpMieHG());
0556   pmanager->AddDiscreteProcess(new G4OpBoundaryProcess());
0557   pmanager->AddDiscreteProcess(new G4OpWLS());
0558   pmanager->AddDiscreteProcess(new G4PhotoElectricEffect());
0559   // pmanager->DumpInfo();
0560
0561   // needs large amount of memory which kills central hijing events
0562   // store generated trajectories
0563   // if( G4TrackingManager* trackingManager = G4EventManager::GetEventManager()->GetTrackingManager() ){
0564   //  trackingManager->SetStoreTrajectory( true );
0565   //}
0566
0567   // quiet some G4 print-outs (EM and Hadronic settings during first event)
0568   G4HadronicProcessStore::Instance()->SetVerbose(0);
0569   G4LossTableManager::Instance()->SetVerbose(1);
0570
0571   if ((Verbosity() < 1) && (m_UISession))
0572   {
0573     m_UISession->Verbosity(1);  // let messages after setup come through
0574   }
0575
0576   // Geometry export to DST
0577   if (m_SaveDstGeometryFlag)
0578   {
0579     const std::string filename = PHGeomUtility::GenerateGeometryFileName("gdml");
0580     std::cout << "PHG4Reco::InitRun - export geometry to DST via tmp file " << filename << std::endl;
0581
0582     Dump_GDML(filename);
0583
0584     PHGeomUtility::ImportGeomFile(topNode, filename);
0585
0586     PHGeomUtility::RemoveGeometryFile(filename);
0587   }
0588
0589   if (Verbosity() > 0)
0590   {
0591     std::cout << "===========================================================================" << std::endl;
0592   }
0593
0594   // dump geometry to root file
0595   if (m_ExportGeometry)
0596   {
0597     std::cout << "PHG4Reco::InitRun - writing geometry to " << m_ExportGeomFilename << std::endl;
0598     PHGeomUtility::ExportGeomtry(topNode, m_ExportGeomFilename);
0599   }
0600
0601   if (PHRandomSeed::Verbosity() >= 2)
0602   {
0603     // at high verbosity, to save the random number to file
0604     G4RunManager::GetRunManager()->SetRandomNumberStore(true);
0605   }
0606   return 0;
0607 }
0608
0609 //________________________________________________________________
0610 // Dump TGeo File
0611 void PHG4Reco::Dump_GDML(const std::string &filename)
0612 {
0613   PHG4GDMLUtility ::Dump_GDML(filename, m_Detector->GetPhysicalVolume());
0614 }
0615
0616 //________________________________________________________________
0617 // Dump TGeo File using native Geant4 tools
0618 void PHG4Reco::Dump_G4_GDML(const std::string &filename)
0619 {
0620   PHG4GDMLUtility::Dump_G4_GDML(filename, m_Detector->GetPhysicalVolume());
0621 }
0622
0623 //_________________________________________________________________
0624 int PHG4Reco::ApplyCommand(const std::string &cmd)
0625 {
0626   InitUImanager();
0627   int iret = m_UImanager->ApplyCommand(cmd.c_str());
0628   return iret;
0629 }
0630 //_________________________________________________________________
0631
0632 int PHG4Reco::StartGui()
0633 {
0634   // kludge, using boost::dll::program_location().string().c_str() for the
0635   // program name and putting it into args lead to invalid reads in G4String
0636   char *args[] = {(char *) ("root.exe")};
0637   G4UIExecutive *ui = new G4UIExecutive(1, args, "qt");
0638   InitUImanager();
0639   m_UImanager->ApplyCommand("/control/execute init_gui_vis.mac");
0640   ui->SessionStart();
0641   delete ui;
0642   return 0;
0643 }
0644
0645 int PHG4Reco::InitUImanager()
0646 {
0647   if (!m_UImanager)
0648   {
0649     // Get the pointer to the User Interface manager
0650     // Initialize visualization
0651     m_VisManager = new G4VisExecutive;
0652     m_VisManager->Initialize();
0653     m_UImanager = G4UImanager::GetUIpointer();
0654   }
0655   return 0;
0656 }
0657
0658 //_________________________________________________________________
0659 int PHG4Reco::process_event(PHCompositeNode *topNode)
0660 {
0661   if (PHRandomSeed::Verbosity() >= 2)
0662   {
0663     G4Random::showEngineStatus();
0664   }
0665
0666   // make sure Actions and subsystems have the relevant pointers set
0667   PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0668   m_GeneratorAction->SetInEvent(ineve);
0669
0670   for (SubsysReco *reco : m_SubsystemList)
0671   {
0672     if (Verbosity() >= 2)
0673     {
0674       std::cout << "PHG4Reco::process_event - " << reco->Name() << "->process_event" << std::endl;
0675     }
0676
0677     try
0678     {
0679       reco->process_event(topNode);
0680     }
0681     catch (const std::exception &e)
0682     {
0683       std::cout << PHWHERE << " caught exception thrown during process_event from "
0684                 << reco->Name() << std::endl;
0685       std::cout << "error: " << e.what() << std::endl;
0686       return Fun4AllReturnCodes::ABORTEVENT;
0687     }
0688   }
0689
0690   // run one event
0691   if (Verbosity() >= 2)
0692   {
0693     std::cout << " PHG4Reco::process_event - "
0694               << "run one event :" << std::endl;
0695     ineve->identify();
0696   }
0697   m_RunManager->BeamOn(1);
0698
0699   for (PHG4Subsystem *g4sub : m_SubsystemList)
0700   {
0701     if (Verbosity() >= 2)
0702     {
0703       std::cout << " PHG4Reco::process_event - " << g4sub->Name() << "->process_after_geant" << std::endl;
0704     }
0705     try
0706     {
0707       g4sub->process_after_geant(topNode);
0708     }
0709     catch (const std::exception &e)
0710     {
0711       std::cout << PHWHERE << " caught exception thrown during process_after_geant from "
0712                 << g4sub->Name() << std::endl;
0713       std::cout << "error: " << e.what() << std::endl;
0714       return Fun4AllReturnCodes::ABORTEVENT;
0715     }
0716   }
0717   return 0;
0718 }
0719
0720 int PHG4Reco::ResetEvent(PHCompositeNode *topNode)
0721 {
0722   for (SubsysReco *reco : m_SubsystemList)
0723   {
0724     reco->ResetEvent(topNode);
0725   }
0726   return 0;
0727 }
0728
0729 void PHG4Reco::Print(const std::string &what) const
0730 {
0731   for (SubsysReco *reco : m_SubsystemList)
0732   {
0733     if (what.empty() || what == "ALL" || (reco->Name()).find(what) != std::string::npos)
0734     {
0735       std::cout << "Printing " << reco->Name() << std::endl;
0736       reco->Print(what);
0737       std::cout << "---------------------------" << std::endl;
0738     }
0739   }
0740   return;
0741 }
0742
0743 int PHG4Reco::setupInputEventNodeReader(PHCompositeNode *topNode)
0744 {
0745   PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0746   if (!ineve)
0747   {
0748     PHNodeIterator iter(topNode);
0749     PHCompositeNode *dstNode;
0750     dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0751
0752     ineve = new PHG4InEvent();
0753     PHDataNode<PHObject> *newNode = new PHDataNode<PHObject>(ineve, "PHG4INEVENT", "PHObject");
0754     dstNode->addNode(newNode);
0755   }
0756   // check if we have already registered a generator before creating the default which uses PHG4InEvent Node
0757   if (!m_GeneratorAction)
0758   {
0759     m_GeneratorAction = new PHG4PrimaryGeneratorAction();
0760   }
0761   m_RunManager->SetUserAction(m_GeneratorAction);
0762   return 0;
0763 }
0764
0765 void PHG4Reco::set_rapidity_coverage(const double eta)
0766 {
0767   m_EtaCoverage = eta;
0768   PHG4Utils::SetPseudoRapidityCoverage(eta);
0769 }
0770
0771 void PHG4Reco::G4Seed(const unsigned int i)
0772 {
0773   CLHEP::HepRandom::setTheSeed(i);
0774
0775   if (PHRandomSeed::Verbosity() >= 2)
0776   {
0777     G4Random::showEngineStatus();
0778   }
0779
0780   return;
0781 }
0782
0783 //____________________________________________________________________________
0784 void PHG4Reco::DefineMaterials()
0785 {
0786   G4String symbol, name;  // a=mass of a mole;
0787   G4double density;       // z=mean number of protons;
0788   G4double fractionmass, a;
0789   G4int ncomponents, natoms, z;
0790   // this is for FTFP_BERT_HP where the neutron code barfs
0791   // if the z difference to the last known element (U) is too large
0792   // home made compounds
0793   // this is a legacy implementation but if they are used in multiple
0794   // subsystems put them here
0795   // otherwise implement locally used ones now in the DefineMaterials()
0796   // method in your subsystem
0797
0798   // making quartz
0799   G4Material *quartz = new G4Material("Quartz", density = 2.200 * g / cm3, ncomponents = 2);
0800   quartz->AddElement(G4NistManager::Instance()->FindOrBuildElement("Si"), 1);
0801   quartz->AddElement(G4NistManager::Instance()->FindOrBuildElement("O"), 2);
0802
0803   // making carbon fiber epoxy
0804   G4Material *cfrp_intt = new G4Material("CFRP_INTT", density = 1.69 * g / cm3, ncomponents = 3);
0805   cfrp_intt->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 10);
0806   cfrp_intt->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), 6);
0807   cfrp_intt->AddElement(G4NistManager::Instance()->FindOrBuildElement("O"), 1);
0808
0809   // water glycol mixture for the INTT endcap rings
0810   G4Material *PropyleneGlycol = new G4Material("Propyleneglycol", 1.036 * g / cm3, 3);
0811   PropyleneGlycol->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 3);
0812   PropyleneGlycol->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), 8);
0813   PropyleneGlycol->AddElement(G4NistManager::Instance()->FindOrBuildElement("O"), 2);
0814
0815   G4Material *WaterGlycol_INTT = new G4Material("WaterGlycol_INTT", density = (0.997 * 0.7 + 1.036 * 0.3) * g / cm3, ncomponents = 2);
0816   WaterGlycol_INTT->AddMaterial(PropyleneGlycol, fractionmass = 0.30811936);
0817   WaterGlycol_INTT->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER"), fractionmass = 0.69188064);
0818
0819   // making Rohacell foam 110
0820   G4Material *rohacell_foam_110 = new G4Material("ROHACELL_FOAM_110", density = 0.110 * g / cm3, ncomponents = 4);
0821   rohacell_foam_110->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 8);
0822   rohacell_foam_110->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), 11);
0823   rohacell_foam_110->AddElement(G4NistManager::Instance()->FindOrBuildElement("O"), 2);
0824   rohacell_foam_110->AddElement(G4NistManager::Instance()->FindOrBuildElement("N"), 1);
0825
0826   // making Rohacell foam ROHACELL 51 WF
0827   // Source of density: https://www.rohacell.com/product/peek-industrial/downloads/rohacell%20wf%20product%20information.pdf
0828   G4Material *rohacell_foam_51 = new G4Material("ROHACELL_FOAM_51", density = 0.052 * g / cm3, ncomponents = 4);
0829   rohacell_foam_51->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 8);
0830   rohacell_foam_51->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), 11);
0831   rohacell_foam_51->AddElement(G4NistManager::Instance()->FindOrBuildElement("O"), 2);
0832   rohacell_foam_51->AddElement(G4NistManager::Instance()->FindOrBuildElement("N"), 1);
0833
0834   // making Carbon PEEK : 30 - 70 Vf.
0835   // https://www.quantum-polymers.com/wp-content/uploads/2017/03/QuantaPEEK-CF30.pdf
0836   G4Material *peek = new G4Material("PEEK", density = 1.32 * g / cm3, ncomponents = 3);
0837   peek->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 19);
0838   peek->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), 12);
0839   peek->AddElement(G4NistManager::Instance()->FindOrBuildElement("O"), 3);
0840
0841   G4Material *cf = new G4Material("CF", density = 1.62 * g / cm3, ncomponents = 1);
0842   cf->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 1.);
0843
0844   G4Material *cf30_peek70 = new G4Material("CF30_PEEK70", density = (1.32 * 0.7 + 1.62 * 0.3) * g / cm3, ncomponents = 2);
0845   cf30_peek70->AddMaterial(cf, fractionmass = 0.34468085);
0846   cf30_peek70->AddMaterial(peek, fractionmass = 0.65531915);
0847
0848   // that seems to be the composition of 304 Stainless steel
0849   G4Material *StainlessSteel =
0850       new G4Material("SS304", density = 7.9 * g / cm3, ncomponents = 8);
0851   StainlessSteel->AddElement(G4NistManager::Instance()->FindOrBuildElement("Fe"), 0.70105);
0852   StainlessSteel->AddElement(G4NistManager::Instance()->FindOrBuildElement("Cr"), 0.18);
0853   StainlessSteel->AddElement(G4NistManager::Instance()->FindOrBuildElement("Ni"), 0.09);
0854   StainlessSteel->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mn"), 0.02);
0855   StainlessSteel->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 0.0007);
0856   StainlessSteel->AddElement(G4NistManager::Instance()->FindOrBuildElement("S"), 0.0003);
0857   StainlessSteel->AddElement(G4NistManager::Instance()->FindOrBuildElement("Si"), 0.0075);
0858   StainlessSteel->AddElement(G4NistManager::Instance()->FindOrBuildElement("P"), 0.00045);
0859
0860   G4Material *SS310 =
0861       new G4Material("SS310", density = 8.0 * g / cm3, ncomponents = 8);
0862   SS310->AddElement(G4NistManager::Instance()->FindOrBuildElement("Fe"), 0.50455);
0863   SS310->AddElement(G4NistManager::Instance()->FindOrBuildElement("Cr"), 0.25);
0864   SS310->AddElement(G4NistManager::Instance()->FindOrBuildElement("Ni"), 0.20);
0865   SS310->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mn"), 0.02);
0866   SS310->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 0.0025);
0867   SS310->AddElement(G4NistManager::Instance()->FindOrBuildElement("S"), 0.015);
0868   SS310->AddElement(G4NistManager::Instance()->FindOrBuildElement("Si"), 0.0075);
0869   SS310->AddElement(G4NistManager::Instance()->FindOrBuildElement("P"), 0.00045);
0870
0871   // SS316 from https://www.azom.com
0872   G4Material *SS316 =
0873       new G4Material("SS316", density = 8.0 * g / cm3, ncomponents = 9);
0874   SS316->AddElement(G4NistManager::Instance()->FindOrBuildElement("Fe"), 0.68095);
0875   SS316->AddElement(G4NistManager::Instance()->FindOrBuildElement("Cr"), 0.16);
0876   SS316->AddElement(G4NistManager::Instance()->FindOrBuildElement("Ni"), 0.11);
0877   SS316->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mn"), 0.02);
0878   SS316->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mo"), 0.02);
0879   SS316->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 0.0008);
0880   SS316->AddElement(G4NistManager::Instance()->FindOrBuildElement("S"), 0.0003);
0881   SS316->AddElement(G4NistManager::Instance()->FindOrBuildElement("Si"), 0.0075);
0882   SS316->AddElement(G4NistManager::Instance()->FindOrBuildElement("P"), 0.00045);
0883
0884   G4Material *Steel =
0885       new G4Material("Steel", density = 7.86 * g / cm3, ncomponents = 5);
0886   Steel->AddElement(G4NistManager::Instance()->FindOrBuildElement("Fe"), 0.9834);
0887   Steel->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mn"), 0.014);
0888   Steel->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 0.0017);
0889   Steel->AddElement(G4NistManager::Instance()->FindOrBuildElement("S"), 0.00045);
0890   Steel->AddElement(G4NistManager::Instance()->FindOrBuildElement("P"), 0.00045);
0891
0892   // a36 steel from http://www.matweb.com
0893   G4Material *a36 = new G4Material("Steel_A36", density = 7.85 * g / cm3, ncomponents = 5);
0894   a36->AddElement(G4NistManager::Instance()->FindOrBuildElement("Fe"), 0.9824);
0895   a36->AddElement(G4NistManager::Instance()->FindOrBuildElement("Cu"), 0.002);
0896   a36->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 0.0025);
0897   a36->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mn"), 0.0103);
0898   a36->AddElement(G4NistManager::Instance()->FindOrBuildElement("Si"), 0.0028);
0899
0900   // 1006 steel from http://www.matweb.com
0901   G4Material *steel_1006 = new G4Material("Steel_1006", density = 7.872 * g / cm3, ncomponents = 2);
0902   steel_1006->AddElement(G4NistManager::Instance()->FindOrBuildElement("Fe"), 0.996);
0903   steel_1006->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mn"), 0.004);
0904
0905   // from www.aalco.co.uk
0906   G4Material *Al5083 = new G4Material("Al5083", density = 2.65 * g / cm3, ncomponents = 3);
0907   Al5083->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mn"), 0.004);
0908   Al5083->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mg"), 0.04);
0909   Al5083->AddElement(G4NistManager::Instance()->FindOrBuildElement("Al"), 0.956);
0910
0911   // Al 4046 from http://www.matweb.com
0912   G4Material *Al4046 = new G4Material("Al4046", density = 2.66 * g / cm3, ncomponents = 3);
0913   Al4046->AddElement(G4NistManager::Instance()->FindOrBuildElement("Al"), 0.897);
0914   Al4046->AddElement(G4NistManager::Instance()->FindOrBuildElement("Si"), 0.1);
0915   Al4046->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mg"), 0.003);
0916
0917   // Al 6061T6 from http://www.matweb.com
0918   G4Material *Al6061T6 = new G4Material("Al6061T6", density = 2.70 * g / cm3, ncomponents = 4);
0919   Al6061T6->AddElement(G4NistManager::Instance()->FindOrBuildElement("Al"), 0.975);
0920   Al6061T6->AddElement(G4NistManager::Instance()->FindOrBuildElement("Si"), 0.008);
0921   Al6061T6->AddElement(G4NistManager::Instance()->FindOrBuildElement("Mg"), 0.01);
0922   Al6061T6->AddElement(G4NistManager::Instance()->FindOrBuildElement("Fe"), 0.007);
0923
0924   // E864 Pb-Scifi calorimeter
0925   // E864 Calorimeter is 99% Pb, 1% Antimony
0926   // Nuclear Instruments and Methods in Physics Research A 406 (1998) 227 258
0927   G4double density_e864 = (0.99 * 11.34 + 0.01 * 6.697) * g / cm3;
0928   G4Material *absorber_e864 = new G4Material("E864_Absorber", density_e864, 2);
0929   absorber_e864->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_Pb"), 0.99);
0930   absorber_e864->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_Sb"), 0.01);
0931
0932   G4Material *FPC = new G4Material("FPC", 1.542 * g / cm3, 2);
0933   FPC->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_Cu"), 0.0162);
0934   FPC->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_KAPTON"), 0.9838);
0935
0936   // This is an approximation for the W saturated epoxy of the EMCal.
0937   G4Material *W_Epoxy = new G4Material("W_Epoxy", density = 10.2 * g / cm3, ncomponents = 2);
0938   W_Epoxy->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_W"), fractionmass = 0.5);
0939   W_Epoxy->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_POLYSTYRENE"), fractionmass = 0.5);
0940
0941   // from http://www.physi.uni-heidelberg.de/~adler/TRD/TRDunterlagen/RadiatonLength/tgc2.htm
0942   // Epoxy (for FR4 )
0943   // density = 1.2*g/cm3;
0944   G4Material *Epoxy = new G4Material("Epoxy", 1.2 * g / cm3, ncomponents = 2);
0945   Epoxy->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), natoms = 2);
0946   Epoxy->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), natoms = 2);
0947
0948   // FR4 (Glass + Epoxy)
0949   density = 1.86 * g / cm3;
0950   G4Material *FR4 = new G4Material("FR4", density, ncomponents = 2);
0951   FR4->AddMaterial(quartz, fractionmass = 0.528);
0952   FR4->AddMaterial(Epoxy, fractionmass = 0.472);
0953   // NOMEX (HoneyComb)
0954   // density from http://www.fibreglast.com/product/Nomex_Honeycomb_1562/Vacuum_Bagging_Sandwich_Core
0955   // 1562: 29 kg/m^3 <-- I guess it is this one
0956   // 2562: 48 kg/m^3
0957   // chemical composition http://ww2.unime.it/cdlchimind/adm/inviofile/uploads/HP_Pols2.b.pdf
0958   density = 29 * kg / m3;
0959   G4Material *NOMEX = new G4Material("NOMEX", density, ncomponents = 4);
0960   NOMEX->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), natoms = 14);
0961   NOMEX->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), natoms = 10);
0962   NOMEX->AddElement(G4NistManager::Instance()->FindOrBuildElement("N"), natoms = 2);
0963   NOMEX->AddElement(G4NistManager::Instance()->FindOrBuildElement("O"), natoms = 2);
0964   // spacal material. Source : EICROOT/A. Kiselev
0965   /*
0966   WEpoxyMix          3  12.011 1.008 183.85  6.  1.  74.  12.18  0.029 0.002 0.969
0967          1  1  30.  .00001
0968                      0
0969                      */
0970   G4Material *Spacal_W_Epoxy =
0971       new G4Material("Spacal_W_Epoxy", density = 12.18 * g / cm3, ncomponents = 3);
0972   Spacal_W_Epoxy->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 0.029);
0973   Spacal_W_Epoxy->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), 0.002);
0974   Spacal_W_Epoxy->AddElement(G4NistManager::Instance()->FindOrBuildElement("W"), 0.969);
0975   /*
0976 PMMA      -3  12.01 1.008 15.99  6.  1.  8.  1.19  3.6  5.7  1.4
0977        1  1  20.  .00001
0978                    0
0979                      */
0980   G4Material *PMMA =
0981       new G4Material("PMMA", density = 1.19 * g / cm3, ncomponents = 3);
0982   PMMA->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), 3.6 / (3.6 + 5.7 + 1.4));
0983   PMMA->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), 5.7 / (3.6 + 5.7 + 1.4));
0984   PMMA->AddElement(G4NistManager::Instance()->FindOrBuildElement("O"), 1.4 / (3.6 + 5.7 + 1.4));
0985
0986   // scintillator for HCal, use a new name in order to change the Birks' constant
0987   G4Material *Uniplast_scintillator = new G4Material("Uniplast_scintillator", 1.06 * g / cm3, ncomponents = 1);
0988   Uniplast_scintillator->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_POLYSTYRENE"), fractionmass = 1.);
0989
0990   G4Material *G10 =
0991       new G4Material("G10", density = 1.700 * g / cm3, ncomponents = 4);
0992   G10->AddElement(G4NistManager::Instance()->FindOrBuildElement("Si"), natoms = 1);
0993   G10->AddElement(G4NistManager::Instance()->FindOrBuildElement("O"), natoms = 2);
0994   G10->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), natoms = 3);
0995   G10->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), natoms = 3);
0996
0997   G4Material *CsI =
0998       new G4Material("CsI", density = 4.534 * g / cm3, ncomponents = 2);
0999   CsI->AddElement(G4NistManager::Instance()->FindOrBuildElement("Cs"), natoms = 1);
1000   CsI->AddElement(G4NistManager::Instance()->FindOrBuildElement("I"), natoms = 1);
1001   CsI->GetIonisation()->SetMeanExcitationEnergy(553.1 * eV);
1002
1003   G4Material *C4F10 =
1004       new G4Material("C4F10", density = 0.00973 * g / cm3, ncomponents = 2);
1005   C4F10->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), natoms = 4);
1006   C4F10->AddElement(G4NistManager::Instance()->FindOrBuildElement("F"), natoms = 10);
1007
1008   G4Element *elLu = new G4Element(name = "Lutetium", symbol = "Lu", z = 71., a = 174.97 * g / mole);
1009   G4Material *LSO = new G4Material("LSO",                    // its name
1010                                    density = 7.4 * g / cm3,  // its density
1011                                    ncomponents = 3);         // number of components
1012
1013   LSO->AddElement(G4NistManager::Instance()->FindOrBuildElement("Si"), natoms = 1);
1014   LSO->AddElement(elLu, natoms = 2);
1015   LSO->AddElement(G4NistManager::Instance()->FindOrBuildElement("O"), natoms = 5);
1016
1017   // Silver epoxy glue LOCTITE ABLESTIK 2902 for the silicon sensors and FPHX chips of INTT
1018   G4Material *SilverEpoxyGlue_INTT = new G4Material("SilverEpoxyGlue_INTT", density = 3.2 * g / cm3, ncomponents = 2);
1019   SilverEpoxyGlue_INTT->AddMaterial(Epoxy, fractionmass = 0.79);
1020   SilverEpoxyGlue_INTT->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_Ag"), fractionmass = 0.21);
1021
1022   // define P10 Gas which will be used for TPC Benchmarking
1023   G4Material *P10 =
1024       new G4Material("P10", density = 1.74 * mg / cm3, ncomponents = 3);  // @ 0K, 1atm
1025   P10->AddElement(G4NistManager::Instance()->FindOrBuildElement("Ar"), fractionmass = 0.9222);
1026   P10->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), fractionmass = 0.0623);
1027   P10->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), fractionmass = 0.0155);
1028 }
1029
1030 void PHG4Reco::DefineRegions()
1031 {
1032   // the PAI model does not work anymore in G4 10.06
1033   //   const G4RegionStore *theRegionStore = G4RegionStore::GetInstance();
1034   //   G4ProductionCuts *gcuts = new G4ProductionCuts(*(theRegionStore->GetRegion("DefaultRegionForTheWorld")->GetProductionCuts()));
1035   //   G4Region *tpcregion = new G4Region("REGION_TPCGAS");
1036   //   tpcregion->SetProductionCuts(gcuts);
1037   // #if G4VERSION_NUMBER >= 1033
1038   //   // Use this from the new G4 version 10.03 on
1039   //   // was commented out, crashes in 10.06 I think
1040   //   // add the PAI model to the TPCGAS region
1041   //   // undocumented, painfully digged out with debugger by tracing what
1042   //   // is done for command "/process/em/AddPAIRegion all TPCGAS PAI"
1043   // //  G4EmParameters *g4emparams = G4EmParameters::Instance();
1044   // //  g4emparams->AddPAIModel("all", "REGION_TPCGAS", "PAI");
1045   // #endif
1046   return;
1047 }
1048
1049 PHG4Subsystem *
1050 PHG4Reco::getSubsystem(const std::string &name)
1051 {
1052   for (PHG4Subsystem *subsys : m_SubsystemList)
1053   {
1054     if (subsys->Name() == name)
1055     {
1056       if (Verbosity() > 0)
1057       {
1058         std::cout << "Found Subsystem " << name << std::endl;
1059       }
1060       return subsys;
1061     }
1062   }
1063   std::cout << "Could not find Subsystem " << name << std::endl;
1064   return nullptr;
1065 }
1066
1067 void PHG4Reco::G4Verbosity(const int i)
1068 {
1069   if (m_RunManager)
1070   {
1071     m_RunManager->SetVerboseLevel(i);
1072   }
1073 }
1074
1075 void PHG4Reco::ApplyDisplayAction()
1076 {
1077   if (!m_Detector)
1078   {
1079     return;
1080   }
1081   G4VPhysicalVolume *physworld = m_Detector->GetPhysicalVolume();
1082   m_DisplayAction->ApplyDisplayAction(physworld);
1083   for (PHG4Subsystem *g4sub : m_SubsystemList)
1084   {
1085     PHG4DisplayAction *action = g4sub->GetDisplayAction();
1086     if (action)
1087     {
1088       if (Verbosity() > 1)
1089       {
1090         std::cout << "Adding steppingaction for " << g4sub->Name() << std::endl;
1091       }
1092       action->ApplyDisplayAction(physworld);
1093     }
1094   }
1095 }