CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
RunManagerMTWorker.cc
Go to the documentation of this file.
12 
14 
23 
29 
33 
36 
38 
44 
45 #include "G4Timer.hh"
46 #include "G4Event.hh"
47 #include "G4Run.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4Threading.hh"
50 #include "G4UImanager.hh"
51 #include "G4WorkerThread.hh"
52 #include "G4WorkerRunManagerKernel.hh"
53 #include "G4StateManager.hh"
54 #include "G4TransportationManager.hh"
55 #include "G4Field.hh"
56 #include "G4FieldManager.hh"
57 
58 #include <atomic>
59 #include <memory>
60 
61 #include <thread>
62 #include <sstream>
63 #include <vector>
64 
65 static std::once_flag applyOnce;
66 
67 // from https://hypernews.cern.ch/HyperNews/CMS/get/edmFramework/3302/2.html
68 namespace {
69  std::atomic<int> thread_counter{0};
70 
71  int get_new_thread_index() { return thread_counter++; }
72 
73  bool createWatchers(const edm::ParameterSet& iP,
74  SimActivityRegistry* iReg,
75  std::vector<std::shared_ptr<SimWatcher>>& oWatchers,
76  std::vector<std::shared_ptr<SimProducer>>& oProds,
77  int threadID) {
78  std::vector<edm::ParameterSet> ws = iP.getParameter<std::vector<edm::ParameterSet>>("Watchers");
79  if (ws.empty()) {
80  return false;
81  }
82 
83  for (auto& watcher : ws) {
84  std::unique_ptr<SimWatcherMakerBase> maker(
85  SimWatcherFactory::get()->create(watcher.getParameter<std::string>("type")));
86  if (maker == nullptr) {
87  throw cms::Exception("Configuration")
88  << "RunManagerMTWorker::createWatchers: "
89  << "Unable to find the requested Watcher " << watcher.getParameter<std::string>("type");
90  } else {
91  std::shared_ptr<SimWatcher> watcherTemp;
92  std::shared_ptr<SimProducer> producerTemp;
93  maker->make(watcher, *(iReg), watcherTemp, producerTemp);
94  if (nullptr != watcherTemp) {
95  if (!watcherTemp->isMT() && 0 < threadID) {
96  throw cms::Exception("Configuration")
97  << "RunManagerMTWorker::createWatchers: "
98  << "Unable to use Watcher " << watcher.getParameter<std::string>("type") << " if number of threads > 1";
99  } else {
100  oWatchers.push_back(watcherTemp);
101  if (nullptr != producerTemp) {
102  oProds.push_back(producerTemp);
103  }
104  }
105  }
106  }
107  }
108  return (!oWatchers.empty());
109  }
110 } // namespace
111 
113  std::unique_ptr<G4RunManagerKernel> kernel; //must be deleted last
114  std::unique_ptr<RunAction> userRunAction;
115  std::unique_ptr<SimRunInterface> runInterface;
116  std::unique_ptr<SimActivityRegistry> registry;
117  std::unique_ptr<SimTrackManager> trackManager;
118  std::vector<SensitiveTkDetector*> sensTkDets;
119  std::vector<SensitiveCaloDetector*> sensCaloDets;
120  std::vector<std::shared_ptr<SimWatcher>> watchers;
121  std::vector<std::shared_ptr<SimProducer>> producers;
122  //G4Run can only be deleted if there is a G4RunManager
123  // on the thread where the G4Run is being deleted,
124  // else it causes a segmentation fault
125  G4Run* currentRun = nullptr;
126  std::unique_ptr<G4Event> currentEvent;
128  bool threadInitialized = false;
129  bool runTerminated = false;
130 
131  TLSData() {}
132 
133  ~TLSData() {}
134 };
135 
136 //This can not be a smart pointer since we must delete some of the members
137 // before leaving main() else we get a segmentation fault caused by accessing
138 // other 'singletons' after those singletons have been deleted. Instead we
139 // atempt to delete all TLS at RunManagerMTWorker destructor. If that fails for
140 // some reason, it is better to leak than cause a crash.
141 //thread_local RunManagerMTWorker::TLSData* RunManagerMTWorker::m_tls{nullptr};
142 
144  : m_generator(iConfig.getParameter<edm::ParameterSet>("Generator")),
145  m_InToken(iC.consumes<edm::HepMCProduct>(
146  iConfig.getParameter<edm::ParameterSet>("Generator").getParameter<edm::InputTag>("HepMCProductLabel"))),
147  m_theLHCTlinkToken(
148  iC.consumes<edm::LHCTransportLinkContainer>(iConfig.getParameter<edm::InputTag>("theLHCTlinkTag"))),
149  m_nonBeam(iConfig.getParameter<bool>("NonBeamEvent")),
150  m_pUseMagneticField(iConfig.getParameter<bool>("UseMagneticField")),
151  m_LHCTransport(iConfig.getParameter<bool>("LHCTransport")),
152  m_thread_index{get_new_thread_index()},
153  m_EvtMgrVerbosity(iConfig.getUntrackedParameter<int>("G4EventManagerVerbosity", 0)),
154  m_pField(iConfig.getParameter<edm::ParameterSet>("MagneticField")),
155  m_pRunAction(iConfig.getParameter<edm::ParameterSet>("RunAction")),
156  m_pEventAction(iConfig.getParameter<edm::ParameterSet>("EventAction")),
157  m_pStackingAction(iConfig.getParameter<edm::ParameterSet>("StackingAction")),
158  m_pTrackingAction(iConfig.getParameter<edm::ParameterSet>("TrackingAction")),
159  m_pSteppingAction(iConfig.getParameter<edm::ParameterSet>("SteppingAction")),
160  m_pCustomUIsession(iConfig.getUntrackedParameter<edm::ParameterSet>("CustomUIsession")),
161  m_p(iConfig) {
162  int thisID = getThreadIndex();
163  edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker for the thread " << thisID;
164 
165  // sensitive detectors
166  std::vector<std::string> onlySDs = iConfig.getParameter<std::vector<std::string>>("OnlySDs");
167  m_sdMakers = sim::sensitiveDetectorMakers(m_p, iC, onlySDs);
168 
169  // TLS and watchers
170  initializeTLS();
171  if (m_hasWatchers) {
172  for (auto& watcher : m_tls->watchers) {
173  watcher->registerConsumes(iC);
174  }
175  }
176 
177  if (m_LHCTransport) {
178  m_LHCToken = iC.consumes<edm::HepMCProduct>(edm::InputTag("LHCTransport"));
179  }
180  if (m_pUseMagneticField) {
181  m_MagField = iC.esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>();
182  }
183  edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker is constructed for the thread " << thisID;
184  unsigned int k = 0;
185  for (std::unordered_map<std::string, std::unique_ptr<SensitiveDetectorMakerBase>>::const_iterator itr =
186  m_sdMakers.begin();
187  itr != m_sdMakers.end();
188  ++itr, ++k)
189  edm::LogVerbatim("SimG4CoreApplication") << "SD[" << k << "] " << itr->first;
190 }
191 
193  m_tls = nullptr;
194  delete m_UIsession;
195 }
196 
198  for (auto& maker : m_sdMakers) {
199  maker.second->beginRun(es);
200  }
201  if (m_pUseMagneticField) {
203  }
204  if (m_hasWatchers) {
205  for (auto& watcher : m_tls->watchers) {
206  watcher->beginRun(es);
207  }
208  }
209 }
210 
212  int thisID = getThreadIndex();
213  edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::endRun for the thread " << thisID;
214  terminateRun();
215 }
216 
218  if (nullptr != m_tls) {
219  return;
220  }
221 
222  m_tls = new TLSData();
223  m_tls->registry = std::make_unique<SimActivityRegistry>();
224 
225  edm::Service<SimActivityRegistry> otherRegistry;
226  //Look for an outside SimActivityRegistry
227  // this is used by the visualization code
228  int thisID = getThreadIndex();
229  if (otherRegistry) {
230  m_tls->registry->connect(*otherRegistry);
231  if (thisID > 0) {
232  throw cms::Exception("Configuration")
233  << "RunManagerMTWorker::initializeTLS : "
234  << "SimActivityRegistry service (i.e. visualization) is not supported for more than 1 thread. "
235  << " \n If this use case is needed, RunManagerMTWorker has to be updated.";
236  }
237  }
239 }
240 
243  return;
244 
245  G4Timer timer;
246  timer.Start();
247 
248  // I guess everything initialized here should be in thread_local storage
249  initializeTLS();
250 
251  int thisID = getThreadIndex();
252  edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::initializeG4 in thread " << thisID << " is started";
253 
254  // Initialize per-thread output
255  G4Threading::G4SetThreadId(thisID);
256  G4UImanager::GetUIpointer()->SetUpForAThread(thisID);
257  const std::string& uitype = m_pCustomUIsession.getUntrackedParameter<std::string>("Type", "MessageLogger");
258  if (uitype == "MessageLogger") {
260  } else if (uitype == "MessageLoggerThreadPrefix") {
262  m_pCustomUIsession.getUntrackedParameter<std::string>("ThreadPrefix", ""), thisID);
263  } else if (uitype == "FilePerThread") {
264  m_UIsession =
266  } else {
267  throw cms::Exception("Configuration")
268  << "RunManagerMTWorker::initializeG4: Invalid value of CustomUIsession.Type '" << uitype
269  << "', valid are MessageLogger, MessageLoggerThreadPrefix, FilePerThread";
270  }
271  G4UImanager::GetUIpointer()->SetCoutDestination(m_UIsession);
272 
273  // Initialize worker part of shared resources (geometry, physics)
274  G4WorkerThread::BuildGeometryAndPhysicsVector();
275 
276  // Create worker run manager
277  m_tls->kernel.reset(G4WorkerRunManagerKernel::GetRunManagerKernel());
278  if (nullptr == m_tls->kernel) {
279  m_tls->kernel = std::make_unique<G4WorkerRunManagerKernel>();
280  }
281 
282  // Define G4 exception handler
283  double th = m_p.getParameter<double>("ThresholdForGeometryExceptions") * CLHEP::GeV;
284  G4StateManager::GetStateManager()->SetExceptionHandler(new ExceptionHandler(th, false));
285 
286  // Set the geometry for the worker, share from master
287  auto worldPV = runManagerMaster->world().GetWorldVolume();
288  m_tls->kernel->WorkerDefineWorldVolume(worldPV);
289  G4TransportationManager* tM = G4TransportationManager::GetTransportationManager();
290  tM->SetWorldForTracking(worldPV);
291 
292  // we need the track manager now
293  m_tls->trackManager = std::make_unique<SimTrackManager>();
294 
295  // setup the magnetic field
296  if (m_pUseMagneticField) {
297  const GlobalPoint g(0.f, 0.f, 0.f);
298 
299  sim::FieldBuilder fieldBuilder(m_pMagField, m_pField);
300 
301  CMSFieldManager* fieldManager = new CMSFieldManager();
302  tM->SetFieldManager(fieldManager);
303  fieldBuilder.build(fieldManager, tM->GetPropagatorInField());
304 
305  std::string fieldFile = m_p.getUntrackedParameter<std::string>("FileNameField", "");
306  if (!fieldFile.empty()) {
307  std::call_once(applyOnce, [this]() { m_dumpMF = true; });
308  if (m_dumpMF) {
309  edm::LogVerbatim("SimG4CoreApplication")
310  << "RunManagerMTWorker::InitializeG4: Dump magnetic field to file " << fieldFile;
311  DumpMagneticField(tM->GetFieldManager()->GetDetectorField(), fieldFile);
312  }
313  }
314  }
315 
316  // attach sensitive detector
317  auto sensDets = sim::attachSD(
318  m_sdMakers, es, runManagerMaster->catalog(), m_p, m_tls->trackManager.get(), *(m_tls->registry.get()));
319 
320  m_tls->sensTkDets.swap(sensDets.first);
321  m_tls->sensCaloDets.swap(sensDets.second);
322 
323  edm::LogVerbatim("SimG4CoreApplication")
324  << "RunManagerMTWorker::InitializeG4: Sensitive Detectors are built in thread " << thisID << " found "
325  << m_tls->sensTkDets.size() << " Tk type SD, and " << m_tls->sensCaloDets.size() << " Calo type SD";
326 
327  // Set the physics list for the worker, share from master
328  PhysicsList* physicsList = runManagerMaster->physicsListForWorker();
329 
330  edm::LogVerbatim("SimG4CoreApplication")
331  << "RunManagerMTWorker::InitializeG4: start initialisation of PhysicsList for the thread " << thisID;
332 
333  // Geant4 UI commands in PreInit state
334  if (!runManagerMaster->G4Commands().empty()) {
335  G4cout << "RunManagerMTWorker::InitializeG4: Requested UI commands: " << G4endl;
336  for (const std::string& command : runManagerMaster->G4Commands()) {
337  G4cout << " " << command << G4endl;
338  G4UImanager::GetUIpointer()->ApplyCommand(command);
339  }
340  }
341  G4StateManager::GetStateManager()->SetNewState(G4State_Init);
342 
343  physicsList->InitializeWorker();
344  m_tls->kernel->SetPhysics(physicsList);
345  m_tls->kernel->InitializePhysics();
346 
347  if (!m_tls->kernel->RunInitialization()) {
348  throw cms::Exception("Configuration")
349  << "RunManagerMTWorker::InitializeG4: Geant4 kernel initialization failed in thread " << thisID;
350  }
351  //tell all interesting parties that we are beginning the job
352  BeginOfJob aBeginOfJob(&es);
353  m_tls->registry->beginOfJobSignal_(&aBeginOfJob);
354 
355  G4int sv = m_p.getUntrackedParameter<int>("SteppingVerbosity", 0);
356  G4double elim = m_p.getUntrackedParameter<double>("StepVerboseThreshold", 0.1) * CLHEP::GeV;
357  std::vector<int> ve = m_p.getUntrackedParameter<std::vector<int>>("VerboseEvents");
358  std::vector<int> vn = m_p.getUntrackedParameter<std::vector<int>>("VertexNumber");
359  std::vector<int> vt = m_p.getUntrackedParameter<std::vector<int>>("VerboseTracks");
360 
361  if (sv > 0) {
362  m_sVerbose = std::make_unique<CMSSteppingVerbose>(sv, elim, ve, vn, vt);
363  }
365 
366  G4StateManager::GetStateManager()->SetNewState(G4State_Idle);
367 
368  timer.Stop();
369  edm::LogVerbatim("SimG4CoreApplication")
370  << "RunManagerMTWorker::initializeG4 done for the thread " << thisID << " " << timer;
371  m_tls->threadInitialized = true;
372 }
373 
375  m_tls->runInterface = std::make_unique<SimRunInterface>(this, false);
376  m_tls->userRunAction = std::make_unique<RunAction>(m_pRunAction, m_tls->runInterface.get(), false);
377  m_tls->userRunAction->SetMaster(false);
378  Connect(m_tls->userRunAction.get());
379 
380  G4EventManager* eventManager = m_tls->kernel->GetEventManager();
381  eventManager->SetVerboseLevel(m_EvtMgrVerbosity);
382 
383  EventAction* userEventAction =
385  Connect(userEventAction);
386  eventManager->SetUserAction(userEventAction);
387 
388  TrackingAction* userTrackingAction = new TrackingAction(userEventAction, m_pTrackingAction, m_sVerbose.get());
389  Connect(userTrackingAction);
390  eventManager->SetUserAction(userTrackingAction);
391 
392  SteppingAction* userSteppingAction =
393  new SteppingAction(userEventAction, m_pSteppingAction, m_sVerbose.get(), m_hasWatchers);
394  Connect(userSteppingAction);
395  eventManager->SetUserAction(userSteppingAction);
396 
397  eventManager->SetUserAction(new StackingAction(userTrackingAction, m_pStackingAction, m_sVerbose.get()));
398 }
399 
401  runAction->m_beginOfRunSignal.connect(m_tls->registry->beginOfRunSignal_);
402  runAction->m_endOfRunSignal.connect(m_tls->registry->endOfRunSignal_);
403 }
404 
406  eventAction->m_beginOfEventSignal.connect(m_tls->registry->beginOfEventSignal_);
407  eventAction->m_endOfEventSignal.connect(m_tls->registry->endOfEventSignal_);
408 }
409 
411  trackingAction->m_beginOfTrackSignal.connect(m_tls->registry->beginOfTrackSignal_);
412  trackingAction->m_endOfTrackSignal.connect(m_tls->registry->endOfTrackSignal_);
413 }
414 
416  steppingAction->m_g4StepSignal.connect(m_tls->registry->g4StepSignal_);
417 }
418 
420  initializeTLS();
421  return m_tls->trackManager.get();
422 }
423 std::vector<SensitiveTkDetector*>& RunManagerMTWorker::sensTkDetectors() {
424  initializeTLS();
425  return m_tls->sensTkDets;
426 }
427 std::vector<SensitiveCaloDetector*>& RunManagerMTWorker::sensCaloDetectors() {
428  initializeTLS();
429  return m_tls->sensCaloDets;
430 }
431 std::vector<std::shared_ptr<SimProducer>>& RunManagerMTWorker::producers() {
432  initializeTLS();
433  return m_tls->producers;
434 }
435 
437  int thisID = getThreadIndex();
438  edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::initializeRun " << thisID << " is started";
439  m_tls->currentRun = new G4Run();
440  G4StateManager::GetStateManager()->SetNewState(G4State_GeomClosed);
441  if (nullptr != m_tls->userRunAction) {
442  m_tls->userRunAction->BeginOfRunAction(m_tls->currentRun);
443  }
444 }
445 
447  edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::terminateRun ";
448  if (nullptr == m_tls || m_tls->runTerminated) {
449  return;
450  }
451  int thisID = getThreadIndex();
452  edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::terminateRun " << thisID << " is started";
453  if (m_tls->userRunAction) {
454  m_tls->userRunAction->EndOfRunAction(m_tls->currentRun);
455  m_tls->userRunAction.reset();
456  }
457  m_tls->currentEvent.reset();
458  m_simEvent = nullptr;
459 
460  if (m_tls->kernel) {
461  m_tls->kernel->RunTermination();
462  }
463 
464  m_tls->runTerminated = true;
465 }
466 
467 std::unique_ptr<G4SimEvent> RunManagerMTWorker::produce(const edm::Event& inpevt,
468  const edm::EventSetup& es,
469  RunManagerMT& runManagerMaster) {
470  // The initialization and begin/end run is a bit convoluted due to
471  // - Geant4 deals per-thread
472  // - OscarMTProducer deals per-stream
473  // and framework/TBB is free to schedule work in streams to the
474  // threads as it likes.
475  //
476  // We have to do the per-thread initialization, and per-thread
477  // per-run initialization here by ourselves.
478 
479  assert(m_tls != nullptr and m_tls->threadInitialized);
480  // Initialize run
481  if (inpevt.id().run() != m_tls->currentRunNumber) {
482  edm::LogVerbatim("SimG4CoreApplication")
483  << "RunManagerMTWorker::produce: RunID= " << inpevt.id().run() << " TLS RunID= " << m_tls->currentRunNumber;
484  if (m_tls->currentRunNumber != 0 && !m_tls->runTerminated) {
485  // If previous run in this thread was not terminated via endRun() call,
486  // terminate it now
487  terminateRun();
488  }
489  initializeRun();
490  m_tls->currentRunNumber = inpevt.id().run();
491  }
492  m_tls->runInterface->setRunManagerMTWorker(this); // For UserActions
493 
494  m_tls->currentEvent.reset(generateEvent(inpevt));
495 
496  auto simEvent = std::make_unique<G4SimEvent>();
497  m_simEvent = simEvent.get();
498  m_simEvent->hepEvent(m_generator.genEvent());
499  m_simEvent->weight(m_generator.eventWeight());
500  if (m_generator.genVertex() != nullptr) {
501  auto genVertex = m_generator.genVertex();
502  m_simEvent->collisionPoint(math::XYZTLorentzVectorD(genVertex->x() / CLHEP::cm,
503  genVertex->y() / CLHEP::cm,
504  genVertex->z() / CLHEP::cm,
505  genVertex->t() / CLHEP::second));
506  }
507  if (m_tls->currentEvent->GetNumberOfPrimaryVertex() == 0) {
508  throw cms::Exception("EventCorruption")
509  << "RunManagerMTWorker::produce: event " << inpevt.id().event() << " with no G4PrimaryVertices"
510  << " StreamID=" << inpevt.streamID() << " threadIndex=" << getThreadIndex();
511 
512  } else {
513  edm::LogVerbatim("SimG4CoreApplication")
514  << "RunManagerMTWorker::produce: start EventID=" << inpevt.id().event() << " StreamID=" << inpevt.streamID()
515  << " threadIndex=" << getThreadIndex() << " weight=" << m_simEvent->weight() << "; "
516  << m_tls->currentEvent->GetNumberOfPrimaryVertex() << " vertices for Geant4; generator produced "
517  << m_simEvent->nGenParts() << " particles.";
518 
519  m_tls->kernel->GetEventManager()->ProcessOneEvent(m_tls->currentEvent.get());
520  }
521 
522  //remove memory only needed during event processing
523  m_tls->currentEvent.reset();
524 
525  for (auto& sd : m_tls->sensCaloDets) {
526  sd->reset();
527  }
528  edm::LogVerbatim("SimG4CoreApplication") << "RunManagerMTWorker::produce: ended Event " << inpevt.id().event();
529  m_simEvent = nullptr;
530  return simEvent;
531 }
532 
534  if (m_tls->runTerminated) {
535  return;
536  }
537  G4Track* t = m_tls->kernel->GetEventManager()->GetTrackingManager()->GetTrack();
538  t->SetTrackStatus(fStopAndKill);
539 
540  // CMS-specific act
541  //
542  TrackingAction* uta = static_cast<TrackingAction*>(m_tls->kernel->GetEventManager()->GetUserTrackingAction());
543  uta->PostUserTrackingAction(t);
544 
545  m_tls->currentEvent->SetEventAborted();
546  m_tls->kernel->GetEventManager()->GetStackManager()->clear();
547  m_tls->kernel->GetEventManager()->GetTrackingManager()->EventAborted();
548 }
549 
550 void RunManagerMTWorker::abortRun(bool softAbort) {
551  if (!softAbort) {
552  abortEvent();
553  }
554  m_tls->currentRun = nullptr;
555  terminateRun();
556 }
557 
559  m_tls->currentEvent.reset();
560  m_simEvent = nullptr;
561 
562  // 64 bits event ID in CMSSW converted into Geant4 event ID
563  G4int evtid = (G4int)inpevt.id().event();
564  G4Event* evt = new G4Event(evtid);
565 
567  inpevt.getByToken(m_InToken, HepMCEvt);
568 
569  m_generator.setGenEvent(HepMCEvt->GetEvent());
570 
571  // required to reset the GenParticle Id for particles transported
572  // along the beam pipe
573  // to their original value for SimTrack creation
574  resetGenParticleId(inpevt);
575 
576  if (!m_nonBeam) {
577  m_generator.HepMC2G4(HepMCEvt->GetEvent(), evt);
578  if (m_LHCTransport) {
580  inpevt.getByToken(m_LHCToken, LHCMCEvt);
581  m_generator.nonCentralEvent2G4(LHCMCEvt->GetEvent(), evt);
582  }
583  } else {
584  m_generator.nonCentralEvent2G4(HepMCEvt->GetEvent(), evt);
585  }
586 
587  return evt;
588 }
589 
592  inpevt.getByToken(m_theLHCTlinkToken, theLHCTlink);
593  if (theLHCTlink.isValid()) {
594  m_tls->trackManager->setLHCTransportLink(theLHCTlink.product());
595  }
596 }
597 
598 void RunManagerMTWorker::DumpMagneticField(const G4Field* field, const std::string& file) const {
599  std::ofstream fout(file.c_str(), std::ios::out);
600  if (fout.fail()) {
601  edm::LogWarning("SimG4CoreApplication")
602  << "MTWorker::DumpMagneticField: error opening file <" << file << "> for magnetic field";
603  } else {
604  // CMS magnetic field volume
605  double rmax = 9000 * mm;
606  double zmax = 24000 * mm;
607 
608  double dr = 1 * cm;
609  double dz = 5 * cm;
610 
611  int nr = (int)(rmax / dr);
612  int nz = 2 * (int)(zmax / dz);
613 
614  double r = 0.0;
615  double z0 = -zmax;
616  double z;
617 
618  double phi = 0.0;
619  double cosf = cos(phi);
620  double sinf = sin(phi);
621 
622  double point[4] = {0.0, 0.0, 0.0, 0.0};
623  double bfield[3] = {0.0, 0.0, 0.0};
624 
625  fout << std::setprecision(6);
626  for (int i = 0; i <= nr; ++i) {
627  z = z0;
628  for (int j = 0; j <= nz; ++j) {
629  point[0] = r * cosf;
630  point[1] = r * sinf;
631  point[2] = z;
632  field->GetFieldValue(point, bfield);
633  fout << "R(mm)= " << r / mm << " phi(deg)= " << phi / degree << " Z(mm)= " << z / mm
634  << " Bz(tesla)= " << bfield[2] / tesla << " Br(tesla)= " << (bfield[0] * cosf + bfield[1] * sinf) / tesla
635  << " Bphi(tesla)= " << (bfield[0] * sinf - bfield[1] * cosf) / tesla << G4endl;
636  z += dz;
637  }
638  r += dr;
639  }
640 
641  fout.close();
642  }
643 }
RunNumber_t run() const
Definition: EventID.h:38
edm::ParameterSet m_pSteppingAction
Log< level::Info, true > LogVerbatim
EventNumber_t event() const
Definition: EventID.h:40
std::vector< std::shared_ptr< SimWatcher > > watchers
const SensitiveDetectorCatalog & catalog() const
Definition: RunManagerMT.h:72
T getUntrackedParameter(std::string const &, T const &) const
virtual const math::XYZTLorentzVector * genVertex() const
Definition: Generator.h:31
CustomUIsession * m_UIsession
void nonCentralEvent2G4(const HepMC::GenEvent *g, G4Event *e)
Definition: Generator.cc:549
const double GeV
Definition: MathUtil.h:16
std::unique_ptr< G4Event > currentEvent
virtual const double eventWeight() const
Definition: Generator.h:32
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
virtual const HepMC::GenEvent * genEvent() const
Definition: Generator.h:30
SimActivityRegistry::G4StepSignal m_g4StepSignal
SimTrackManager * GetSimTrackManager()
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
void HepMC2G4(const HepMC::GenEvent *g, G4Event *e)
Definition: Generator.cc:112
std::unique_ptr< SimRunInterface > runInterface
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::unique_ptr< CMSSteppingVerbose > m_sVerbose
SimActivityRegistry::EndOfRunSignal m_endOfRunSignal
Definition: RunAction.h:24
std::vector< SensitiveTkDetector * > sensTkDets
SimActivityRegistry::EndOfEventSignal m_endOfEventSignal
Definition: EventAction.h:47
assert(be >=bs)
void beginRun(const edm::EventSetup &)
edm::EDGetTokenT< edm::HepMCProduct > m_InToken
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
std::unique_ptr< G4RunManagerKernel > kernel
edm::EDGetTokenT< edm::HepMCProduct > m_LHCToken
auto &__restrict__ ws
bool getData(T &iHolder) const
Definition: EventSetup.h:128
U second(std::pair< T, U > const &p)
G4SimEvent * simEvent()
std::unique_ptr< SimActivityRegistry > registry
edm::ParameterSet m_pEventAction
const DDDWorld & world() const
Definition: RunManagerMT.h:70
std::unique_ptr< SimTrackManager > trackManager
edm::ParameterSet m_pRunAction
std::vector< SensitiveCaloDetector * > sensCaloDets
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::vector< SensitiveTkDetector * > & sensTkDetectors()
void resetGenParticleId(const edm::Event &inpevt)
RunManagerMTWorker(const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iC)
std::unique_ptr< G4SimEvent > produce(const edm::Event &inpevt, const edm::EventSetup &es, RunManagerMT &runManagerMaster)
std::unique_ptr< RunAction > userRunAction
int getThreadIndex() const
void hepEvent(const HepMC::GenEvent *r)
Definition: G4SimEvent.h:23
void setGenEvent(const HepMC::GenEvent *inpevt)
Definition: Generator.h:24
G4VPhysicalVolume * GetWorldVolume() const
Definition: DDDWorld.h:24
SimActivityRegistry::EndOfTrackSignal m_endOfTrackSignal
bool isValid() const
Definition: HandleBase.h:70
std::vector< std::shared_ptr< SimProducer > > producers
void abortRun(bool softAbort=false)
void connect(Observer< const T * > *iObs)
does not take ownership of memory
Definition: Signaler.h:55
void build(CMSFieldManager *fM, G4PropagatorInField *fP)
Definition: FieldBuilder.cc:31
std::pair< std::vector< SensitiveTkDetector * >, std::vector< SensitiveCaloDetector * > > attachSD(const std::unordered_map< std::string, std::unique_ptr< SensitiveDetectorMakerBase >> &, const edm::EventSetup &, const SensitiveDetectorCatalog &, edm::ParameterSet const &, const SimTrackManager *, SimActivityRegistry &reg)
static std::once_flag applyOnce
PhysicsList * physicsListForWorker() const
Definition: RunManagerMT.h:78
T const * product() const
Definition: Handle.h:70
std::vector< SensitiveCaloDetector * > & sensCaloDetectors()
edm::ParameterSet m_pTrackingAction
double sd
edm::ParameterSet m_pField
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
G4Event * generateEvent(const edm::Event &inpevt)
std::vector< std::shared_ptr< SimProducer > > & producers()
edm::EDGetTokenT< edm::LHCTransportLinkContainer > m_theLHCTlinkToken
list command
Definition: mps_check.py:25
SimActivityRegistry::BeginOfRunSignal m_beginOfRunSignal
Definition: RunAction.h:23
const std::vector< std::string > & G4Commands() const
Definition: RunManagerMT.h:74
edm::EventID id() const
Definition: EventBase.h:59
edm::ParameterSet m_pCustomUIsession
const MagneticField * m_pMagField
std::unordered_map< std::string, std::unique_ptr< SensitiveDetectorMakerBase > > sensitiveDetectorMakers(edm::ParameterSet const &, edm::ConsumesCollector, std::vector< std::string > const &chosenMakers)
StreamID streamID() const
Definition: Event.h:98
std::vector< LHCTransportLink > LHCTransportLinkContainer
static void createWatchers(const edm::ParameterSet &iP, SimActivityRegistry &iReg, std::vector< std::shared_ptr< SimWatcher >> &oWatchers, std::vector< std::shared_ptr< SimProducer >> &oProds)
void PostUserTrackingAction(const G4Track *aTrack) override
unsigned int RunNumber_t
void initializeG4(RunManagerMT *runManagerMaster, const edm::EventSetup &es)
#define get
Log< level::Warning, false > LogWarning
edm::ParameterSet m_p
edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > m_MagField
SimActivityRegistry::BeginOfEventSignal m_beginOfEventSignal
Definition: EventAction.h:46
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
void DumpMagneticField(const G4Field *, const std::string &) const
edm::ParameterSet m_pStackingAction
void Connect(RunAction *)
SimActivityRegistry::BeginOfTrackSignal m_beginOfTrackSignal