Application/src/RunManagerMT.cc

0001 #include "SimG4Core/Application/interface/RunManagerMT.h"
0002 #include "SimG4Core/Application/interface/PrimaryTransformer.h"
0003 #include "SimG4Core/Application/interface/SimRunInterface.h"
0004 #include "SimG4Core/Application/interface/RunAction.h"
0005 #include "SimG4Core/Application/interface/ParametrisedEMPhysics.h"
0006 #include "SimG4Core/Application/interface/ExceptionHandler.h"
0007
0008 #include "SimG4Core/Geometry/interface/DDDWorld.h"
0009 #include "SimG4Core/Geometry/interface/CustomUIsession.h"
0010
0011 #include "SimG4Core/Physics/interface/PhysicsListFactory.h"
0012 #include "SimG4Core/PhysicsLists/interface/CMSMonopolePhysics.h"
0013 #include "SimG4Core/CustomPhysics/interface/CMSExoticaPhysics.h"
0014
0015 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0016
0017 #include "SimG4Core/Notification/interface/SimTrackManager.h"
0018 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0019 #include "SimG4Core/Notification/interface/CurrentG4Track.h"
0020 #include "SimG4Core/Geometry/interface/CMSG4CheckOverlap.h"
0021
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/Utilities/interface/Exception.h"
0024
0025 #include "DetectorDescription/Core/interface/DDCompactView.h"
0026 #include "DetectorDescription/DDCMS/interface/DDCompactView.h"
0027
0028 #include "SimDataFormats/Forward/interface/LHCTransportLinkContainer.h"
0029
0030 #include "HepPDT/ParticleDataTable.hh"
0031
0032 #include "G4Timer.hh"
0033 #include "G4GeometryManager.hh"
0034 #include "G4ScoringManager.hh"
0035 #include "G4StateManager.hh"
0036 #include "G4ApplicationState.hh"
0037 #include "G4MTRunManagerKernel.hh"
0038 #include "G4UImanager.hh"
0039
0040 #include "G4Run.hh"
0041 #include "G4Event.hh"
0042 #include "G4TransportationManager.hh"
0043 #include "G4ParticleTable.hh"
0044 #include "G4CascadeInterface.hh"
0045 #include "G4EmParameters.hh"
0046 #include "G4HadronicParameters.hh"
0047 #include "G4NuclearLevelData.hh"
0048 #include <CLHEP/Units/SystemOfUnits.h>
0049 #include "G4LogicalVolume.hh"
0050 #include "G4LogicalVolumeStore.hh"
0051 #include "G4PhysicalVolumeStore.hh"
0052 #include "G4Region.hh"
0053 #include "G4RegionStore.hh"
0054
0055 #include <iostream>
0056 #include <sstream>
0057 #include <fstream>
0058 #include <memory>
0059
0060 RunManagerMT::RunManagerMT(edm::ParameterSet const& p)
0061     : m_PhysicsTablesDir(p.getUntrackedParameter<std::string>("PhysicsTablesDirectory", "")),
0062       m_StorePhysicsTables(p.getUntrackedParameter<bool>("StorePhysicsTables", false)),
0063       m_RestorePhysicsTables(p.getUntrackedParameter<bool>("RestorePhysicsTables", false)),
0064       m_check(p.getUntrackedParameter<bool>("CheckGeometry")),
0065       m_geoFromDD4hep(p.getParameter<bool>("g4GeometryDD4hepSource")),
0066       m_score(p.getParameter<bool>("UseCommandBaseScorer")),
0067       m_stepverb(p.getUntrackedParameter<int>("SteppingVerbosity", 0)),
0068       m_regionFile(p.getUntrackedParameter<std::string>("FileNameRegions", "")),
0069       m_pPhysics(p.getParameter<edm::ParameterSet>("Physics")),
0070       m_pRunAction(p.getParameter<edm::ParameterSet>("RunAction")),
0071       m_Init(p.getParameter<edm::ParameterSet>("Init")),
0072       m_G4Commands(p.getParameter<std::vector<std::string> >("G4Commands")) {
0073   m_physicsList.reset(nullptr);
0074   m_world.reset(nullptr);
0075   m_runInterface.reset(nullptr);
0076
0077   m_kernel = new G4MTRunManagerKernel();
0078   m_stateManager = G4StateManager::GetStateManager();
0079   double th = p.getParameter<double>("ThresholdForGeometryExceptions") * CLHEP::GeV;
0080   bool tr = p.getParameter<bool>("TraceExceptions");
0081   m_stateManager->SetExceptionHandler(new ExceptionHandler(th, tr));
0082   if (m_check) {
0083     m_CheckOverlap = p.getUntrackedParameter<edm::ParameterSet>("G4CheckOverlap");
0084   }
0085   m_UIsession = new CustomUIsession();
0086   G4UImanager::GetUIpointer()->SetCoutDestination(m_UIsession);
0087   G4UImanager::GetUIpointer()->SetMasterUIManager(true);
0088 }
0089
0090 RunManagerMT::~RunManagerMT() { delete m_UIsession; }
0091
0092 void RunManagerMT::initG4(const DDCompactView* pDD,
0093                           const cms::DDCompactView* pDD4hep,
0094                           const HepPDT::ParticleDataTable* fPDGTable) {
0095   if (m_managerInitialized) {
0096     edm::LogWarning("SimG4CoreApplication") << "RunManagerMT::initG4 was already done - exit";
0097     return;
0098   }
0099   bool cuts = m_pPhysics.getParameter<bool>("CutsPerRegion");
0100   bool protonCut = m_pPhysics.getParameter<bool>("CutsOnProton");
0101   int verb = m_pPhysics.getUntrackedParameter<int>("Verbosity", 0);
0102   edm::LogVerbatim("SimG4CoreApplication")
0103       << "RunManagerMT: start initialising of geometry DD4hep: " << m_geoFromDD4hep << "\n"
0104       << "              cutsPerRegion: " << cuts << " cutForProton: " << protonCut << "\n"
0105       << "              G4 verbosity: " << verb;
0106
0107   G4Timer timer;
0108   timer.Start();
0109
0110   m_world = std::make_unique<DDDWorld>(pDD, pDD4hep, m_catalog, verb, cuts, protonCut);
0111   G4VPhysicalVolume* world = m_world.get()->GetWorldVolume();
0112
0113   m_kernel->SetVerboseLevel(verb);
0114   edm::LogVerbatim("SimG4CoreApplication")
0115       << "RunManagerMT: Define cuts: " << cuts << " Geant4 run manager verbosity: " << verb;
0116
0117   const G4RegionStore* regStore = G4RegionStore::GetInstance();
0118   const G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
0119   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0120   unsigned int numPV = pvs->size();
0121   unsigned int numLV = lvs->size();
0122   unsigned int nn = regStore->size();
0123   edm::LogVerbatim("SimG4CoreApplication")
0124       << "RunManagerMT: " << numPV << " physical volumes; " << numLV << " logical volumes; " << nn << " regions.";
0125
0126   if (m_check) {
0127     m_kernel->SetVerboseLevel(2);
0128   }
0129   m_kernel->DefineWorldVolume(world, true);
0130   m_registry.dddWorldSignal_(m_world.get());
0131   m_stateManager->SetNewState(G4State_PreInit);
0132
0133   // Create physics list
0134   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT: create PhysicsList";
0135
0136   std::unique_ptr<PhysicsListMakerBase> physicsMaker(
0137       PhysicsListFactory::get()->create(m_pPhysics.getParameter<std::string>("type")));
0138   if (physicsMaker.get() == nullptr) {
0139     throw cms::Exception("Configuration") << "Unable to find the Physics list requested";
0140   }
0141   m_physicsList = physicsMaker->make(m_pPhysics, m_registry);
0142
0143   PhysicsList* phys = m_physicsList.get();
0144   if (phys == nullptr) {
0145     throw cms::Exception("Configuration") << "Physics list construction failed!";
0146   }
0147   if (m_stepverb > 0) {
0148     verb = std::max(verb, 1);
0149   }
0150   G4HadronicParameters::Instance()->SetVerboseLevel(verb);
0151   G4EmParameters::Instance()->SetVerbose(verb);
0152   G4EmParameters::Instance()->SetWorkerVerbose(std::max(verb - 1, 0));
0153
0154   // exotic particle physics
0155   double monopoleMass = m_pPhysics.getUntrackedParameter<double>("MonopoleMass", 0);
0156   if (monopoleMass > 0.0) {
0157     phys->RegisterPhysics(new CMSMonopolePhysics(fPDGTable, m_pPhysics));
0158   }
0159   bool exotica = m_pPhysics.getUntrackedParameter<bool>("ExoticaTransport", false);
0160   if (exotica) {
0161     CMSExoticaPhysics exo(phys, m_pPhysics);
0162   }
0163
0164   // adding GFlash, Russian Roulette for eletrons and gamma,
0165   // step limiters on top of any Physics Lists
0166   phys->RegisterPhysics(new ParametrisedEMPhysics("EMoptions", m_pPhysics));
0167
0168   if (m_RestorePhysicsTables) {
0169     m_physicsList->SetPhysicsTableRetrieved(m_PhysicsTablesDir);
0170   }
0171   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT: start initialisation of PhysicsList for master";
0172
0173   m_physicsList->SetDefaultCutValue(m_pPhysics.getParameter<double>("DefaultCutValue") * CLHEP::cm);
0174   m_physicsList->SetCutsWithDefault();
0175   m_kernel->SetPhysics(phys);
0176
0177   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT: PhysicsList and cuts are defined";
0178
0179   // Enable couple transportation
0180   if (m_score) {
0181     G4ScoringManager* scManager = G4ScoringManager::GetScoringManager();
0182     scManager->SetVerboseLevel(1);
0183   }
0184   // Geant4 UI commands before initialisation of physics
0185   if (!m_G4Commands.empty()) {
0186     edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT: Requested UI commands: ";
0187     for (const std::string& command : m_G4Commands) {
0188       edm::LogVerbatim("SimG4CoreApplication") << "    " << command;
0189       G4UImanager::GetUIpointer()->ApplyCommand(command);
0190     }
0191   }
0192
0193   setupVoxels();
0194
0195   m_stateManager->SetNewState(G4State_Init);
0196   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT: G4State is Init";
0197   m_kernel->InitializePhysics();
0198   if (verb > 0) {
0199     G4EmParameters::Instance()->Dump();
0200   }
0201   m_kernel->SetUpDecayChannels();
0202
0203   if (m_kernel->RunInitialization()) {
0204     m_managerInitialized = true;
0205   } else {
0206     throw cms::Exception("LogicError") << "G4RunManagerKernel initialization failed!";
0207   }
0208
0209   if (m_check) {
0210     checkVoxels();
0211   }
0212
0213   if (m_StorePhysicsTables) {
0214     std::ostringstream dir;
0215     dir << m_PhysicsTablesDir << '\0';
0216     std::string cmd = std::string("/control/shell mkdir -p ") + m_PhysicsTablesDir;
0217     if (!std::ifstream(dir.str().c_str(), std::ios::in))
0218       G4UImanager::GetUIpointer()->ApplyCommand(cmd);
0219     m_physicsList->StorePhysicsTable(m_PhysicsTablesDir);
0220   }
0221   // Appload nuclear level data up to Z=84
0222   G4NuclearLevelData::GetInstance()->UploadNuclearLevelData(84);
0223
0224   if (verb > 1) {
0225     m_physicsList->DumpCutValuesTable();
0226   }
0227   edm::LogVerbatim("SimG4CoreApplication")
0228       << "RunManagerMT: Physics is initilized, now initialise user actions, verb=" << verb;
0229
0230   initializeUserActions();
0231
0232   // G4Region dump file name
0233   runForPhase2();
0234
0235   // Geometry checks
0236   if (m_check || !m_regionFile.empty()) {
0237     CMSG4CheckOverlap check(m_CheckOverlap, m_regionFile, m_UIsession, world);
0238   }
0239
0240   m_stateManager->SetNewState(G4State_PreInit);
0241   G4HadronicParameters::Instance()->SetVerboseLevel(std::max(verb - 1, 0));
0242
0243   // If the Geant4 particle table is needed, decomment the lines below
0244   //
0245   //G4ParticleTable::GetParticleTable()->DumpTable("ALL");
0246   //
0247   m_stateManager->SetNewState(G4State_GeomClosed);
0248   m_currentRun = new G4Run();
0249   m_userRunAction->BeginOfRunAction(m_currentRun);
0250   timer.Stop();
0251   G4cout.precision(4);
0252   G4cout << "RunManagerMT: initG4 done " << timer << G4endl;
0253 }
0254
0255 void RunManagerMT::initializeUserActions() {
0256   m_runInterface = std::make_unique<SimRunInterface>(this, true);
0257   m_userRunAction = new RunAction(m_pRunAction, m_runInterface.get(), true);
0258   Connect(m_userRunAction);
0259 }
0260
0261 void RunManagerMT::Connect(RunAction* runAction) {
0262   runAction->m_beginOfRunSignal.connect(m_registry.beginOfRunSignal_);
0263   runAction->m_endOfRunSignal.connect(m_registry.endOfRunSignal_);
0264 }
0265
0266 void RunManagerMT::stopG4() {
0267   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT::stopG4";
0268   G4GeometryManager::GetInstance()->OpenGeometry();
0269   m_stateManager->SetNewState(G4State_Quit);
0270   if (!m_runTerminated) {
0271     terminateRun();
0272   }
0273   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT::stopG4 done";
0274 }
0275
0276 void RunManagerMT::terminateRun() {
0277   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT::terminateRun";
0278   if (nullptr != m_userRunAction) {
0279     m_userRunAction->EndOfRunAction(m_currentRun);
0280     delete m_userRunAction;
0281     m_userRunAction = nullptr;
0282   }
0283   if (!m_runTerminated) {
0284     m_kernel->RunTermination();
0285   }
0286   m_runTerminated = true;
0287   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT::terminateRun done";
0288 }
0289
0290 void RunManagerMT::checkVoxels() {
0291   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0292   int numLV = lvs->size();
0293   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMT: nLV=" << numLV;
0294   int nvox = 0;
0295   int nslice = 0;
0296   for (int i = 0; i < numLV; ++i) {
0297     auto lv = (*lvs)[i];
0298     auto nd = lv->GetNoDaughters();
0299     auto vox = lv->GetVoxelHeader();
0300     auto sma = lv->GetSmartless();
0301     auto reg = lv->GetRegion();
0302     size_t nsl = (nullptr == vox) ? 0 : vox->GetNoSlices();
0303     if (0 < nsl) {
0304       nslice += nsl;
0305       std::string rname = (nullptr != reg) ? reg->GetName() : "";
0306       edm::LogVerbatim("Voxels") << " " << i << ". Nd=" << nd << " Nsl=" << nsl << " Smartless=" << sma << " "
0307                                  << lv->GetName() << " Region: " << rname;
0308       ++nvox;
0309     }
0310   }
0311   edm::LogVerbatim("SimG4CoreApplication")
0312       << "RunManagerMT: nLV=" << numLV << " NlvVox=" << nvox << " Nslices=" << nslice;
0313 }
0314
0315 void RunManagerMT::setupVoxels() {
0316   double density = m_Init.getParameter<double>("DefaultVoxelDensity");
0317   std::vector<std::string> rnames = m_Init.getParameter<std::vector<std::string> >("VoxelRegions");
0318   std::vector<double> rdensities = m_Init.getParameter<std::vector<double> >("VoxelDensityPerRegion");
0319   int nr = 0;
0320   std::size_t n = rnames.size();
0321   if (n == rdensities.size()) {
0322     nr = (int)n;
0323   }
0324   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0325   for (auto const& lv : *lvs) {
0326     double den = density;
0327     if (0 < nr) {
0328       std::string nam = lv->GetRegion()->GetName();
0329       for (int i = 0; i < nr; ++i) {
0330         if (nam == rnames[i]) {
0331           den = rdensities[i];
0332           break;
0333         }
0334       }
0335     }
0336     lv->SetSmartless(den);
0337   }
0338   edm::LogVerbatim("SimG4CoreApplication")
0339       << "RunManagerMT: default voxel density=" << density << "; number of regions with special density " << nr;
0340 }
0341
0342 void RunManagerMT::runForPhase2() {
0343   const G4RegionStore* regStore = G4RegionStore::GetInstance();
0344   for (auto const& r : *regStore) {
0345     const G4String& name = r->GetName();
0346     if (name == "HGCalRegion" || name == "FastTimerRegionETL" || name == "FastTimerRegionBTL") {
0347       m_isPhase2 = true;
0348       break;
0349     }
0350   }
0351 }