FairRoot
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FairMCApplication.cxx
Go to the documentation of this file.
1 /********************************************************************************
2  * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
3  * *
4  * This software is distributed under the terms of the *
5  * GNU Lesser General Public Licence (LGPL) version 3, *
6  * copied verbatim in the file "LICENSE" *
7  ********************************************************************************/
8 // -------------------------------------------------------------------------
9 // ----- FairMCApplication source file -----
10 // ----- Created 06/01/04 by M. Al-Turany -----
11 // -------------------------------------------------------------------------
12 #include "FairMCApplication.h"
13 
14 #include "FairDetector.h" // for FairDetector
15 #include "FairField.h" // for FairField
16 #include "FairGenericStack.h" // for FairGenericStack
17 #include "FairGeoInterface.h" // for FairGeoInterface
18 #include "FairGeoLoader.h" // for FairGeoLoader
19 #include "FairGeoMedia.h" // for FairGeoMedia
20 #include "FairGeoMedium.h" // for FairGeoMedium
21 #include "FairIon.h" // for FairIon
22 #include "FairLogger.h" // for FairLogger, MESSAGE_ORIGIN
23 #include "FairMCEventHeader.h" // for FairMCEventHeader
24 #include "FairMesh.h" // for FairMesh
25 #include "FairModule.h" // for FairModule, etc
26 #include "FairParticle.h" // for FairParticle
27 #include "FairPrimaryGenerator.h" // for FairPrimaryGenerator
28 #include "FairRadGridManager.h" // for FairRadGridManager
29 #include "FairRadLenManager.h" // for FairRadLenManager
30 #include "FairRadMapManager.h" // for FairRadMapManager
31 #include "FairRootFileSink.h" // for CloneForWorker (in MT mode only)
32 #include "FairRootManager.h" // for FairRootManager
33 #include "FairRun.h" // for FairRun
34 #include "FairRunInfo.h" // for FairRunInfo
35 #include "FairRunSim.h" // for FairRunSim
36 #include "FairRuntimeDb.h" // for FairRuntimeDb
37 #include "FairTask.h" // for FairTask
38 #include "FairTrajFilter.h" // for FairTrajFilter
39 #include "FairVolume.h" // for FairVolume
40 
41 #include <TDatabasePDG.h> // for TDatabasePDG
42 #include <TFile.h> // for TFile
43 #include <TGeoManager.h> // for gGeoManager, TGeoManager
44 #include <TGeoMedium.h> // for TGeoMedium
45 #include <TGeoNode.h> // for TGeoNode
46 #include <TGeoPhysicalNode.h> // for TGeoPhysicalNode
47 #include <TGeoTrack.h> // for TGeoTrack
48 #include <TGeoVolume.h> // for TGeoVolume
49 #include <TH2.h> // for TH2D
50 #include <THashList.h> // for THashList
51 #include <TInterpreter.h> // for TInterpreter, gInterpreter
52 #include <TIterator.h> // for TIterator
53 #include <TList.h> // for TList, TListIter
54 #include <TObjArray.h> // for TObjArray
55 #include <TObject.h> // for TObject
56 #include <TParticlePDG.h> // for TParticlePDG
57 #include <TROOT.h> // for TROOT, gROOT
58 #include <TRefArray.h> // for TRefArray
59 #include <TSystem.h> // for TSystem, gSystem
60 #include <TVirtualMC.h> // for TVirtualMC
61 #include <TVirtualMCStack.h> // for TVirtualMCStack
62 class TParticle;
63 
64 #include <float.h> // for DBL_MAX
65 #include <mutex> // std::mutex
66 #include <stdlib.h> // for getenv, exit
67 #include <utility> // for pair
68 std::mutex mtx; // mutex for critical section
69 
70 using std::pair;
71 
72 FairMCApplication* FairMCApplication::fgMasterInstance = 0;
73 
74 //_____________________________________________________________________________
75 FairMCApplication::FairMCApplication(const char* name, const char* title, TObjArray* ModList, const char*)
76  : TVirtualMCApplication(name, title)
77  , fActiveDetectors(nullptr)
78  , fFairTaskList(nullptr)
79  , fDetectors(nullptr)
80  , fDetMap(nullptr)
81  , fModIter(nullptr)
82  , fModules(nullptr)
83  , fNoSenVolumes(0)
84  , fPythiaDecayer(kFALSE)
85  , fPythiaDecayerConfig("")
86  , fStack(nullptr)
87  , fRootManager(nullptr)
88  , fSenVolumes(nullptr)
89  , fxField(nullptr)
90  , fEvGen(nullptr)
91  , fMcVersion(-1)
92  , fTrajFilter(nullptr)
93  , fTrajAccepted(kFALSE)
94  , fUserDecay(kFALSE)
95  , fUserDecayConfig("")
96  , fDebug(kFALSE)
97  , fDisVol(nullptr)
98  , fDisDet(nullptr)
99  , fVolMap()
100  , fVolIter()
101  , fModVolMap()
102  , fModVolIter()
103  , fTrkPos(TLorentzVector(0, 0, 0, 0))
104  , fRadLength(kFALSE)
105  , fRadLenMan(nullptr)
106  , fRadMap(kFALSE)
107  , fRadMapMan(nullptr)
108  , fRadGridMan(nullptr)
109  , fEventHeader(nullptr)
110  , fMCEventHeader(nullptr)
111  , listActiveDetectors()
112  , listDetectors()
113  , fMC(nullptr)
114  , fRun(nullptr)
115  , fSaveCurrentEvent(kTRUE)
117  , fRunInfo()
118  , fGeometryIsInitialized(kFALSE)
119 {
120  // Standard Simulation constructor
121 
122  // Create an ObjArray of Modules and its iterator
123 
124  LOG(debug) << "FairMCApplication-ctor " << this;
125 
127  fModules = ModList;
128  fModIter = fModules->MakeIterator();
129  // Create and fill a list of active detectors
130  fDetectors = new TRefArray;
131  fActiveDetectors = new TRefArray();
132  fModIter->Reset();
133  FairDetector* detector;
134  TObject* obj;
135 
136  while ((obj = fModIter->Next())) {
137  detector = dynamic_cast<FairDetector*>(obj);
138  if (detector) {
139  fDetectors->Add(detector);
140  listDetectors.push_back(detector);
141  if (detector->IsActive()) {
142  fActiveDetectors->Add(detector);
143  listActiveDetectors.push_back(detector);
144  }
145  } else if (!dynamic_cast<FairModule*>(obj)) {
146  LOG(error) << "Dynamic cast fails. Object neither FairDetector nor FairModule in module list";
147  }
148  }
149 
150  // Create a Task list
151  fFairTaskList = new FairTask("Task List", 1);
152  gROOT->GetListOfBrowsables()->Add(fFairTaskList);
153  fMcVersion = -1;
154  // Initialise fTrajFilter pointer
155  fTrajFilter = nullptr;
156  fDetMap = new TRefArray(1000);
157  fDisVol = 0;
158  fDisDet = 0;
159 
160  // This ctor is used to construct the application on master
161  fgMasterInstance = this;
162 }
163 
164 //_____________________________________________________________________________
166  : TVirtualMCApplication(rhs.GetName(), rhs.GetTitle())
167  , fActiveDetectors(nullptr)
168  , fFairTaskList(nullptr)
169  , fDetectors(nullptr)
170  , fDetMap(nullptr)
171  , fModIter(nullptr)
172  , fModules(nullptr)
173  , fNoSenVolumes(0)
174  , fPythiaDecayer(kFALSE)
175  , fPythiaDecayerConfig(rhs.fPythiaDecayerConfig)
176  , fStack(nullptr)
177  , fRootManager(nullptr)
178  , fSenVolumes(nullptr)
179  , fxField(rhs.fxField)
180  , fEvGen(nullptr)
181  , fMcVersion(rhs.fMcVersion)
182  , fTrajFilter(nullptr)
183  , fTrajAccepted(kFALSE)
184  , fUserDecay(kFALSE)
185  , fUserDecayConfig(rhs.fUserDecayConfig)
186  , fDebug(rhs.fDebug)
187  , fDisVol(nullptr)
188  , fDisDet(nullptr)
189  , fVolMap()
190  , fVolIter()
191  , fModVolMap()
192  , fModVolIter()
193  , fTrkPos(rhs.fTrkPos)
194  , fRadLength(kFALSE)
195  , fRadLenMan(nullptr)
196  , fRadMap(kFALSE)
197  , fRadMapMan(nullptr)
198  , fRadGridMan(nullptr)
199  , fEventHeader(nullptr)
200  , fMCEventHeader(nullptr)
201  , listActiveDetectors()
202  , listDetectors()
203  , fMC(nullptr)
204  , fRun(nullptr)
205  , fSaveCurrentEvent(kTRUE)
207  , fRunInfo()
208  , fGeometryIsInitialized(kFALSE)
209 {
210  // Copy constructor
211  // Do not create Root manager
212 
213  LOG(debug) << "FairMCApplication-copy-ctor " << this;
214 
215  // Create an ObjArray of Modules and its iterator
216  fModules = new TObjArray();
217  fModIter = fModules->MakeIterator();
218  // Clone modules
219  TObject* obj;
220  rhs.fModIter->Reset();
221  while ((obj = rhs.fModIter->Next())) {
222  LOG(debug) << "cloning " << (static_cast<FairModule*>(obj))->GetName();
223  fModules->Add(static_cast<FairModule*>(obj)->CloneModule());
224  }
225 
226  // Create and fill a list of active detectors
227  fDetectors = new TRefArray;
228  fActiveDetectors = new TRefArray();
229  fModIter->Reset();
230  FairDetector* detector;
231  while ((obj = fModIter->Next())) {
232  if (obj->InheritsFrom("FairDetector")) {
233  detector = dynamic_cast<FairDetector*>(obj);
234  if (detector) {
235  fDetectors->Add(detector);
236  listDetectors.push_back(detector);
237  if (detector->IsActive()) {
238  fActiveDetectors->Add(detector);
239  listActiveDetectors.push_back(detector);
240  }
241  } else {
242  LOG(error) << "Dynamic cast fails.";
243  }
244  }
245  }
246 
247  // Clone stack
248  fStack = rhs.fStack->CloneStack();
249 
250  // Create a Task list
251  // Let's try without it
252  // fFairTaskList= new FairTask("Task List", 1);
253  // gROOT->GetListOfBrowsables()->Add(fFairTaskList);
254 
255  fDetMap = new TRefArray(1000);
256 }
257 //_____________________________________________________________________________
259  : TVirtualMCApplication()
260  , fActiveDetectors(0)
261  , fFairTaskList(0)
262  , fDetectors(0)
263  , fDetMap(0)
264  , fModIter(0)
265  , fModules(0)
266  , fNoSenVolumes(0)
267  , fPythiaDecayer(kFALSE)
268  , fPythiaDecayerConfig("")
269  , fStack(0)
270  , fRootManager(0)
271  , fSenVolumes(0)
272  , fxField(0)
273  , fEvGen(0)
274  , fMcVersion(-1)
275  , fTrajFilter(nullptr)
276  , fTrajAccepted(kFALSE)
277  , fUserDecay(kFALSE)
278  , fUserDecayConfig("")
279  , fDebug(kFALSE)
280  , fDisVol(0)
281  , fDisDet(0)
282  , fVolMap()
283  , fVolIter()
284  , fModVolMap()
285  , fModVolIter()
286  , fTrkPos(TLorentzVector(0, 0, 0, 0))
287  , fRadLength(kFALSE)
288  , fRadLenMan(nullptr)
289  , fRadMap(kFALSE)
290  , fRadMapMan(nullptr)
291  , fRadGridMan(nullptr)
292  , fEventHeader(nullptr)
293  , fMCEventHeader(nullptr)
294  , listActiveDetectors()
295  , listDetectors()
296  , fMC(nullptr)
297  , fRun(nullptr)
298  , fSaveCurrentEvent(kTRUE)
300  , fRunInfo()
301  , fGeometryIsInitialized(kFALSE)
302 {
303  // Default constructor
304 }
305 //_____________________________________________________________________________
307 {
308  // Destructor
309  // LOG(debug3) << "Enter Destructor of FairMCApplication";
310  delete fStack;
311  delete fActiveDetectors; // don't do fActiveDetectors->Delete() here
312  // the modules are already deleted in FairRunSim
313  delete fDetectors;
314  delete fModIter;
315  // LOG(debug3) << "Leave Destructor of FairMCApplication";
316  delete fMC;
317 }
318 
319 //_____________________________________________________________________________
320 FairMCApplication& FairMCApplication::operator=(const FairMCApplication& rhs)
321 {
322  // Assignment operator
323 
324  // check assignment to self
325  if (this != &rhs) {
326 
327  // base class assignment
328  TVirtualMCApplication::operator=(rhs);
329 
330  fActiveDetectors = nullptr;
331  fFairTaskList = nullptr;
332  fDetectors = nullptr;
333  fDetMap = nullptr;
334  fModIter = nullptr;
335  fModules = nullptr;
336  fNoSenVolumes = 0;
337  fPythiaDecayer = kFALSE;
339  fStack = nullptr;
340  fRootManager = nullptr;
341  fSenVolumes = nullptr;
342  fxField = rhs.fxField;
343  fEvGen = nullptr;
344  fMcVersion = rhs.fMcVersion;
345  fTrajFilter = nullptr;
346  fTrajAccepted = kFALSE;
347  fUserDecay = kFALSE;
349  fDebug = rhs.fDebug;
350  fDisVol = nullptr;
351  fDisDet = nullptr;
352  fTrkPos = rhs.fTrkPos;
353  fRadLength = kFALSE;
354  fRadLenMan = nullptr;
355  fRadMap = kFALSE;
356  fRadMapMan = nullptr;
357  fRadGridMan = nullptr;
358  fEventHeader = nullptr;
359  fMCEventHeader = nullptr;
360  fGeometryIsInitialized = kFALSE;
361 
362  // Do not create Root manager
363 
364  // Create an ObjArray of Modules and its iterator
365  fModules = new TObjArray();
366  fModIter = fModules->MakeIterator();
367  // Clone modules
368  TObject* obj;
369  while ((obj = rhs.fModIter->Next())) {
370  fModules->Add(static_cast<FairModule*>(obj)->CloneModule());
371  }
372 
373  // Create and fill a list of active detectors
374  fDetectors = new TRefArray;
375  fActiveDetectors = new TRefArray();
376  fModIter->Reset();
377  FairDetector* detector;
378  while ((obj = fModIter->Next())) {
379  if (obj->InheritsFrom("FairDetector")) {
380  detector = dynamic_cast<FairDetector*>(obj);
381  if (detector) {
382  fDetectors->Add(detector);
383  listDetectors.push_back(detector);
384  if (detector->IsActive()) {
385  fActiveDetectors->Add(detector);
386  listActiveDetectors.push_back(detector);
387  }
388  } else {
389  LOG(error) << "Dynamic cast fails.";
390  }
391  }
392  }
393 
394  // Clone stack
395  fStack = rhs.fStack->CloneStack();
396 
397  // Create a Task list
398  // Let's try without it
399  // fFairTaskList= new FairTask("Task List", 1);
400  // gROOT->GetListOfBrowsables()->Add(fFairTaskList);
401 
402  fDetMap = new TRefArray(1000);
403 
404  fState = rhs.fState;
405  }
406 
407  return *this;
408 }
409 
410 //_____________________________________________________________________________
411 void FairMCApplication::InitMC(const char*, const char*)
412 {
413  // Initialize MC.
414  // ---
415  // This methode is called from FairRunSim::SetMCConfig which excucute first the gconfig
416  // macro that creates the MC instance (G3 or G4)
417 
418  fMC = TVirtualMC::GetMC();
419 
420  if (fMC == 0) {
421  LOG(fatal) << "No MC engine defined";
422  }
423 
424  fStack = dynamic_cast<FairGenericStack*>(fMC->GetStack());
425  if (fStack == nullptr) {
426  LOG(fatal) << "No Stack defined.";
427  }
428  fMC->SetMagField(fxField);
429 
431  // fRootManager->SetDebug(true);
432 
433  fMC->Init();
434  fMC->BuildPhysics();
435  TString MCName = fMC->GetName();
436  if (MCName == "TGeant3" || MCName == "TGeant3TGeo") {
437  fMcVersion = 0;
438  } else if (MCName == "TGeant4") {
439  fMcVersion = 1;
440  } else if (MCName == "TFluka") {
441  fMcVersion = 2;
442  } else {
443  fMcVersion = 3; // Geane
444  }
446 
447  LOG(info) << "Monte Carlo Engine Initialisation with: " << MCName.Data();
448 }
449 
450 //_____________________________________________________________________________
451 void FairMCApplication::RunMC(Int_t nofEvents)
452 {
453  // Reset the time for FairRunInfo. Otherwise the time of the
454  // first event will include the time needed for initilization.
455  fRunInfo.Reset();
456 
459 
460  // MC run.
461  fMC->ProcessRun(nofEvents);
462  // finish run
463  FinishRun();
464 }
465 
466 //____________________________________________________________________________
468 {
469  // Finish MC run.
470  // ---
471 
472  LOG(debug) << "FairMCMCApplication::FinishRun() start";
473 
474  for (auto detectorPtr : listActiveDetectors) {
475  detectorPtr->FinishRun();
476  }
477 
478  // TO DO: clone this in MT mode ?
479  if (fFairTaskList) {
481  }
482  // fRootManager->Fill();
483 
485  // FairMCEventHeader* header = gen->GetEvent();
486  Int_t nprimary = gen->GetTotPrimary();
487  TObjArray* meshlist = nullptr;
488 
489  // Only in sequential mode
490  if (fRadGridMan) {
491 
492  meshlist = fRadGridMan->GetMeshList();
493 
494  TH2D* tid = nullptr;
495  TH2D* flu = nullptr;
496  TH2D* seu = nullptr;
497 
498  LOG(info) << "=======================================================";
499  LOG(info) << " Dosimetry histos saving in \"" << fRadGridMan->GetOutputFileName() << "\"";
500  LOG(info) << "=======================================================";
501 
502  TFile* radGridFile = TFile::Open(fRadGridMan->GetOutputFileName(), "recreate");
503  radGridFile->mkdir("Dosimetry");
504  radGridFile->cd("Dosimetry");
505 
506  for (Int_t i = 0; i < meshlist->GetEntriesFast(); i++) {
507  FairMesh* aMesh = dynamic_cast<FairMesh*>(meshlist->At(i));
508  if (aMesh) {
509  aMesh->Scale(1. / nprimary);
510  tid = aMesh->GetMeshTid();
511  flu = aMesh->GetMeshFlu();
512  seu = aMesh->GetMeshSEU();
513  tid->Write();
514  flu->Write();
515  seu->Write();
516  }
517  }
518  radGridFile->Write();
519  radGridFile->Close();
520  delete radGridFile;
521  }
522 
523  // Save histograms with memory and runtime information in the output file
524  if (FairRunSim::Instance()->IsRunInfoGenerated()) {
525  fRunInfo.WriteInfo();
526  }
527 
528  if (!fRadGridMan && fRootManager) {
529  fRootManager->Write();
531  }
532 
533  if (!fMC->IsMT()) {
534  UndoGeometryModifications();
535  }
536 
537  LOG(debug) << "Done FairMCMCApplication::FinishRun()";
538 }
539 
540 //_____________________________________________________________________________
542 {
543  // User actions at beginning of event
544  // ---
545 
546  for (auto detectorPtr : listActiveDetectors) {
547  detectorPtr->BeginEvent();
548  }
549 }
550 
551 //_____________________________________________________________________________
553 {
554  // User actions at beginning of a primary track
555  // ---
556  for (auto detectorPtr : listActiveDetectors) {
557  detectorPtr->BeginPrimary();
558  }
559 }
560 
561 //_____________________________________________________________________________
563 {
564  // User actions at beginning of each track
565  // ---
566 
567  for (auto detectorPtr : listActiveDetectors) {
568  detectorPtr->PreTrack();
569  }
570 
571  fTrajAccepted = kFALSE;
572  if (nullptr != fTrajFilter) {
573  // Get the pointer to current track
574  TParticle* particle = fStack->GetCurrentTrack();
575  // LOG(debug) << " FairMCApplication::PreTrack(): " << particle;
576  // Apply cuts
577  fTrajAccepted = fTrajFilter->IsAccepted(particle);
578  if (fTrajAccepted) {
579  // Add trajectory to geo manager
580  // Int_t trackId = fStack->GetCurrentTrackNumber();
581  TGeoTrack* fTrack = fTrajFilter->CheckAddTrack(fMC->GetStack()->GetCurrentTrackNumber(), particle);
582  // TLorentzVector pos;
583  fMC->TrackPosition(fTrkPos);
584  fTrack->AddPoint(fTrkPos.X(), fTrkPos.Y(), fTrkPos.Z(), fTrkPos.T());
585  }
586  }
587 }
588 
589 //_____________________________________________________________________________
590 TVirtualMCApplication* FairMCApplication::CloneForWorker() const
591 {
592  mtx.lock();
593  LOG(info) << "FairMCApplication::CloneForWorker ";
594 
595  // Create new FairRunSim object on worker
596  // and pass some data from master FairRunSim object
597  FairRunSim* workerRun = new FairRunSim(kFALSE);
598  workerRun->SetName(fRun->GetName()); // Transport engine
599  workerRun->SetSink(fRun->GetSink()->CloneSink());
600 
601  // Trajectories filter is created explicitly as we do not call
602  // FairRunSim::Init on workers
603  if (fRun->GetStoreTraj()) {
604  new FairTrajFilter();
605  }
606 
607  // Create new FairMCApplication object on worker
608  FairMCApplication* workerApplication = new FairMCApplication(*this);
609  workerApplication->SetGenerator(fEvGen->ClonePrimaryGenerator());
610 
611  mtx.unlock();
612 
613  return workerApplication;
614 }
615 
616 //_____________________________________________________________________________
618 {
619  // Create Root manager
621 
622  LOG(info) << "FairMCApplication::InitForWorker " << fRootManager->GetInstanceId() << " " << this;
623 
624  // Set FairRunSim worker(just for consistency, not needed on worker)
626 
627  // Generate per-thread file name
628  // and create a new sink on worker
629  // moved to CloneForWorker and specific sink->CloneSink();
630 
632 
633  // Cache thread-local gMC
634  fMC = gMC;
635 
636  // Set data to MC
637  fMC->SetStack(fStack);
638  fMC->SetMagField(fxField);
639 
640  LOG(info) << "Monte Carlo Engine Worker Initialisation with: " << fMC->GetName();
641 }
642 
643 //_____________________________________________________________________________
645 {
646  LOG(debug) << "FairMCApplication::FinishRunOnWorker: ";
647 
648  FinishRun();
649 }
650 
651 //_____________________________________________________________________________
653 {
654  // User actions at each step
655  // ---
656 
657  // Work around for Fluka VMC, which does not call
658  // MCApplication::PreTrack()
659  static Int_t TrackId = 0;
660  if (fMcVersion == 2 && fMC->GetStack()->GetCurrentTrackNumber() != TrackId) {
661  PreTrack();
662  TrackId = fMC->GetStack()->GetCurrentTrackNumber();
663  }
664 
665  // Check if the volume with id is in the volume multimap.
666  // If it is not in the map the volume is not a sensitive volume
667  // and we do not call nay of our ProcessHits functions.
668 
669  // If the volume is in the multimap, check in second step if the current
670  // copy is alredy inside the multimap.
671  // If the volume is not in the multimap add the copy of the volume to the
672  // multimap.
673  // In any case call the ProcessHits function for this specific detector.
674  Int_t copyNo;
675  Int_t id = fMC->CurrentVolID(copyNo);
676  Bool_t InMap = kFALSE;
677  fDisVol = 0;
678  fDisDet = 0;
679  Int_t fCopyNo = 0;
680  fVolIter = fVolMap.find(id);
681 
682  if (fVolIter != fVolMap.end()) {
683 
684  // Call Process hits for FairVolume with this id, copyNo
685  do {
686  fDisVol = fVolIter->second;
687  fCopyNo = fDisVol->getCopyNo();
688  if (copyNo == fCopyNo) {
690  if (fDisDet) {
692  }
693  InMap = kTRUE;
694  break;
695  }
696  ++fVolIter;
697  } while (fVolIter != fVolMap.upper_bound(id));
698 
699  // if (fDisVol && !InMap) { // fDisVolume is set previously, no check needed
700 
701  // Create new FairVolume with this id, copyNo.
702  // Use the FairVolume with the same id found in the map to get
703  // the link to the detector.
704  // Seems that this never happens (?)
705  if (!InMap) {
706  // cout << "Volume not in map; fDisVol ? " << fDisVol << endl
707  FairVolume* fNewV = new FairVolume(fMC->CurrentVolName(), id);
708  fNewV->setMCid(id);
709  fNewV->setModId(fDisVol->getModId());
710  fNewV->SetModule(fDisVol->GetModule());
711  fNewV->setCopyNo(copyNo);
712  fVolMap.insert(pair<Int_t, FairVolume*>(id, fNewV));
714 
715  // LOG(info) << "FairMCApplication::Stepping: new fair volume"
716  // << id << " " << copyNo << " " << fDisDet;
717  if (fDisDet) {
718  fDisDet->ProcessHits(fNewV);
719  }
720  }
721  }
722 
723  // If information about the tracks should be stored the information as to be
724  // stored for any step.
725  // Information about each single step has also to be stored for the other
726  // special run modes of the simulation which are used to store information
727  // about
728  // 1.) Radiation length in each volume
729  // 2.) Energy deposition in each volume
730  // 3.) Fluence of particles through a defined plane which can be anywhere
731  // in the geometry. This plane has not to be correlated with any real
732  // volume
733  if (fTrajAccepted) {
734  if (fMC->TrackStep() > fTrajFilter->GetStepSizeCut()) {
735  fMC->TrackPosition(fTrkPos);
736  fTrajFilter->GetCurrentTrk()->AddPoint(fTrkPos.X(), fTrkPos.Y(), fTrkPos.Z(), fTrkPos.T());
737  }
738  }
739  if (fRadLenMan) {
740  id = fMC->CurrentVolID(copyNo);
741  fModVolIter = fgMasterInstance->fModVolMap.find(id);
742  fRadLenMan->AddPoint(fModVolIter->second);
743  }
744  if (fRadMapMan) {
745  id = fMC->CurrentVolID(copyNo);
746  fModVolIter = fgMasterInstance->fModVolMap.find(id);
747  fRadMapMan->AddPoint(fModVolIter->second);
748  }
749  if (fRadGridMan) {
751  }
752 }
753 
754 //_____________________________________________________________________________
756 {
757  // User actions after finishing of each track
758  // ---
759  for (auto detectorPtr : listActiveDetectors) {
760  detectorPtr->PostTrack();
761  }
762 }
763 
764 //_____________________________________________________________________________
766 {
767  // User actions after finishing of a primary track
768  // ---
769  for (auto detectorPtr : listActiveDetectors) {
770  detectorPtr->FinishPrimary();
771  }
772 
774 }
775 
776 //_____________________________________________________________________________
778 {
779  FinishEvent();
780  FinishRun();
781  if (fRootManager) {
782  fRootManager->Write();
784  }
785  LOG(warn) << "StopRun() exiting not safetly oopps !!!@@@!!!";
786  exit(0);
787 }
788 
789 //_____________________________________________________________________________
791 {
792  if (fMC)
793  fMC->StopRun();
794 }
795 
796 //_____________________________________________________________________________
798 {
799  // User actions after finishing of an event
800  // ---
801  LOG(debug) << "[" << fRootManager->GetInstanceId()
802  << " FairMCMCApplication::FinishEvent: " << fMCEventHeader->GetEventID() << " (MC "
803  << gMC->CurrentEvent() << ")";
804  if (gMC->IsMT()
805  && fRun->GetSink()->GetSinkType() == kONLINESINK) { // fix the rare case when running G4 multithreaded on MQ
806  fMCEventHeader->SetEventID(gMC->CurrentEvent() + 1);
807  }
808 
809  // --> Fill the stack output array
811  // --> Update track indizes in MCTracks and MCPoints
813  // --> Screen output of stack
814  // fStack->Print();
815 
816  if (fFairTaskList) {
819  }
820 
821  for (auto detectorPtr : listActiveDetectors) {
822  detectorPtr->FinishEvent();
823  }
824 
826  fRootManager->Fill();
827  } else {
828  fSaveCurrentEvent = kTRUE;
829  }
830 
831  for (auto detectorPtr : listActiveDetectors) {
832  detectorPtr->EndOfEvent();
833  }
834 
835  fStack->Reset();
836  if (nullptr != fTrajFilter) {
837  fTrajFilter->Reset();
838  // TObjArray* fListOfTracks=gGeoManager->GetListOfTracks();
839  // fListOfTracks->Delete();
840  gGeoManager->GetListOfTracks()->Delete();
841  }
842  if (nullptr != fRadLenMan) {
843  fRadLenMan->Reset();
844  }
845  if (nullptr != fRadMapMan) {
846  fRadMapMan->Reset();
847  }
848 
849  // Store information about runtime for one event and memory consuption
850  // for later usage.
851  if ((FairRunSim::Instance()->IsRunInfoGenerated()) && !gMC->IsMT()) {
852  fRunInfo.StoreInfo();
853  }
854 }
855 //_____________________________________________________________________________
857 {
858  // No limit
859  // ---
860  return DBL_MAX;
861 }
862 //_____________________________________________________________________________
864 {
865  // No limit
866  // ---
867  return DBL_MAX;
868 }
869 //_____________________________________________________________________________
871 //_____________________________________________________________________________
873 {
875  FairGeoInterface* GeoInterface = loader->getGeoInterface();
876  FairGeoMedia* media = GeoInterface->getMedia();
877  TList* MediaList = media->getListOfMedia();
878  TListIter iter(MediaList);
879  FairGeoMedium* medium;
880  Int_t NK = 0;
881  Double_t p[4];
882  while ((medium = dynamic_cast<FairGeoMedium*>(iter.Next()))) {
883  NK = medium->getNpckov();
884  if (NK > 0) {
885  Int_t Mid = 0;
886  TGeoMedium* Med = 0;
887  if (gGeoManager && (Med = gGeoManager->GetMedium(medium->GetName()))) {
888  Mid = Med->GetId();
889  } else {
890  Mid = medium->getMediumIndex();
891  if (Mid <= 0) {
892  continue;
893  }
894  }
895  Double_t ppckov[NK], absco[NK], effic[NK], rindex[NK];
896  for (Int_t i = 0; i < NK; i++) {
897  medium->getCerenkovPar(i, p);
898  ppckov[i] = p[0] * 1E-9;
899  absco[i] = p[1];
900  effic[i] = p[2];
901  rindex[i] = p[3];
902  }
903  TVirtualMC::GetMC()->SetCerenkov(Mid, NK, ppckov, absco, effic, rindex);
904  }
905  }
906  fModIter->Reset();
907  FairModule* Mod = nullptr;
908  while ((Mod = dynamic_cast<FairModule*>(fModIter->Next()))) {
909  Mod->ConstructOpGeometry();
910  }
911 }
912 
913 //_____________________________________________________________________________
915 {
916  // Construct geometry and also fill following member data:
917  // - fModVolMap: (volId,moduleId)
918  // - fSenVolumes: list of sensitive volumes
919  if (!gGeoManager) {
920  LOG(fatal) << "gGeoManager not initialized at FairMCApplication::ConstructGeometry\n";
921  }
922 
924 
925  fModIter->Reset();
926  FairModule* Mod = nullptr;
927  Int_t NoOfVolumes = 0;
928  Int_t NoOfVolumesBefore = 0;
929  Int_t ModId = 0;
930 
931  TObjArray* tgeovolumelist = gGeoManager->GetListOfVolumes();
932 
933  while ((Mod = dynamic_cast<FairModule*>(fModIter->Next()))) {
934  NoOfVolumesBefore = tgeovolumelist->GetEntriesFast();
935  Mod->InitParContainers();
936  Mod->ConstructGeometry();
937  ModId = Mod->GetModId();
938  NoOfVolumes = tgeovolumelist->GetEntriesFast();
939  for (Int_t n = NoOfVolumesBefore; n < NoOfVolumes; n++) {
940  TGeoVolume* v = (TGeoVolume*)tgeovolumelist->At(n);
941  fModVolMap.insert(pair<Int_t, Int_t>(v->GetNumber(), ModId));
942  }
943  }
944 
945  // LOG(debug) << "FairMCApplication::ConstructGeometry() : Now closing the geometry";
946  int NoOfVolumesBeforeClose = tgeovolumelist->GetEntries();
947  gGeoManager->CloseGeometry(); // close geometry
948  int NoOfVolumesAfterClose = tgeovolumelist->GetEntries();
949 
950  // Check if CloseGeometry has modified the volume list which might happen for
951  // runtime shapes (parametrizations,etc). If this is the case, our lookup structures
952  // built above are likely out of date. Issue at least a warning here.
953  if (NoOfVolumesBeforeClose != NoOfVolumesAfterClose) {
954  LOG(error) << "TGeoManager::CloseGeometry() modified the volume list from " << NoOfVolumesBeforeClose << " to "
955  << NoOfVolumesAfterClose << "\n"
956  << "This almost certainly means inconsistent lookup structures used in simulation/stepping.\n";
957  }
958 
959  if (fRun->IsImportTGeoToVMC()) {
960  TVirtualMC::GetMC()->SetRootGeometry(); // notify VMC about Root geometry
961  LOG(info) << "TGeometry will be imported to VMC"
962  << "\n";
963  } else {
964  LOG(info) << "TGeometry will not be imported to VMC"
965  << "\n";
966  }
967  Int_t Counter = 0;
968  TDatabasePDG* pdgDatabase = TDatabasePDG::Instance();
969  const THashList* list = pdgDatabase->ParticleList();
970  if (list == 0)
971  pdgDatabase->ReadPDGTable();
972  list = pdgDatabase->ParticleList();
973  if (list != 0) {
974  TIterator* particleIter = list->MakeIterator();
975  TParticlePDG* Particle = 0;
976  while ((Particle = dynamic_cast<TParticlePDG*>(particleIter->Next())) && (Counter <= 256)) {
977  TString Name = gGeoManager->GetPdgName(Particle->PdgCode());
978  // LOG(info) << Counter <<" : Particle name: "<< Name.Data() << " PDG " << Particle->PdgCode();
979  if (Name == "XXX")
980  gGeoManager->SetPdgName(Particle->PdgCode(), Particle->GetName());
981  Counter++;
982  }
983  delete particleIter;
984  }
985  fModIter->Reset();
986  while ((Mod = dynamic_cast<FairModule*>(fModIter->Next()))) {
988  }
989 
990  // dont use this here anymore, use FairMCApplication::MisalignGeometry
991  // fRun->AlignGeometry()
992 
993  gGeoManager->RefreshPhysicalNodes(kFALSE);
994 
996 }
997 
998 // ____________________________________________________________________________
1000 {
1001  // call this only here
1002  fRun->AlignGeometry();
1003  return true;
1004 }
1005 
1006 //_____________________________________________________________________________
1008 {
1010 
1011  LOG(info) << "FairMCApplication::InitGeometry: " << fRootManager->GetInstanceId();
1012 
1013  // ToBeDone
1014  // Recently the InitGeometry is called twice from the G4VMC, This is a work around tell the problem get fixed in
1015  // G4VMC
1016  if (fGeometryIsInitialized)
1017  return;
1018 
1020 
1022  FairVolume* fv = 0;
1023  Int_t id = 0;
1024  fModIter->Reset();
1025 
1026  // Register stack
1027  if (fEvGen && fStack && fRootManager) {
1028  fStack->Register();
1029  } else {
1030  LOG(warn) << "Stack is not registered ";
1031  }
1032 
1034  // initialize and register FairDetector objects in addition
1035  // Note: listActiveDetectors or fActiveDetectors not used to include passive modules in the same loop.
1036  FairModule* module;
1037  FairDetector* detector;
1038  TObject* obj;
1039  fModIter->Reset();
1040  while ((obj = fModIter->Next())) {
1041  detector = dynamic_cast<FairDetector*>(obj);
1042  module = dynamic_cast<FairModule*>(obj);
1043  if (module) {
1044  module->SetSpecialPhysicsCuts();
1045  }
1046  if (detector) {
1047  // check whether detector is active
1048  if (detector->IsActive()) {
1049  detector->Initialize();
1050  detector->Register();
1051  }
1052  }
1053  }
1054  fModIter->Reset();
1055 
1058  if (fFairTaskList) {
1061  InitTasks();
1062  }
1063 
1064  // store the EventHeader Info
1065  // Get and register EventHeader
1066  UInt_t runId = FairRunSim::Instance()->GetRunId();
1067 
1068  LOG(info) << "Simulation RunID: " << runId;
1069 
1070  // Get and register the MCEventHeader
1072  fMCEventHeader->SetRunID(runId);
1073  if (fRootManager) {
1075  }
1076 
1077  if (fEvGen) {
1079  }
1081  if (nullptr != fTrajFilter) {
1082  fTrajFilter->Init();
1083  }
1084  if (nullptr != fRadLenMan) {
1085  fRadLenMan->Init();
1086  }
1087  if (nullptr != fRadMapMan) {
1088  fRadMapMan->Init();
1089  }
1090  if (nullptr != fRadGridMan) {
1091  fRadGridMan->Init();
1092  }
1093 
1095  if (fRootManager) {
1097  }
1098 
1099  // Get static thread local svList
1101  if (fSenVolumes) {
1102  fNoSenVolumes = fSenVolumes->GetEntries();
1103  }
1104 
1105  // Fill sensitive volumes in fVolMap
1106  for (Int_t i = 0; i < fNoSenVolumes; i++) {
1107 
1108  fv = dynamic_cast<FairVolume*>(fSenVolumes->At(i));
1109  if (!fv) {
1110  LOG(error) << "No FairVolume in fSenVolumes at position " << i;
1111  continue;
1112  }
1113  id = fv->getMCid();
1114  if (fv->getGeoNode() == 0) {
1115  TGeoNode* fN = 0;
1116  TGeoVolume* v = gGeoManager->GetVolume(fv->GetName());
1117  TObjArray* fNs = 0;
1118  if (v) {
1119  fNs = v->GetNodes();
1120  }
1121  if (fNs) {
1122  for (Int_t k = 0; k < fNs->GetEntriesFast(); k++) {
1123  fN = dynamic_cast<TGeoNode*>(fNs->At(k));
1124  if (!fN) {
1125  LOG(error) << "No TGeoNode in fNs at position " << k;
1126  continue;
1127  }
1128  FairVolume* fNewV = new FairVolume(fv->GetName(), id);
1129  fNewV->setModId(fv->getModId());
1130  fNewV->SetModule(fv->GetModule());
1131  fNewV->setCopyNo(fN->GetNumber());
1132  fNewV->setMCid(id);
1133  fVolMap.insert(pair<Int_t, FairVolume*>(id, fNewV));
1134  }
1135  } else {
1136  FairVolume* fNewV = new FairVolume(fv->GetName(), id);
1137  fNewV->setModId(fv->getModId());
1138  fNewV->SetModule(fv->GetModule());
1139  fNewV->setCopyNo(1);
1140  fNewV->setMCid(id);
1141  fVolMap.insert(pair<Int_t, FairVolume*>(id, fNewV));
1142  }
1143  } else {
1144  fVolMap.insert(pair<Int_t, FairVolume*>(id, fv));
1145  }
1146  } // end off loop Fill sensitive volumes
1147 
1148  fGeometryIsInitialized = kTRUE;
1149 
1151 }
1152 
1153 //_____________________________________________________________________________
1155 {
1156  // Fill the user stack (derived from TVirtualMCStack) with primary particles.
1157  // ---
1158  LOG(debug) << "FairMCApplication::GeneratePrimaries: " << fEvGen;
1159 
1160  if (fEvGen) {
1161  // LOG(debug) << "FairMCApplication::GeneratePrimaries()";
1162  if (!fEvGen->GenerateEvent(fStack)) {
1163  StopRun();
1164  }
1165  }
1166 }
1167 //_____________________________________________________________________________
1169 {
1170  return dynamic_cast<FairDetector*>(fModules->FindObject(DetName));
1171 }
1172 //_____________________________________________________________________________
1174 {
1175  TDatabasePDG* pdgDatabase = TDatabasePDG::Instance();
1176  TObjArray* NewIons = fRun->GetUserDefIons();
1177  TIterator* Iter = NewIons->MakeIterator();
1178  Iter->Reset();
1179  TObject* obj = 0;
1180  FairIon* ion = 0;
1181  while ((obj = Iter->Next())) {
1182  ion = dynamic_cast<FairIon*>(obj);
1183  if (ion) {
1184  // Check if an ion with the calculated pdg code already exists in the
1185  // TDatabasePDG.
1186  // If the ion already exist don't create a new one because this fails for
1187  // Geant4. Instead modify the FairIon to use the already existing ion
1188  // from the TDatabasePDG.
1189  // The problem occured for example for Alphas which exist already.
1190  Int_t ionPdg = GetIonPdg(ion->GetZ(), ion->GetA());
1191  if (!pdgDatabase->GetParticle(ionPdg)) {
1192  fMC->DefineIon(
1193  ion->GetName(), ion->GetZ(), ion->GetA(), ion->GetQ(), ion->GetExcEnergy(), ion->GetMass());
1194 
1195  } else {
1196  ion->SetName(pdgDatabase->GetParticle(ionPdg)->GetName());
1197  }
1198  // Add Ion to gGeoManager visualization
1199  if (gGeoManager) {
1200  gGeoManager->SetPdgName(pdgDatabase->GetParticle(ion->GetName())->PdgCode(), ion->GetName());
1201  }
1202  LOG(info) << "Add Ion: " << ion->GetName() << " with PDG "
1203  << pdgDatabase->GetParticle(ion->GetName())->PdgCode();
1204  }
1205  }
1206  delete Iter;
1208  if (fEvGen) {
1209  fEvGen->Init();
1210  }
1211 }
1212 //_____________________________________________________________________________
1214 {
1215  TObjArray* NewPart = fRun->GetUserDefParticles();
1216  TIterator* parIter = NewPart->MakeIterator();
1217  parIter->Reset();
1218 
1219  // check MC engine is not null (fMC is 0x00 at this line in case of Geant4)
1220  TVirtualMC* curMC = TVirtualMC::GetMC();
1221  if (!curMC)
1222  LOG(fatal) << "No MC engine was defined before AddParticles()";
1223 
1224  TObject* obj = 0;
1225  FairParticle* particle = 0;
1226  while ((obj = parIter->Next())) {
1227  particle = dynamic_cast<FairParticle*>(obj);
1228  if (particle) { // (Int_t pdg, const char* name, TMCParticleType type, Double_t mass, Double_t charge,
1229  // Double_t lifetime);
1230  LOG(info) << "Add Particle: " << particle->GetName() << " with PDG " << particle->GetPDG() << "\n"
1231  << particle->GetName() << " // const TString& name \n"
1232  << particle->GetMCType() << " // TMCParticleType mcType \n"
1233  << particle->GetMass() << " // Double_t mass \n"
1234  << particle->GetCharge() << " // Double_t charge \n"
1235  << particle->GetDecayTime() << " // Double_t lifetime \n"
1236  << particle->GetPType() << " // const TString& pType, \n"
1237  << particle->GetWidth() << " // Double_t width \n"
1238  << particle->GetSpin() << " // Int_t iSpin \n"
1239  << particle->GetiParity() << " // Int_t iParity \n"
1240  << particle->GetConjugation() << " // Int_t iConjugation \n"
1241  << particle->GetIsospin() << " // Int_t iIsospin \n"
1242  << particle->GetIsospinZ() << " // Int_t iIsospinZ \n"
1243  << particle->GetgParity() << " // Int_t gParity \n"
1244  << particle->GetLepton() << " // Int_t lepton \n"
1245  << particle->GetBaryon() << " // Int_t baryon \n"
1246  << particle->IsStable() << " // Bool_t stable \n";
1247 
1248  curMC->DefineParticle(particle->GetPDG(), // Int_t pdg
1249  particle->GetName(), // const TString& name
1250  particle->GetMCType(), // TMCParticleType mcType
1251  particle->GetMass(), // Double_t mass
1252  particle->GetCharge(), // Double_t charge
1253  particle->GetDecayTime(), // Double_t lifetime
1254  particle->GetPType(), // const TString& pType,
1255  particle->GetWidth(), // Double_t width
1256  particle->GetSpin(), // Int_t iSpin
1257  particle->GetiParity(), // Int_t iParity
1258  particle->GetConjugation(), // Int_t iConjugation
1259  particle->GetIsospin(), // Int_t iIsospin
1260  particle->GetIsospinZ(), // Int_t iIsospinZ
1261  particle->GetgParity(), // Int_t gParity
1262  particle->GetLepton(), // Int_t lepton
1263  particle->GetBaryon(), // Int_t baryon
1264  particle->IsStable() // Bool_t stable
1265  );
1266 
1267  // Add Particle to gGeoManager visualization
1268  if (gGeoManager)
1269  gGeoManager->SetPdgName(particle->GetPDG(), particle->GetName());
1270  }
1271  }
1272 
1273  delete parIter;
1274  AddDecayModes();
1275 }
1276 
1277 //_____________________________________________________________________________
1279 {
1280  TString work = getenv("VMCWORKDIR");
1281  TString work_config = work + "/gconfig/";
1282  work_config.ReplaceAll("//", "/");
1283 
1284  TString config_dir = getenv("CONFIG_DIR");
1285  config_dir.ReplaceAll("//", "/");
1286 
1287  Bool_t AbsPath = kFALSE;
1288 
1289  if (!config_dir.EndsWith("/")) {
1290  config_dir += "/";
1291  }
1292  // set Pythia as external decayer
1293 
1294  if (fPythiaDecayer) {
1295  TString decayConfig;
1296  if (fPythiaDecayerConfig.IsNull()) {
1297  decayConfig = "DecayConfig.C";
1298  fPythiaDecayerConfig = decayConfig;
1299  } else {
1300  if (fPythiaDecayerConfig.Contains("/")) {
1301  AbsPath = kTRUE;
1302  }
1303  decayConfig = fPythiaDecayerConfig;
1304  }
1305 
1306  if (!AbsPath && TString(gSystem->FindFile(config_dir.Data(), decayConfig)) != TString("")) {
1307  LOG(info) << "---User path for Configuration (DecayConfig.C) is used : " << config_dir.Data();
1308  } else {
1309  if (AbsPath) {
1310  decayConfig = fPythiaDecayerConfig;
1311  } else {
1312  decayConfig = work_config + fPythiaDecayerConfig;
1313  }
1314  }
1315  // Add decay modes using an external configuration script
1316  LOG(info) << "External Decay Modes with script \n " << decayConfig.Data();
1317  // Load configuration script and execute it
1318  Int_t pyt = gROOT->LoadMacro(decayConfig.Data());
1319  if (pyt == 0) {
1320  gInterpreter->ProcessLine("DecayConfig()");
1321  }
1322  }
1323  // set user defined phase space decay for particles (ions)
1324  AbsPath = kFALSE;
1325  if (fUserDecay) {
1326  TString Userdecay;
1327  if (fUserDecayConfig.IsNull()) {
1328  Userdecay = "UserDecay.C";
1329  fUserDecayConfig = Userdecay;
1330  } else {
1331  if (fUserDecayConfig.Contains("/")) {
1332  AbsPath = kTRUE;
1333  }
1334  Userdecay = fUserDecayConfig;
1335  }
1336 
1337  if (!AbsPath && TString(gSystem->FindFile(config_dir.Data(), Userdecay)) != TString("")) {
1338  LOG(info) << "---User path for Configuration (UserDecay.C) is used : " << config_dir.Data();
1339  } else {
1340  if (AbsPath) {
1341  Userdecay = fUserDecayConfig;
1342  } else {
1343  Userdecay = work_config + fUserDecayConfig;
1344  }
1345  }
1346  LOG(info) << "User Decay Modes with script \n " << Userdecay.Data();
1347  Int_t dec = gROOT->LoadMacro(Userdecay.Data());
1348  if (dec == 0) {
1349  gInterpreter->ProcessLine("UserDecayConfig()");
1350  }
1351  }
1352 }
1353 
1354 //_____________________________________________________________________________
1356 
1357 //_____________________________________________________________________________
1359 
1360 //_____________________________________________________________________________
1361 void FairMCApplication::AddTask(TTask* fTask)
1362 {
1363  if (!fFairTaskList) {
1364  fFairTaskList = new FairTask("Task List", 1);
1365  gROOT->GetListOfBrowsables()->Add(fFairTaskList);
1366  }
1367  fFairTaskList->Add(fTask);
1368  SetParTask();
1369 }
1370 
1371 //_____________________________________________________________________________
1373 
1374 //_____________________________________________________________________________
1376 
1377 //_____________________________________________________________________________
1379 {
1380  // Only RTDB init when more than Main Task list
1381  if (FairRun::Instance()->GetNTasks() >= 1) {
1383  }
1384  fModIter->Reset();
1385  FairModule* Mod = nullptr;
1386  while ((Mod = dynamic_cast<FairModule*>(fModIter->Next()))) {
1387  Mod->SetParContainers();
1388  }
1390  fRTdb->initContainers(FairRunSim::Instance()->GetRunId());
1391 }
1392 //_____________________________________________________________________________
1394 {
1395 
1396  // Only RTDB init when more than Main Task list
1397  if (FairRun::Instance()->GetNTasks() >= 1) {
1398  LOG(info) << "Initialize Tasks--------------------------";
1400  }
1401 }
1402 
1403 //_____________________________________________________________________________
1405 {
1406  if (fRootManager) {
1407  return fRootManager->GetInChain();
1408  } else {
1409  LOG(warn) << "The function is not available in MT mode";
1410  return 0;
1411  }
1412 }
1413 
1414 //_____________________________________________________________________________
1416 {
1417  fRadLength = RadLen;
1418  if (fRadLength) {
1419  fRadLenMan = new FairRadLenManager();
1420  }
1421 }
1422 
1423 //_____________________________________________________________________________
1425 {
1426  fRadMap = RadMap;
1427  if (fRadMap) {
1428  fRadMapMan = new FairRadMapManager();
1429  }
1430 }
1431 
1432 //_____________________________________________________________________________
1433 void FairMCApplication::AddMeshList(TObjArray* meshList)
1434 {
1435  if (!fRadGridMan) {
1437  }
1438  fRadGridMan->AddMeshList(meshList);
1439 }
1440 
1441 //_____________________________________________________________________________
1442 Int_t FairMCApplication::GetIonPdg(Int_t z, Int_t a) const
1443 {
1444  // Acording to
1445  // http://pdg.lbl.gov/2012/reviews/rpp2012-rev-monte-carlo-numbering.pdf
1446 
1447  return 1000000000 + 10 * 1000 * z + 10 * a;
1448 }
1449 
1450 //_____________________________________________________________________________
1451 void FairMCApplication::UndoGeometryModifications()
1452 {
1453  // Undo all misalignment done in the MisalignGeometry methods of the
1454  // several FairModuls.
1455  // In the output (parameter container and separate geometry file)
1456  // only the ideal geometry is stored.
1457  // I don't know any better way than to loop over all physical nodes
1458  // and to set the matrix back to the original one.
1459  // TODO: Check if it is more easy to write the ideal geometry before
1460  // the geometry is misaligned. In this case one does not have
1461  // to revert the misalignment.
1462 
1463  TObjArray* physNodes = gGeoManager->GetListOfPhysicalNodes();
1464  Int_t numPhysNodes = physNodes->GetEntriesFast();
1465 
1466  if (0 == numPhysNodes)
1467  return;
1468 
1469  // fRootManager->CreateGeometryFile("misaligned_geometry.root");
1470  LOG(info) << "Undo all misalignment";
1471 
1472  TGeoPhysicalNode* node = nullptr;
1473  TGeoHMatrix* ng3 = nullptr;
1474  for (Int_t k = 0; k < numPhysNodes; k++) {
1475  node = static_cast<TGeoPhysicalNode*>(physNodes->At(k));
1476  ng3 = node->GetOriginalMatrix(); //"real" global matrix, what survey sees
1477  node->Align(ng3);
1478  }
1479 
1480  gGeoManager->ClearPhysicalNodes(kFALSE);
1481 }
1482 
virtual void FinishRunOnWorker()
FairMCApplicationState fState
virtual FairGenericStack * CloneStack() const
virtual const char * GetName() const
Definition: FairParticle.h:74
void AlignGeometry() const
Definition: FairRun.cxx:181
virtual void ConstructGeometry()
Definition: FairModule.cxx:60
list of container factories
Definition: FairRuntimeDb.h:24
TList * getListOfMedia()
Definition: FairGeoMedia.h:37
FairGeoMedia * getMedia()
TRefArray * fActiveDetectors
virtual void AddParticles()
Int_t GetBaryon()
Definition: FairParticle.h:88
static thread_local TRefArray * svList
Definition: FairModule.h:137
Double_t GetMass() const
Definition: FairIon.h:78
void WriteInfo()
Definition: FairRunInfo.cxx:79
virtual void ExecuteTask(Option_t *option="0")
Definition: FairTask.cxx:118
virtual void ConstructOpGeometry()
Definition: FairModule.cxx:66
Bool_t IsStable()
Definition: FairParticle.h:89
TChain * GetInChain()
Int_t GetA() const
Definition: FairIon.h:66
virtual void FinishEvent()
virtual void RegisterAlignmentMatrices()
Definition: FairModule.h:77
virtual void Register()=0
const TString & GetPType()
Definition: FairParticle.h:79
virtual void SetParContainers()
Definition: FairModule.h:125
virtual void StopRun()
TH2D * GetMeshTid()
Definition: FairMesh.h:65
virtual void SetSpecialPhysicsCuts()
Definition: FairModule.h:83
FairEventHeader * fEventHeader
Int_t GetLepton()
Definition: FairParticle.h:87
std::map< Int_t, Int_t > fModVolMap
FairRootManager * fRootManager
virtual Double_t TrackingRmax() const
Int_t GetQ() const
Definition: FairIon.h:70
Int_t getCopyNo()
Definition: FairVolume.h:58
Int_t GetZ() const
Definition: FairIon.h:62
FairMCEventHeader * GetMCEventHeader()
Definition: FairRunSim.cxx:357
FairGeoInterface * getGeoInterface()
Definition: FairGeoLoader.h:34
virtual void Stepping()
virtual Double_t TrackingZmax() const
void SetDetArrayList(TRefArray *detArray)
static FairGeoLoader * Instance()
static FairRun * Instance()
Definition: FairRun.cxx:31
TObjArray * GetMeshList()
FairGeoNode * getGeoNode()
Definition: FairVolume.h:60
virtual void InitGeometry()
virtual void UpdateTrackIndex(TRefArray *)
void SetRadiationMapReg(Bool_t RadMap)
Double_t GetStepSizeCut() const
void getCerenkovPar(Int_t, Double_t *)
virtual Bool_t MisalignGeometry()
Bool_t IsImportTGeoToVMC() const
Definition: FairRunSim.h:188
std::map< Int_t, Int_t >::iterator fModVolIter
static FairRootManager * Instance()
Int_t GetgParity()
Definition: FairParticle.h:86
static FairRunSim * Instance()
Definition: FairRunSim.cxx:116
FairSink * GetSink()
Definition: FairRun.h:93
Bool_t IsAccepted(const TParticle *p) const
ClassImp(FairEventBuilder)
FairDetector * fDisDet
FairGenericStack * GetStack()
FairRadGridManager * fRadGridMan
virtual void FinishPrimary()
Int_t getModId()
Definition: FairVolume.h:39
std::list< FairDetector * > listActiveDetectors
FairMCApplicationState
void SetSink(FairSink *tempSink)
Definition: FairRun.h:84
FairDetector * GetDetector(const char *DetName)
virtual void ConstructGeometry()
virtual void Register()
FairRadLenManager * fRadLenMan
std::mutex mtx
Double_t GetCharge()
Definition: FairParticle.h:77
Int_t GetPDG() const
std::multimap< Int_t, FairVolume * > fVolMap
void SetRadiationLengthReg(Bool_t RadLen)
void SetModule(FairModule *mod)
Definition: FairVolume.h:49
Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
void AddMeshList(TObjArray *meshList)
virtual void Initialize()
Bool_t initContainers(Int_t runId, Int_t refId=-1, const Text_t *fileName="")
FairModule * GetModule()
Definition: FairVolume.h:47
virtual TVirtualMCApplication * CloneForWorker() const
FairPrimaryGenerator * GetPrimaryGenerator()
Definition: FairRunSim.h:121
TLorentzVector fTrkPos
void Scale(Double_t fac)
Definition: FairMesh.h:77
TGeoTrack * GetCurrentTrk()
TH2D * GetMeshFlu()
Definition: FairMesh.h:66
Int_t getNpckov()
Definition: FairGeoMedium.h:62
void StoreInfo()
Definition: FairRunInfo.cxx:36
virtual void PreTrack()
TObjArray * GetUserDefIons()
Definition: FairRunSim.cxx:131
void SetField(FairField *field)
virtual void FillTrackArray()
Int_t GetSpin()
Definition: FairParticle.h:81
virtual FairPrimaryGenerator * ClonePrimaryGenerator() const
TGeoTrack * CheckAddTrack(Int_t trackId, TParticle *p)
virtual void Reset()
virtual void ConstructOpGeometry()
std::list< FairDetector * > listDetectors
virtual void FinishTask()
Definition: FairTask.cxx:90
virtual void GeneratePrimaries()
virtual void Register()
virtual void AddDecayModes()
FairMCEventHeader * fMCEventHeader
Double_t GetWidth()
Definition: FairParticle.h:80
FairPrimaryGenerator * GetGenerator()
Double_t GetMass()
Definition: FairParticle.h:76
static FairTrajFilter * Instance()
FairDetector * GetDetector()
Definition: FairVolume.h:48
Int_t getMCid()
Definition: FairVolume.h:57
Int_t GetInstanceId() const
void InitTask()
Definition: FairTask.cxx:41
virtual Sink_Type GetSinkType()=0
virtual void InitOnWorker()
virtual void InitParContainers()
Definition: FairModule.h:128
void SetRunID(UInt_t runId)
virtual Bool_t GenerateEvent(FairGenericStack *pStack)
FairRuntimeDb * GetRuntimeDb(void)
Definition: FairRun.h:80
void SetEvent(FairMCEventHeader *event)
void AddPoint(Int_t &ModuleId)
Bool_t GetStoreTraj() const
Definition: FairRunSim.h:95
virtual FairSink * CloneSink()=0
void SetEventID(UInt_t eventId)
FairPrimaryGenerator * fEvGen
TH2D * GetMeshSEU()
Definition: FairMesh.h:67
FairGenericStack * fStack
TMCParticleType GetMCType()
Definition: FairParticle.h:75
virtual void BeginPrimary()
Int_t GetConjugation()
Definition: FairParticle.h:83
virtual void StopMCRun()
virtual Bool_t ProcessHits(FairVolume *v=0)=0
virtual void FinishPrimary()
virtual void BeginEvent()
void AddTask(TTask *fTask)
void setCopyNo(Int_t id)
Definition: FairVolume.h:41
std::multimap< Int_t, FairVolume * >::iterator fVolIter
virtual void FinishEvent()
Definition: FairTask.cxx:82
void setModId(Int_t id)
Definition: FairVolume.h:40
virtual void PostTrack()
Int_t GetIsospin()
Definition: FairParticle.h:84
void setMCid(Int_t id)
Definition: FairVolume.h:59
Double_t GetDecayTime()
Definition: FairParticle.h:78
void AddPoint(Int_t &ModuleId)
Int_t GetiParity()
Definition: FairParticle.h:82
Int_t getMediumIndex()
Definition: FairGeoMedium.h:49
Bool_t IsActive()
Definition: FairModule.h:115
void InitMC(const char *setup, const char *cuts)
void SetGenerator(FairPrimaryGenerator *fxGenerator)
FairTrajFilter * fTrajFilter
Int_t GetRunId()
Definition: FairRun.h:97
UInt_t GetEventID() const
run identifier
Int_t GetModId()
Definition: FairModule.h:111
void Init(TString brName="GeoTracks", TString folderName="MCGeoTrack")
void RunMC(Int_t nofEvents)
void AddMeshList(TObjArray *list)
void SetParTask()
Definition: FairTask.cxx:73
FairRadMapManager * fRadMapMan
Int_t GetIsospinZ()
Definition: FairParticle.h:85
TObjArray * GetUserDefParticles()
Definition: FairRunSim.cxx:137
Double_t GetExcEnergy() const
Definition: FairIon.h:74