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