Application/src/RunManagerMTWorker.cc

0001 #include "SimG4Core/Application/interface/RunManagerMTWorker.h"
0002 #include "SimG4Core/Application/interface/RunManagerMT.h"
0003 #include "SimG4Core/Application/interface/SimRunInterface.h"
0004 #include "SimG4Core/Application/interface/RunAction.h"
0005 #include "SimG4Core/Application/interface/EventAction.h"
0006 #include "SimG4Core/Application/interface/Phase2EventAction.h"
0007 #include "SimG4Core/Application/interface/StackingAction.h"
0008 #include "SimG4Core/Application/interface/TrackingAction.h"
0009 #include "SimG4Core/Application/interface/Phase2TrackingAction.h"
0010 #include "SimG4Core/Application/interface/SteppingAction.h"
0011 #include "SimG4Core/Application/interface/Phase2SteppingAction.h"
0012 #include "SimG4Core/Application/interface/CMSSimEventManager.h"
0013 #include "SimG4Core/Application/interface/CustomUIsessionThreadPrefix.h"
0014 #include "SimG4Core/Application/interface/CustomUIsessionToFile.h"
0015 #include "SimG4Core/Application/interface/ExceptionHandler.h"
0016 #include "SimG4Core/Application/interface/CMSGDMLWriteStructure.h"
0017
0018 #include "SimG4Core/Physics/interface/CMSG4TrackInterface.h"
0019 #include "SimG4Core/Geometry/interface/CustomUIsession.h"
0020
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/EventSetup.h"
0025 #include "FWCore/Framework/interface/ConsumesCollector.h"
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027 #include "FWCore/Utilities/interface/Exception.h"
0028 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0029
0030 #include "SimG4Core/Notification/interface/SimActivityRegistry.h"
0031 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0032 #include "SimG4Core/Notification/interface/CMSSteppingVerbose.h"
0033 #include "SimG4Core/Notification/interface/SimTrackManager.h"
0034 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0035
0036 #include "SimG4Core/Geometry/interface/DDDWorld.h"
0037 #include "SimG4Core/MagneticField/interface/FieldBuilder.h"
0038 #include "SimG4Core/MagneticField/interface/CMSFieldManager.h"
0039
0040 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0041 #include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
0042 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0043 #include "HepMC3/Print.h"
0044
0045 #include "SimG4Core/Physics/interface/PhysicsList.h"
0046 #include "SimG4Core/Physics/interface/CMSG4TrackInterface.h"
0047
0048 #include "SimG4Core/SensitiveDetector/interface/AttachSD.h"
0049 #include "SimG4Core/SensitiveDetector/interface/SensitiveTkDetector.h"
0050 #include "SimG4Core/SensitiveDetector/interface/SensitiveCaloDetector.h"
0051 #include "SimG4Core/SensitiveDetector/interface/SensitiveDetectorMakerBase.h"
0052 #include "SimG4Core/SensitiveDetector/interface/sensitiveDetectorMakers.h"
0053
0054 #include "G4Timer.hh"
0055 #include "G4Event.hh"
0056 #include "G4Run.hh"
0057 #include <CLHEP/Units/SystemOfUnits.h>
0058 #include "G4Threading.hh"
0059 #include "G4UImanager.hh"
0060 #include "G4WorkerThread.hh"
0061 #include "G4WorkerRunManagerKernel.hh"
0062 #include "G4StateManager.hh"
0063 #include "G4TransportationManager.hh"
0064 #include "G4LossTableManager.hh"
0065 #include "G4PhysListUtil.hh"
0066 #include "G4Field.hh"
0067 #include "G4FieldManager.hh"
0068 #include "G4ScoringManager.hh"
0069 #include "G4UserSteppingAction.hh"
0070 #include "G4GDMLParser.hh"
0071 #include "G4Threading.hh"
0072
0073 #include <atomic>
0074 #include <memory>
0075 #include <thread>
0076 #include <sstream>
0077 #include <vector>
0078
0079 static std::once_flag applyOnce;
0080 static std::once_flag applyOnceEnd;
0081 static std::once_flag applyOnceGDML;
0082
0083 // from https://hypernews.cern.ch/HyperNews/CMS/get/edmFramework/3302/2.html
0084 namespace {
0085   std::atomic<int> thread_counter{0};
0086
0087   int get_new_thread_index() { return thread_counter++; }
0088
0089   bool createWatchers(const edm::ParameterSet& iP,
0090                       SimActivityRegistry* iReg,
0091                       std::vector<std::shared_ptr<SimWatcher>>& oWatchers,
0092                       std::vector<std::shared_ptr<SimProducer>>& oProds,
0093                       int threadID) {
0094     std::vector<edm::ParameterSet> ws = iP.getParameter<std::vector<edm::ParameterSet>>("Watchers");
0095     if (ws.empty()) {
0096       return false;
0097     }
0098
0099     for (auto& watcher : ws) {
0100       std::unique_ptr<SimWatcherMakerBase> maker(
0101           SimWatcherFactory::get()->create(watcher.getParameter<std::string>("type")));
0102       if (maker == nullptr) {
0103         throw cms::Exception("Configuration")
0104             << "RunManagerMTWorker::createWatchers: "
0105             << "Unable to find the requested Watcher " << watcher.getParameter<std::string>("type");
0106       } else {
0107         std::shared_ptr<SimWatcher> watcherTemp;
0108         std::shared_ptr<SimProducer> producerTemp;
0109         maker->make(watcher, *(iReg), watcherTemp, producerTemp);
0110         if (nullptr != watcherTemp) {
0111           if (!watcherTemp->isMT() && 0 < threadID) {
0112             throw cms::Exception("Configuration")
0113                 << "RunManagerMTWorker::createWatchers: "
0114                 << "Unable to use Watcher " << watcher.getParameter<std::string>("type") << " if number of threads > 1";
0115           } else {
0116             oWatchers.push_back(watcherTemp);
0117             if (nullptr != producerTemp) {
0118               oProds.push_back(producerTemp);
0119             }
0120           }
0121         }
0122       }
0123     }
0124     return (!oWatchers.empty());
0125   }
0126 }  // namespace
0127
0128 struct RunManagerMTWorker::TLSData {
0129   std::unique_ptr<G4RunManagerKernel> kernel;  //must be deleted last
0130   std::unique_ptr<RunAction> userRunAction;
0131   std::unique_ptr<SimRunInterface> runInterface;
0132   std::unique_ptr<SimActivityRegistry> registry;
0133   std::unique_ptr<SimTrackManager> trackManager;
0134   std::vector<SensitiveTkDetector*> sensTkDets;
0135   std::vector<SensitiveCaloDetector*> sensCaloDets;
0136   std::vector<std::shared_ptr<SimWatcher>> watchers;
0137   std::vector<std::shared_ptr<SimProducer>> producers;
0138   // G4Run can only be deleted if there is a G4RunManager
0139   // on the thread where the G4Run is being deleted,
0140   // else it causes a segmentation fault
0141   G4Run* currentRun = nullptr;
0142   std::unique_ptr<G4Event> currentEvent;
0143   edm::RunNumber_t currentRunNumber = 0;
0144   bool threadInitialized = false;
0145   bool runTerminated = false;
0146
0147   TLSData() {}
0148
0149   ~TLSData() = default;
0150 };
0151
0152 // This can not be a smart pointer since we must delete some of the members
0153 // before leaving main() else we get a segmentation fault caused by accessing
0154 // other 'singletons' after those singletons have been deleted. Instead we
0155 // atempt to delete all TLS at RunManagerMTWorker destructor. If that fails for
0156 // some reason, it is better to leak than cause a crash.
0157 // thread_local RunManagerMTWorker::TLSData* RunManagerMTWorker::m_tls{nullptr};
0158
0159 RunManagerMTWorker::RunManagerMTWorker(const edm::ParameterSet& p, edm::ConsumesCollector&& iC)
0160     : m_generator(p.getParameter<edm::ParameterSet>("Generator")),
0161       m_generator3(p.getParameter<edm::ParameterSet>("Generator")),
0162       m_InToken(iC.consumes<edm::HepMCProduct>(
0163           p.getParameter<edm::ParameterSet>("Generator").getParameter<edm::InputTag>("HepMCProductLabel"))),
0164       m_InToken3(iC.consumes<edm::HepMC3Product>(
0165           p.getParameter<edm::ParameterSet>("Generator").getParameter<edm::InputTag>("HepMCProductLabel"))),
0166       m_theLHCTlinkToken(iC.consumes<edm::LHCTransportLinkContainer>(p.getParameter<edm::InputTag>("theLHCTlinkTag"))),
0167       m_nonBeam(p.getParameter<bool>("NonBeamEvent")),
0168       m_UseG4EventManager(p.getParameter<bool>("UseG4EventManager")),
0169       m_pUseMagneticField(p.getParameter<bool>("UseMagneticField")),
0170       m_LHCTransport(p.getParameter<bool>("LHCTransport")),
0171       m_thread_index{(get_new_thread_index())},
0172       m_pField(p.getParameter<edm::ParameterSet>("MagneticField")),
0173       m_pRunAction(p.getParameter<edm::ParameterSet>("RunAction")),
0174       m_pEventAction(p.getParameter<edm::ParameterSet>("EventAction")),
0175       m_pStackingAction(p.getParameter<edm::ParameterSet>("StackingAction")),
0176       m_pTrackingAction(p.getParameter<edm::ParameterSet>("TrackingAction")),
0177       m_pSteppingAction(p.getParameter<edm::ParameterSet>("SteppingAction")),
0178       m_G4Commands(p.getParameter<std::vector<std::string>>("G4Commands")),
0179       m_G4CommandsEndRun(p.getParameter<std::vector<std::string>>("G4CommandsEndRun")),
0180       m_p(p) {
0181   int id = getThreadIndex();
0182   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker for the thread " << id;
0183
0184   // Initialize per-thread output
0185   CMSG4TrackInterface::instance()->setThreadID(id);
0186   G4Threading::G4SetThreadId(id);
0187   G4UImanager::GetUIpointer()->SetUpForAThread(id);
0188   auto iPset = p.getUntrackedParameter<edm::ParameterSet>("CustomUIsession");
0189   const std::string& uitype = iPset.getUntrackedParameter<std::string>("Type");
0190   if (uitype == "MessageLogger") {
0191     m_UIsession = new CustomUIsession();
0192   } else if (uitype == "MessageLoggerThreadPrefix") {
0193     m_UIsession = new CustomUIsessionThreadPrefix(iPset.getUntrackedParameter<std::string>("ThreadPrefix"), id);
0194   } else if (uitype == "FilePerThread") {
0195     m_UIsession = new CustomUIsessionToFile(iPset.getUntrackedParameter<std::string>("ThreadFile"), id);
0196   } else {
0197     throw cms::Exception("Configuration")
0198         << "RunManagerMTWorker::initializeG4: Invalid value of CustomUIsession.Type '" << uitype
0199         << "', valid are MessageLogger, MessageLoggerThreadPrefix, FilePerThread";
0200   }
0201   G4UImanager::GetUIpointer()->SetCoutDestination(m_UIsession);
0202
0203   // sensitive detectors
0204   std::vector<std::string> onlySDs = p.getParameter<std::vector<std::string>>("OnlySDs");
0205   m_sdMakers = sim::sensitiveDetectorMakers(p, iC, onlySDs);
0206
0207   // TLS and watchers
0208   initializeTLS();
0209   if (m_hasWatchers) {
0210     for (auto& watcher : m_tls->watchers) {
0211       watcher->registerConsumes(iC);
0212     }
0213   }
0214   if (m_LHCTransport) {
0215     m_LHCToken = iC.consumes<edm::HepMCProduct>(edm::InputTag("LHCTransport"));
0216   }
0217   if (m_pUseMagneticField) {
0218     m_MagField = iC.esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>();
0219   }
0220   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker is constructed for the thread " << id;
0221   unsigned int k = 0;
0222   for (std::unordered_map<std::string, std::unique_ptr<SensitiveDetectorMakerBase>>::const_iterator itr =
0223            m_sdMakers.begin();
0224        itr != m_sdMakers.end();
0225        ++itr, ++k) {
0226     edm::LogVerbatim("SimG4CoreApplication") << "SD[" << k << "] " << itr->first;
0227   }
0228 }
0229
0230 RunManagerMTWorker::~RunManagerMTWorker() {
0231   m_tls = nullptr;
0232   delete m_UIsession;
0233 }
0234
0235 void RunManagerMTWorker::beginRun(edm::EventSetup const& es) {
0236   int id = getThreadIndex();
0237   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::beginRun for the thread " << id;
0238   for (auto& maker : m_sdMakers) {
0239     maker.second->beginRun(es);
0240   }
0241   if (m_pUseMagneticField) {
0242     m_pMagField = &es.getData(m_MagField);
0243   }
0244   if (m_hasWatchers) {
0245     for (auto& watcher : m_tls->watchers) {
0246       watcher->beginRun(es);
0247     }
0248   }
0249   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::beginRun done for the thread " << id;
0250 }
0251
0252 void RunManagerMTWorker::endRun() {
0253   int id = getThreadIndex();
0254   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::endRun for the thread " << id;
0255   terminateRun();
0256 }
0257
0258 void RunManagerMTWorker::initializeTLS() {
0259   if (nullptr != m_tls) {
0260     return;
0261   }
0262
0263   m_tls = new TLSData();
0264   m_tls->registry = std::make_unique<SimActivityRegistry>();
0265
0266   edm::Service<SimActivityRegistry> otherRegistry;
0267   //Look for an outside SimActivityRegistry
0268   // this is used by the visualization code
0269   int thisID = getThreadIndex();
0270   if (otherRegistry) {
0271     m_tls->registry->connect(*otherRegistry);
0272     if (thisID > 0) {
0273       throw cms::Exception("Configuration")
0274           << "RunManagerMTWorker::initializeTLS : "
0275           << "SimActivityRegistry service (i.e. visualization) is not supported for more than 1 thread. "
0276           << " \n If this use case is needed, RunManagerMTWorker has to be updated.";
0277     }
0278   }
0279   m_hasWatchers = createWatchers(m_p, m_tls->registry.get(), m_tls->watchers, m_tls->producers, thisID);
0280 }
0281
0282 void RunManagerMTWorker::initializeG4(RunManagerMT* runManagerMaster, const edm::EventSetup& es) {
0283   if (m_tls->threadInitialized)
0284     return;
0285
0286   G4Timer timer;
0287   timer.Start();
0288
0289   // I guess everything initialized here should be in thread_local storage
0290   initializeTLS();
0291
0292   int thisID = getThreadIndex();
0293   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::initializeG4 in thread " << thisID << " is started";
0294
0295   // Initialize worker part of shared resources (geometry, physics)
0296   G4WorkerThread::BuildGeometryAndPhysicsVector();
0297
0298   // Create worker run manager
0299   m_tls->kernel.reset(G4WorkerRunManagerKernel::GetRunManagerKernel());
0300   if (nullptr == m_tls->kernel) {
0301     m_tls->kernel = std::make_unique<G4WorkerRunManagerKernel>();
0302   }
0303
0304   // Define G4 exception handler
0305   double th = m_p.getParameter<double>("ThresholdForGeometryExceptions") * CLHEP::GeV;
0306   bool tr = m_p.getParameter<bool>("TraceExceptions");
0307   G4StateManager::GetStateManager()->SetExceptionHandler(new ExceptionHandler(th, tr));
0308
0309   // Set the geometry for the worker, share from master
0310   auto worldPV = runManagerMaster->world().GetWorldVolume();
0311   m_tls->kernel->WorkerDefineWorldVolume(worldPV);
0312   G4TransportationManager* tM = G4TransportationManager::GetTransportationManager();
0313   tM->SetWorldForTracking(worldPV);
0314
0315   // flag of Phase2
0316   m_isPhase2 = runManagerMaster->isPhase2();
0317
0318   // we need the track manager now
0319   int verbose = m_p.getParameter<int>("EventVerbose");
0320   m_tls->trackManager = std::make_unique<SimTrackManager>(&m_simEvent, verbose);
0321
0322   // setup the magnetic field
0323   if (m_pUseMagneticField) {
0324     const GlobalPoint gg(0.f, 0.f, 0.f);
0325
0326     sim::FieldBuilder fieldBuilder(m_pMagField, m_pField);
0327
0328     CMSFieldManager* fieldManager = new CMSFieldManager();
0329     tM->SetFieldManager(fieldManager);
0330     fieldBuilder.build(fieldManager, tM->GetPropagatorInField());
0331
0332     std::string fieldFile = m_p.getUntrackedParameter<std::string>("FileNameField", "");
0333     if (!fieldFile.empty()) {
0334       std::call_once(applyOnce, [this]() { m_dumpMF = true; });
0335       if (m_dumpMF) {
0336         edm::LogVerbatim("SimG4CoreApplication")
0337             << "RunManagerMTWorker::InitializeG4: Dump magnetic field to file " << fieldFile;
0338         DumpMagneticField(tM->GetFieldManager()->GetDetectorField(), fieldFile);
0339       }
0340     }
0341   }
0342
0343   // attach sensitive detector
0344   auto sensDets = sim::attachSD(
0345       m_sdMakers, es, runManagerMaster->catalog(), m_p, m_tls->trackManager.get(), *(m_tls->registry.get()));
0346
0347   m_tls->sensTkDets = sensDets.first;
0348   m_tls->sensCaloDets = sensDets.second;
0349
0350   edm::LogVerbatim("SimG4CoreApplication")
0351       << "RunManagerMTWorker::InitializeG4: Sensitive Detectors are built in thread " << thisID << " found "
0352       << m_tls->sensTkDets.size() << " Tk type SD, and " << m_tls->sensCaloDets.size() << " Calo type SD";
0353
0354   // geometry dump
0355   G4String writeFile = (G4String)m_p.getUntrackedParameter<std::string>("FileNameGDML");
0356   if (!writeFile.empty()) {
0357     std::call_once(applyOnceGDML, [this]() { m_dumpGDML = true; });
0358     edm::LogVerbatim("SimG4CoreApplication") << "DumpGDML:" << m_dumpGDML;
0359     if (m_dumpGDML) {
0360       G4int thID = G4Threading::G4GetThreadId();
0361       edm::LogVerbatim("SimG4CoreApplication") << "ThreadID=" << thID;
0362       G4Threading::G4SetThreadId(-1);
0363       G4GDMLParser gdml;
0364       gdml.SetRegionExport(true);
0365       gdml.SetEnergyCutsExport(true);
0366       gdml.SetSDExport(true);
0367       gdml.Write(writeFile, worldPV, true);
0368       G4Threading::G4SetThreadId(thID);
0369     }
0370   }
0371
0372   // Enable couple transportation
0373   bool scorer = m_p.getParameter<bool>("UseCommandBaseScorer");
0374   if (scorer) {
0375     G4ScoringManager* scManager = G4ScoringManager::GetScoringManager();
0376     scManager->SetVerboseLevel(1);
0377   }
0378
0379   // Set the physics list for the worker, share from master
0380   PhysicsList* physicsList = runManagerMaster->physicsListForWorker();
0381
0382   edm::LogVerbatim("SimG4CoreApplication")
0383       << "RunManagerMTWorker::InitializeG4: start initialisation of PhysicsList for the thread " << thisID;
0384
0385   // Geant4 UI commands in PreInit state
0386   if (!runManagerMaster->G4Commands().empty()) {
0387     G4cout << "RunManagerMTWorker::InitializeG4: Requested UI commands: " << G4endl;
0388     for (const std::string& command : runManagerMaster->G4Commands()) {
0389       G4cout << "          " << command << G4endl;
0390       G4UImanager::GetUIpointer()->ApplyCommand(command);
0391     }
0392   }
0393   G4StateManager::GetStateManager()->SetNewState(G4State_Init);
0394
0395   physicsList->InitializeWorker();
0396   m_tls->kernel->SetPhysics(physicsList);
0397   m_tls->kernel->InitializePhysics();
0398
0399   if (!m_tls->kernel->RunInitialization()) {
0400     throw cms::Exception("Configuration")
0401         << "RunManagerMTWorker::InitializeG4: Geant4 kernel initialization failed in thread " << thisID;
0402   }
0403
0404   //tell all interesting parties that we are beginning the job
0405   BeginOfJob aBeginOfJob(&es);
0406   m_tls->registry->beginOfJobSignal_(&aBeginOfJob);
0407
0408   G4int sv = m_p.getUntrackedParameter<int>("SteppingVerbosity");
0409
0410   if (sv > 0) {
0411     G4double elim = m_p.getUntrackedParameter<double>("StepVerboseThreshold", 0.1) * CLHEP::GeV;
0412     std::vector<int> ve = m_p.getUntrackedParameter<std::vector<int>>("VerboseEvents");
0413     std::vector<int> vn = m_p.getUntrackedParameter<std::vector<int>>("VertexNumber");
0414     std::vector<int> vt = m_p.getUntrackedParameter<std::vector<int>>("VerboseTracks");
0415     m_sVerbose = std::make_unique<CMSSteppingVerbose>(sv, elim, ve, vn, vt);
0416   }
0417   if (!m_UseG4EventManager)
0418     m_evtManager = std::make_unique<CMSSimEventManager>(m_p);
0419   initializeUserActions();
0420
0421   G4StateManager::GetStateManager()->SetNewState(G4State_Idle);
0422
0423   timer.Stop();
0424   edm::LogVerbatim("SimG4CoreApplication")
0425       << "RunManagerMTWorker::initializeG4 done for the thread " << thisID << "  " << timer;
0426   m_tls->threadInitialized = true;
0427 }
0428
0429 void RunManagerMTWorker::initializeUserActions() {
0430   m_tls->runInterface = std::make_unique<SimRunInterface>(this, false);
0431   m_tls->userRunAction = std::make_unique<RunAction>(m_pRunAction, m_tls->runInterface.get(), false);
0432   m_tls->userRunAction->SetMaster(false);
0433   Connect(m_tls->userRunAction.get());
0434
0435   G4int ver = m_p.getParameter<int>("EventVerbose");
0436   G4EventManager* eventManager = m_tls->kernel->GetEventManager();
0437   eventManager->SetVerboseLevel(ver);
0438
0439   // different event actions for Run2,3 and Phase2
0440   G4UserEventAction* userEventAction;
0441   if (m_isPhase2) {
0442     auto ptr =
0443         new Phase2EventAction(m_pEventAction, m_tls->runInterface.get(), m_tls->trackManager.get(), m_sVerbose.get());
0444     Connect(ptr);
0445     userEventAction = (G4UserEventAction*)ptr;
0446   } else {
0447     auto ptr = new EventAction(m_pEventAction, m_tls->runInterface.get(), m_tls->trackManager.get(), m_sVerbose.get());
0448     Connect(ptr);
0449     userEventAction = (G4UserEventAction*)ptr;
0450   }
0451
0452   // different tracking actions for Run2,3 and Phase2
0453   G4UserTrackingAction* userTrackingAction;
0454   if (m_isPhase2) {
0455     auto ptr = new Phase2TrackingAction(m_tls->trackManager.get(), m_sVerbose.get(), m_pTrackingAction);
0456     Connect(ptr);
0457     userTrackingAction = (G4UserTrackingAction*)ptr;
0458   } else {
0459     auto ptr = new TrackingAction(m_tls->trackManager.get(), m_sVerbose.get(), m_pTrackingAction);
0460     Connect(ptr);
0461     userTrackingAction = (G4UserTrackingAction*)ptr;
0462   }
0463
0464   // different stepping actions for Run2,3 and Phase2
0465   G4UserSteppingAction* userSteppingAction;
0466   bool dd4hep = m_p.getParameter<bool>("g4GeometryDD4hepSource");
0467   if (m_isPhase2) {
0468     auto ptr = new Phase2SteppingAction(m_sVerbose.get(), m_pSteppingAction, m_hasWatchers, dd4hep);
0469     Connect(ptr);
0470     userSteppingAction = (G4UserSteppingAction*)ptr;
0471   } else {
0472     auto ptr = new SteppingAction(m_sVerbose.get(), m_pSteppingAction, m_hasWatchers, dd4hep);
0473     Connect(ptr);
0474     userSteppingAction = (G4UserSteppingAction*)ptr;
0475   }
0476
0477   // staking actions and event manager
0478   auto userStackingAction = new StackingAction(m_pStackingAction, m_sVerbose.get());
0479   if (m_UseG4EventManager) {
0480     eventManager->SetUserAction(userEventAction);
0481     eventManager->SetUserAction(userTrackingAction);
0482     eventManager->SetUserAction(userSteppingAction);
0483     eventManager->SetUserAction(userStackingAction);
0484   } else {
0485     m_evtManager->SetUserAction(userEventAction);
0486     m_evtManager->SetUserAction(userTrackingAction);
0487     m_evtManager->SetUserAction(userSteppingAction);
0488     m_evtManager->SetUserAction(userStackingAction);
0489   }
0490 }
0491
0492 void RunManagerMTWorker::Connect(RunAction* runAction) {
0493   runAction->m_beginOfRunSignal.connect(m_tls->registry->beginOfRunSignal_);
0494   runAction->m_endOfRunSignal.connect(m_tls->registry->endOfRunSignal_);
0495 }
0496
0497 void RunManagerMTWorker::Connect(EventAction* eventAction) {
0498   eventAction->m_beginOfEventSignal.connect(m_tls->registry->beginOfEventSignal_);
0499   eventAction->m_endOfEventSignal.connect(m_tls->registry->endOfEventSignal_);
0500 }
0501
0502 void RunManagerMTWorker::Connect(Phase2EventAction* eventAction) {
0503   eventAction->m_beginOfEventSignal.connect(m_tls->registry->beginOfEventSignal_);
0504   eventAction->m_endOfEventSignal.connect(m_tls->registry->endOfEventSignal_);
0505 }
0506
0507 void RunManagerMTWorker::Connect(TrackingAction* trackingAction) {
0508   trackingAction->m_beginOfTrackSignal.connect(m_tls->registry->beginOfTrackSignal_);
0509   trackingAction->m_endOfTrackSignal.connect(m_tls->registry->endOfTrackSignal_);
0510 }
0511
0512 void RunManagerMTWorker::Connect(Phase2TrackingAction* trackingAction) {
0513   trackingAction->m_beginOfTrackSignal.connect(m_tls->registry->beginOfTrackSignal_);
0514   trackingAction->m_endOfTrackSignal.connect(m_tls->registry->endOfTrackSignal_);
0515 }
0516
0517 void RunManagerMTWorker::Connect(SteppingAction* steppingAction) {
0518   steppingAction->m_g4StepSignal.connect(m_tls->registry->g4StepSignal_);
0519 }
0520
0521 void RunManagerMTWorker::Connect(Phase2SteppingAction* steppingAction) {
0522   steppingAction->m_g4StepSignal.connect(m_tls->registry->g4StepSignal_);
0523 }
0524
0525 SimTrackManager* RunManagerMTWorker::getSimTrackManager() {
0526   initializeTLS();
0527   return m_tls->trackManager.get();
0528 }
0529 std::vector<SensitiveTkDetector*>& RunManagerMTWorker::sensTkDetectors() {
0530   initializeTLS();
0531   return m_tls->sensTkDets;
0532 }
0533 std::vector<SensitiveCaloDetector*>& RunManagerMTWorker::sensCaloDetectors() {
0534   initializeTLS();
0535   return m_tls->sensCaloDets;
0536 }
0537 std::vector<std::shared_ptr<SimProducer>>& RunManagerMTWorker::producers() {
0538   initializeTLS();
0539   return m_tls->producers;
0540 }
0541
0542 void RunManagerMTWorker::initializeRun() {
0543   m_tls->currentRun = new G4Run();
0544   G4StateManager::GetStateManager()->SetNewState(G4State_GeomClosed);
0545   if (nullptr != m_tls->userRunAction) {
0546     m_tls->userRunAction->BeginOfRunAction(m_tls->currentRun);
0547   }
0548   int id = getThreadIndex();
0549   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::initializeRun done for thread " << id;
0550 }
0551
0552 void RunManagerMTWorker::terminateRun() {
0553   int id = getThreadIndex();
0554   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::terminateRun for thread " << id;
0555   if (nullptr == m_tls || m_tls->runTerminated) {
0556     return;
0557   }
0558
0559   // Geant4 UI commands after the run
0560   if (!m_G4CommandsEndRun.empty()) {
0561     std::call_once(applyOnceEnd, [this]() { m_endOfRun = true; });
0562     if (m_endOfRun) {
0563       edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker: Requested end of run UI commands: ";
0564       for (const std::string& command : m_G4CommandsEndRun) {
0565         edm::LogVerbatim("SimG4CoreApplication") << "    " << command;
0566         G4UImanager::GetUIpointer()->ApplyCommand(command);
0567       }
0568     }
0569   }
0570   if (m_tls->userRunAction) {
0571     m_tls->userRunAction->EndOfRunAction(m_tls->currentRun);
0572     m_tls->userRunAction.reset();
0573   }
0574   m_tls->currentEvent.reset();
0575   if (m_tls->kernel) {
0576     m_tls->kernel->RunTermination();
0577   }
0578   m_tls->runTerminated = true;
0579   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::terminateRun done for thread " << id;
0580 }
0581
0582 TmpSimEvent* RunManagerMTWorker::produce(const edm::Event& inpevt,
0583                                          const edm::EventSetup& es,
0584                                          RunManagerMT& runManagerMaster) {
0585   // We have to do the per-thread initialization, and per-thread
0586   // per-run initialization here by ourselves.
0587
0588   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::produce: start EventID=" << inpevt.id().event();
0589
0590   assert(m_tls != nullptr and m_tls->threadInitialized);
0591   // Initialize run
0592   if (inpevt.id().run() != m_tls->currentRunNumber) {
0593     edm::LogVerbatim("SimG4CoreApplication")
0594         << "RunManagerMTWorker::produce: RunID= " << inpevt.id().run() << "  TLS RunID= " << m_tls->currentRunNumber;
0595     if (m_tls->currentRunNumber != 0 && !m_tls->runTerminated) {
0596       // If previous run in this thread was not terminated via endRun() call,
0597       // then terminate it now
0598       terminateRun();
0599     }
0600     initializeRun();
0601     m_tls->currentRunNumber = inpevt.id().run();
0602   }
0603
0604   // event and primary
0605
0606   m_tls->currentEvent.reset(generateEvent(inpevt));
0607   m_simEvent.clear();
0608   m_simEvent.setHepEvent(m_generator.genEvent());
0609   m_simEvent.setWeight(m_generator.eventWeight());
0610
0611   if (m_generator.genVertex() != nullptr) {
0612     auto genVertex = m_generator.genVertex();
0613     m_simEvent.collisionPoint(math::XYZTLorentzVectorD(genVertex->x() / CLHEP::cm,
0614                                                        genVertex->y() / CLHEP::cm,
0615                                                        genVertex->z() / CLHEP::cm,
0616                                                        genVertex->t() / CLHEP::second));
0617   }
0618   if (m_tls->currentEvent->GetNumberOfPrimaryVertex() == 0) {
0619     throw cms::Exception("EventCorruption")
0620         << "RunManagerMTWorker::produce: event " << inpevt.id().event() << " with no G4PrimaryVertices"
0621         << " StreamID=" << inpevt.streamID() << " threadIndex=" << getThreadIndex();
0622
0623   } else {
0624     edm::LogVerbatim("SimG4CoreApplication")
0625         << "RunManagerMTWorker::produce: start EventID=" << inpevt.id().event() << " StreamID=" << inpevt.streamID()
0626         << " threadIndex=" << getThreadIndex() << " weight=" << m_simEvent.weight()
0627         << " Nprimary: " << m_tls->currentEvent->GetNumberOfPrimaryVertex() << " vertices and " << m_simEvent.nTracks()
0628         << " particles";
0629     // process event
0630     if (m_UseG4EventManager) {
0631       m_tls->kernel->GetEventManager()->ProcessOneEvent(m_tls->currentEvent.get());
0632     } else {
0633       m_evtManager->ProcessOneEvent(m_tls->currentEvent.get());
0634     }
0635   }
0636
0637   // remove memory only needed during event processing
0638   m_tls->currentEvent.reset();
0639   for (auto& sd : m_tls->sensCaloDets) {
0640     sd->reset();
0641   }
0642
0643   edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::produce: ended Event " << inpevt.id().event();
0644   return &m_simEvent;
0645 }
0646
0647 void RunManagerMTWorker::abortEvent() {
0648   if (m_tls->runTerminated) {
0649     return;
0650   }
0651   G4Track* t = m_tls->kernel->GetEventManager()->GetTrackingManager()->GetTrack();
0652   t->SetTrackStatus(fStopAndKill);
0653
0654   // CMS-specific act
0655   //
0656   TrackingAction* uta = static_cast<TrackingAction*>(m_tls->kernel->GetEventManager()->GetUserTrackingAction());
0657   uta->PostUserTrackingAction(t);
0658
0659   m_tls->currentEvent->SetEventAborted();
0660   m_tls->kernel->GetEventManager()->GetStackManager()->clear();
0661   m_tls->kernel->GetEventManager()->GetTrackingManager()->EventAborted();
0662 }
0663
0664 void RunManagerMTWorker::abortRun(bool softAbort) {
0665   if (!softAbort) {
0666     abortEvent();
0667   }
0668   m_tls->currentRun = nullptr;
0669   terminateRun();
0670 }
0671
0672 G4Event* RunManagerMTWorker::generateEvent(const edm::Event& inpevt) {
0673   m_tls->currentEvent.reset();
0674
0675   // 64 bits event ID in CMSSW converted into Geant4 event ID
0676   G4int evtid = (G4int)inpevt.id().event();
0677   G4Event* evt = new G4Event(evtid);
0678
0679   // required to reset the GenParticle Id for particles transported
0680   // along the beam pipe to their original value for SimTrack creation
0681   resetGenParticleId(inpevt);
0682
0683   edm::Handle<edm::HepMCProduct> HepMCEvt;
0684   bool found = inpevt.getByToken(m_InToken, HepMCEvt);
0685
0686   if (found) {  // HepMC event exists
0687
0688     m_generator.setGenEvent(HepMCEvt->GetEvent());
0689
0690     if (!m_nonBeam) {
0691       m_generator.HepMC2G4(HepMCEvt->GetEvent(), evt);
0692       if (m_LHCTransport) {
0693         edm::Handle<edm::HepMCProduct> LHCMCEvt;
0694         inpevt.getByToken(m_LHCToken, LHCMCEvt);
0695         m_generator.nonCentralEvent2G4(LHCMCEvt->GetEvent(), evt);
0696       }
0697     } else {
0698       m_generator.nonCentralEvent2G4(HepMCEvt->GetEvent(), evt);
0699     }
0700
0701   } else {  // no HepMC event, try to get HepMC3 event
0702
0703     edm::Handle<edm::HepMC3Product> HepMCEvt3;
0704     inpevt.getByToken(m_InToken3, HepMCEvt3);
0705
0706     HepMC3::GenEvent* genevt3 = new HepMC3::GenEvent();
0707     genevt3->read_data(*HepMCEvt3->GetEvent());
0708     m_generator3.setGenEvent(genevt3);
0709
0710     if (!m_nonBeam) {
0711       m_generator3.HepMC2G4(genevt3, evt);
0712       if (m_LHCTransport) {
0713         edm::Handle<edm::HepMC3Product> LHCMCEvt;
0714         inpevt.getByToken(m_LHCToken, LHCMCEvt);
0715         //m_generator3.nonCentralEvent2G4(LHCMCEvt->GetEvent(), evt);
0716       }
0717     } else {
0718       //m_generator3.nonCentralEvent2G4(HepMCEvt->GetEvent(), evt);
0719     }
0720   }
0721   return evt;
0722 }
0723
0724 void RunManagerMTWorker::resetGenParticleId(const edm::Event& inpevt) {
0725   edm::Handle<edm::LHCTransportLinkContainer> theLHCTlink;
0726   inpevt.getByToken(m_theLHCTlinkToken, theLHCTlink);
0727   if (theLHCTlink.isValid()) {
0728     m_tls->trackManager->setLHCTransportLink(theLHCTlink.product());
0729   }
0730 }
0731
0732 void RunManagerMTWorker::DumpMagneticField(const G4Field* field, const std::string& file) const {
0733   std::ofstream fout(file.c_str(), std::ios::out);
0734   if (fout.fail()) {
0735     edm::LogWarning("SimG4CoreApplication")
0736         << "MTWorker::DumpMagneticField: error opening file <" << file << "> for magnetic field";
0737   } else {
0738     // CMS magnetic field volume
0739     double rmax = 9000 * CLHEP::mm;
0740     double zmax = 24000 * CLHEP::mm;
0741
0742     double dr = 1 * CLHEP::cm;
0743     double dz = 5 * CLHEP::cm;
0744
0745     int nr = (int)(rmax / dr);
0746     int nz = 2 * (int)(zmax / dz);
0747
0748     double r = 0.0;
0749     double z0 = -zmax;
0750     double z;
0751
0752     double phi = 0.0;
0753     double cosf = cos(phi);
0754     double sinf = sin(phi);
0755
0756     double point[4] = {0.0, 0.0, 0.0, 0.0};
0757     double bfield[3] = {0.0, 0.0, 0.0};
0758
0759     fout << std::setprecision(6);
0760     for (int i = 0; i <= nr; ++i) {
0761       z = z0;
0762       for (int j = 0; j <= nz; ++j) {
0763         point[0] = r * cosf;
0764         point[1] = r * sinf;
0765         point[2] = z;
0766         field->GetFieldValue(point, bfield);
0767         fout << "R(mm)= " << r / CLHEP::mm << " phi(deg)= " << phi / CLHEP::degree << " Z(mm)= " << z / CLHEP::mm
0768              << "   Bz(tesla)= " << bfield[2] / CLHEP::tesla
0769              << " Br(tesla)= " << (bfield[0] * cosf + bfield[1] * sinf) / CLHEP::tesla
0770              << " Bphi(tesla)= " << (bfield[0] * sinf - bfield[1] * cosf) / CLHEP::tesla << G4endl;
0771         z += dz;
0772       }
0773       r += dr;
0774     }
0775
0776     fout.close();
0777   }
0778 }