File indexing completed on 2023-10-25 10:04:21
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
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 }
0119
0120 struct RunManagerMTWorker::TLSData {
0121 std::unique_ptr<G4RunManagerKernel> kernel;
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
0131
0132
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
0145
0146
0147
0148
0149
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
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
0192 std::vector<std::string> onlySDs = p.getParameter<std::vector<std::string>>("OnlySDs");
0193 m_sdMakers = sim::sensitiveDetectorMakers(p, iC, onlySDs);
0194
0195
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
0256
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
0278 initializeTLS();
0279
0280 int thisID = getThreadIndex();
0281 edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::initializeG4 in thread " << thisID << " is started";
0282
0283
0284 G4WorkerThread::BuildGeometryAndPhysicsVector();
0285
0286
0287 m_tls->kernel.reset(G4WorkerRunManagerKernel::GetRunManagerKernel());
0288 if (nullptr == m_tls->kernel) {
0289 m_tls->kernel = std::make_unique<G4WorkerRunManagerKernel>();
0290 }
0291
0292
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
0298 auto worldPV = runManagerMaster->world().GetWorldVolume();
0299 m_tls->kernel->WorkerDefineWorldVolume(worldPV);
0300 G4TransportationManager* tM = G4TransportationManager::GetTransportationManager();
0301 tM->SetWorldForTracking(worldPV);
0302
0303
0304 m_tls->trackManager = std::make_unique<SimTrackManager>(&m_simEvent);
0305
0306
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
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
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
0357 bool scorer = m_p.getParameter<bool>("UseCommandBaseScorer");
0358 if (scorer) {
0359 G4ScoringManager* scManager = G4ScoringManager::GetScoringManager();
0360 scManager->SetVerboseLevel(1);
0361 }
0362
0363
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
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
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
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
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
0550
0551
0552
0553
0554
0555
0556
0557
0558 assert(m_tls != nullptr and m_tls->threadInitialized);
0559
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
0565
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
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
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
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
0651
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
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 }