File indexing completed on 2025-04-09 02:02:20
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
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 }
0127
0128 struct RunManagerMTWorker::TLSData {
0129 std::unique_ptr<G4RunManagerKernel> kernel;
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
0139
0140
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
0153
0154
0155
0156
0157
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
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
0204 std::vector<std::string> onlySDs = p.getParameter<std::vector<std::string>>("OnlySDs");
0205 m_sdMakers = sim::sensitiveDetectorMakers(p, iC, onlySDs);
0206
0207
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
0268
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
0290 initializeTLS();
0291
0292 int thisID = getThreadIndex();
0293 edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::initializeG4 in thread " << thisID << " is started";
0294
0295
0296 G4WorkerThread::BuildGeometryAndPhysicsVector();
0297
0298
0299 m_tls->kernel.reset(G4WorkerRunManagerKernel::GetRunManagerKernel());
0300 if (nullptr == m_tls->kernel) {
0301 m_tls->kernel = std::make_unique<G4WorkerRunManagerKernel>();
0302 }
0303
0304
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
0310 auto worldPV = runManagerMaster->world().GetWorldVolume();
0311 m_tls->kernel->WorkerDefineWorldVolume(worldPV);
0312 G4TransportationManager* tM = G4TransportationManager::GetTransportationManager();
0313 tM->SetWorldForTracking(worldPV);
0314
0315
0316 m_isPhase2 = runManagerMaster->isPhase2();
0317
0318
0319 int verbose = m_p.getParameter<int>("EventVerbose");
0320 m_tls->trackManager = std::make_unique<SimTrackManager>(&m_simEvent, verbose);
0321
0322
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
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
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
0373 bool scorer = m_p.getParameter<bool>("UseCommandBaseScorer");
0374 if (scorer) {
0375 G4ScoringManager* scManager = G4ScoringManager::GetScoringManager();
0376 scManager->SetVerboseLevel(1);
0377 }
0378
0379
0380 PhysicsList* physicsList = runManagerMaster->physicsListForWorker();
0381
0382 edm::LogVerbatim("SimG4CoreApplication")
0383 << "RunManagerMTWorker::InitializeG4: start initialisation of PhysicsList for the thread " << thisID;
0384
0385
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
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
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
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
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
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
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
0586
0587
0588 edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::produce: start EventID=" << inpevt.id().event();
0589
0590 assert(m_tls != nullptr and m_tls->threadInitialized);
0591
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
0597
0598 terminateRun();
0599 }
0600 initializeRun();
0601 m_tls->currentRunNumber = inpevt.id().run();
0602 }
0603
0604
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
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
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
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
0676 G4int evtid = (G4int)inpevt.id().event();
0677 G4Event* evt = new G4Event(evtid);
0678
0679
0680
0681 resetGenParticleId(inpevt);
0682
0683 edm::Handle<edm::HepMCProduct> HepMCEvt;
0684 bool found = inpevt.getByToken(m_InToken, HepMCEvt);
0685
0686 if (found) {
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 {
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
0716 }
0717 } else {
0718
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
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 }