8 #include "TDatabasePDG.h"
9 #include "TG4RunConfiguration.h"
11 #include "TGeant3TGeo.h"
13 #include "TGeoManager.h"
14 #include "TPythia6Decayer.h"
16 #include "TStopwatch.h"
20 #include "TVirtualMC.h"
22 #include <boost/filesystem.hpp>
28 #include "FairLogger.h"
29 #include "FairMonitor.h"
30 #include "FairParRootFileIo.h"
31 #include "FairRunSim.h"
32 #include "FairRuntimeDb.h"
33 #include "FairSystemInfo.h"
34 #include "FairUrqmdGenerator.h"
49 using std::stringstream;
54 : TNamed(
"CbmTransport",
"Transport Run")
60 , fRun(new FairRunSim())
70 , fGenerateRunInfo(kFALSE)
71 , fStoreTrajectories(kFALSE) {
87 auto pdgdb = TDatabasePDG::Instance();
88 pdgdb->ReadPDGTable();
105 FairGenerator* generator = NULL;
107 if (gSystem->AccessPathName(fileName)) {
108 LOG(fatal) << GetName() <<
": Input file " << fileName <<
" not found!";
114 case kUrqmd: generator =
new FairUrqmdGenerator(fileName);
break;
135 std::cout << std::endl;
136 LOG(info) << GetName() <<
": Configuring VMC...";
137 TVirtualMC* vmc =
nullptr;
140 TString* gModel =
fRun->GetGeoModel();
141 if (strncmp(gModel->Data(),
"TGeo", 4) == 0) {
142 LOG(info) << GetName() <<
": Create TGeant3TGeo";
143 vmc =
new TGeant3TGeo(
"C++ Interface to Geant3 with TGeo");
146 LOG(info) << GetName() <<
": Create TGeant3";
147 vmc =
new TGeant3(
"C++ Interface to Geant3");
154 LOG(info) << GetName() <<
": Create TGeant4";
156 std::array<std::string, 3> runConfigSettings =
158 TG4RunConfiguration* runConfig =
new TG4RunConfiguration(
159 runConfigSettings[0], runConfigSettings[1], runConfigSettings[2]);
160 vmc =
new TGeant4(
"TGeant4",
"C++ Interface to Geant4", runConfig);
165 LOG(fatal) << GetName() <<
": unknown transport engine!";
168 std::unique_ptr<CbmStack> stack(
new CbmStack());
170 if (vmc) vmc->SetStack(stack.release());
195 auto pdgdb = TDatabasePDG::Instance();
199 LOG(fatal) << GetName()
200 <<
": Forcing decay modes is not possible with TGeant4!";
204 Int_t pdg = decay.first;
205 UInt_t nDaughters = decay.second.size();
207 log << GetName() <<
": Force decay " << pdgdb->GetParticle(pdg)->GetName()
219 if (gMC->ParticleMCType(pdg)
221 LOG(info) <<
log.str();
222 LOG(fatal) << GetName() <<
": PDG " << pdg <<
" not in VMC!";
227 if (nDaughters <= 3) {
228 Float_t branch[6] = {100., 0., 0., 0., 0., 0.};
229 Int_t mode[6][3] = {{0, 0, 0},
235 for (UInt_t iDaughter = 0; iDaughter < nDaughters; iDaughter++) {
236 mode[0][iDaughter] = decay.second[iDaughter];
237 log << pdgdb->GetParticle(decay.second[iDaughter])->GetName() <<
" ";
239 Bool_t success = gMC->SetDecayMode(pdg, branch, mode);
241 LOG(info) <<
log.str();
242 LOG(fatal) << GetName() <<
": Setting decay mode failed!";
244 log <<
", using native decayer.";
249 auto p6decayer = TPythia6Decayer::Instance();
250 Int_t daughterPdg[nDaughters];
251 Int_t multiplicity[nDaughters];
252 for (UInt_t iDaughter = 0; iDaughter < nDaughters; iDaughter++) {
253 daughterPdg[iDaughter] = decay.second[iDaughter];
254 multiplicity[iDaughter] = 1;
255 log << pdgdb->GetParticle(decay.second[iDaughter])->GetName() <<
" ";
257 p6decayer->ForceParticleDecay(pdg, daughterPdg, multiplicity, nDaughters);
259 gMC->SetUserDecay(pdg);
262 LOG(info) <<
log.str();
274 std::cout << std::endl;
275 LOG(info) <<
"----- Settings for event generator";
277 LOG(info) <<
"----- End settings for event generator";
278 std::cout << std::endl;
284 LOG(fatal) << GetName()
285 <<
": Beam profile is not consistent with target!";
319 TDatabasePDG* pdgdb = TDatabasePDG::Instance();
320 const char* name =
"";
323 Bool_t stable = kTRUE;
324 Double_t charge = 0.;
332 pdgdb->AddParticle(name, name, mass, stable, 0., charge,
"Ion", code);
333 pdgdb->AddAntiParticle(
"d-", -1 * code);
341 pdgdb->AddParticle(name, name, mass, stable, 0., charge,
"Ion", code);
342 pdgdb->AddAntiParticle(
"t-", -1 * code);
350 pdgdb->AddParticle(name, name, mass, stable, 0., charge,
"Ion", code);
351 pdgdb->AddAntiParticle(
"He3-", -1 * code);
359 pdgdb->AddParticle(name, name, mass, stable, 0., charge,
"Ion", code);
360 pdgdb->AddAntiParticle(
"He3-", -1 * code);
368 fRun->SetRadLenRegister(choice);
369 LOG(info) << GetName() <<
": Radiation length register is enabled";
383 LOG(info) << GetName() <<
": Set decay modes for pi0 and eta";
386 std::cout << std::endl << std::endl;
387 LOG(warn) <<
"***********************************";
388 LOG(warn) <<
"***********************************";
389 LOG(warn) << GetName() <<
": User decay modes cannot be set with TGeant4!";
390 LOG(warn) << GetName()
391 <<
": Built-in decay modes for pi0 and eta will be used!";
392 LOG(warn) <<
"***********************************";
393 LOG(warn) <<
"***********************************";
394 std::cout << std::endl << std::endl;
398 TGeant3* gMC3 =
static_cast<TGeant3*
>(vmc);
430 gMC3->Gsdk(17, brEta, modeEta);
445 for (Int_t iMode = 2; iMode < 6; iMode++) {
451 gMC3->Gsdk(7, brPi, modePi);
464 LOG(fatal) << GetName() <<
": No output file specified!";
466 LOG(fatal) << GetName() <<
": No parameter file specified!";
467 std::cout << std::endl << std::endl;
475 const char* engineName =
"";
477 case kGeant3: engineName =
"TGeant3";
break;
478 case kGeant4: engineName =
"TGeant4";
break;
480 LOG(warn) << GetName() <<
": Unknown transport engine ";
481 engineName =
"TGeant3";
485 LOG(info) << GetName() <<
": Using engine " << engineName;
486 fRun->SetName(engineName);
497 LOG(info) <<
fTarget->ToString();
500 LOG(warn) << GetName() <<
": No target defined!";
504 LOG(info) << GetName() <<
": Register magnetic field";
525 fRun->SetSimSetup(
f);
547 FairRuntimeDb* rtdb =
fRun->GetRuntimeDb();
549 static_cast<CbmFieldPar*
>(rtdb->getContainer(
"CbmFieldPar"));
551 fieldPar->setChanged();
552 FairParRootFileIo* parOut =
new FairParRootFileIo(kTRUE);
554 rtdb->setOutput(parOut);
564 TDirectory* oldDir = gDirectory;
567 TFile* outfile =
fRun->GetOutputFile();
570 LOG(info) <<
"Here I am";
572 LOG(info) <<
"Write Geant3Settings";
575 LOG(info) <<
"Write Geant4Settings";
593 std::cout << std::endl;
594 LOG(info) << GetName() <<
": Run finished successfully.";
595 LOG(info) << GetName() <<
": Wall time for Init : " <<
fRealTimeInit <<
" s ";
596 LOG(info) << GetName() <<
": Wall time for Run : " <<
fRealTimeRun <<
" s ("
599 LOG(info) << GetName() <<
": Output file : " <<
fOutFileName;
600 LOG(info) << GetName() <<
": Parameter file : " <<
fParFileName;
602 LOG(info) << GetName() <<
": Geometry file : " <<
fGeoFileName;
603 std::cout << std::endl << std::endl;
609 if (gROOT->GetVersionInt() >= 60602) {
610 gGeoManager->GetListOfVolumes()->Delete();
611 gGeoManager->GetListOfShapes()->Delete();
645 Int_t* daughterPdg) {
648 LOG(fatal) << GetName() <<
": User decay mode for PDG " << pdg
649 <<
" is already defined!";
653 for (UInt_t iDaughter = 0; iDaughter < nDaughters; iDaughter++) {
654 fDecayModes[pdg].push_back(daughterPdg[iDaughter]);
664 std::string name = fileName.Data();
665 Int_t found = name.find_last_of(
"/");
667 TString geoDir = name.substr(0, found);
668 if (gSystem->AccessPathName(geoDir.Data())) {
669 LOG(error) << GetName() <<
": Directory for geometry file " << geoDir
670 <<
" does not exist; the file will not be created.";
684 if (gSystem->AccessPathName(fileName)) {
685 std::string name = fileName.Data();
686 Int_t found = name.find_last_of(
"/");
688 TString parDir = name.substr(0, found);
689 if (gSystem->AccessPathName(parDir.Data())) {
690 LOG(fatal) << GetName() <<
": Parameter directory " << parDir
691 <<
" does not exist!";
704 fEventGen->SetEventPlane(phiMin, phiMax);
713 std::string name = fileName.Data();
714 Int_t found = name.find_last_of(
"/");
716 TString outDir = name.substr(0, found);
717 if (gSystem->AccessPathName(outDir.Data())) {
718 LOG(fatal) << GetName() <<
": Output directory " << outDir
719 <<
" does not exist!";