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