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