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