added runMC_v2.C cleared from ifdef directives

generators/CMakeLists.txt

@@ -8,6 +8,8 @@ STRING(REGEX REPLACE "^[0-9]+\\.([0-9]+)\\.[0-9]" "\\1" req_gcc_minor_vers "${GC
 IF(${req_gcc_major_vers} MATCHES "4" AND NOT ${req_gcc_minor_vers} MATCHES "0")
 ENDIF(${req_gcc_major_vers} MATCHES "4" AND NOT ${req_gcc_minor_vers} MATCHES "0")
 
+add_subdirectory(genFactory)
+
 set(INCLUDE_DIRECTORIES
  ${BASE_INCLUDE_DIRECTORIES}
  ${CMAKE_SOURCE_DIR}/config

generators/genFactory/CMakeLists.txt

+set(INCLUDE_DIRECTORIES
+${BASE_INCLUDE_DIRECTORIES}
+ ${CMAKE_SOURCE_DIR}/generators
+ ${CMAKE_SOURCE_DIR}/generators/genFactory
+ ${CMAKE_SOURCE_DIR}/shield_pack/THadgen/inc
+ ${CMAKE_SOURCE_DIR}/shield_pack/hadgen/inc
+ ${CMAKE_SOURCE_DIR}/macro/mpd
+ ${CMAKE_SOURCE_DIR}
+)
+
+Set(SYSTEM_INCLUDE_DIRECTORIES
+ ${ROOT_INCLUDE_DIR}
+ ${FAIRROOT_LIBRARY_DIR}
+)
+
+include_directories(${INCLUDE_DIRECTORIES})
+Include_Directories(SYSTEM ${SYSTEM_INCLUDE_DIRECTORIES})
+
+set(LINK_DIRECTORIES
+ ${ROOT_LIBRARY_DIR}
+ ${FAIRROOT_LIBRARY_DIR}
+)
+ 
+link_directories(${LINK_DIRECTORIES})
+
+# List of source files
+set(SRCS
+MpdGeneratorsFactory.cxx
+)
+
+Set(HEADERS)
+Set(LINKDEF LinkDefFactory.h)
+Set(LIBRARY_NAME MpdGenFactory)
+Set(DEPENDENCIES Core Base MpdGen)
+
+GENERATE_LIBRARY()

generators/genFactory/LinkDefFactory.h

+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ defined_in MpdGeneratorType.h;
+#pragma link C++ defined_in MpdGeneratorsFactory.h;
+
+#endif

generators/genFactory/MpdGeneratorType.h

+#ifndef MPDGENERATORTYPE_H
+#define MPDGENERATORTYPE_H
+#include <TString.h>
+
+enum class MpdGeneratorType
+{
+	URQMD,
+	VHLLE,
+	FLUID,
+	PART,
+	ION,
+	BOX,
+	HSD,
+	LAQGSM,
+	HADGEN,
+	CUSTOM = 1000
+};
+
+struct MpdGenerator
+{
+	MpdGeneratorType genType;
+	TString inFile;
+	Int_t startEvent;
+	Int_t nEvents;
+};
+#endif

generators/genFactory/MpdGeneratorsFactory.cxx

+#include "MpdGeneratorsFactory.h"
+
+#include <TMath.h>
+#include <TRandom.h>
+#include <TDatabasePDG.h>
+
+#include "MpdUrqmdGenerator.h"
+#include "MpdVHLLEGenerator.h"
+#include "Mpd3fdGenerator.h"
+#include "FairParticleGenerator.h"
+#include "FairIonGenerator.h"
+#include "FairBoxGenerator.h"
+#include "MpdPHSDGenerator.h"
+#include "MpdLAQGSMGenerator.h"
+#include "THadgen.h"
+#include "MpdGetNumEvents.h"
+#include "mpdloadlibs.C"
+
+
+MpdGeneratorsFactory::MpdGeneratorsFactory()
+{
+	add<MpdUrqmdGenCreator>(MpdGeneratorType::URQMD);
+	add<MpdVHLLEGenCreator>(MpdGeneratorType::VHLLE);
+	add<MpdFLUIDGenCreator>(MpdGeneratorType::FLUID);
+	//add<MpdPARTGenCreator>(MpdGeneratorType::PART);
+	//add<MpdIONGenCreator>(MpdGeneratorType::ION);
+	//add<MpdBOXGenCreator>(MpdGeneratorType::BOX);
+	add<MpdHSDGenCreator>(MpdGeneratorType::HSD);
+	add<MpdLAQGSMGenCreator>(MpdGeneratorType::LAQGSM);
+	//add<MpdHADGENGenCreator>(MpdGeneratorType::HADGEN);
+}
+	
+MpdGeneratorsFactory::~MpdGeneratorsFactory()
+{
+	for(CreatorsMap::iterator it = creators.begin(); it != creators.end(); ++it)
+		delete it->second;
+}
+
+std::shared_ptr<MpdFactoryMadeGenerator> MpdGeneratorsFactory::create(MpdGenerator & Gen)
+{
+	typename CreatorsMap::iterator it = creators.find(Gen.genType);
+	if (it != creators.end())
+	{
+		std::shared_ptr<FairGenerator> gen = std::shared_ptr<FairGenerator>(it->second->create(Gen.inFile, Gen.startEvent, Gen.nEvents));
+		return std::make_shared<MpdFactoryMadeGenerator>(gen, it->second->postActions(Gen.inFile), Gen.inFile);
+	}
+	return std::shared_ptr<MpdFactoryMadeGenerator>(NULL);
+}
+
+template <class C> 
+void MpdGeneratorsFactory::add(const MpdGeneratorType & GenType)
+{
+	typename CreatorsMap::iterator it = creators.find(GenType);
+	if (it == creators.end())
+		creators[GenType] = new C();
+}
+
+FairGenerator * MpdUrqmdGenCreator::create (TString & inFile, Int_t & nStartEvent, Int_t & nEvents)
+{
+		// ------- Urqmd  Generator
+	if (!CheckFileExist(inFile)) return NULL;
+
+	MpdUrqmdGenerator* urqmdGen = new MpdUrqmdGenerator(inFile);
+	// Event plane angle will be generated by uniform distribution from min to max.
+	// Angles are in degrees
+	Float_t min = 0.0;
+	Float_t max = 30.0;
+	urqmdGen->SetEventPlane(min * TMath::DegToRad(), max * TMath::DegToRad());
+	if (nStartEvent > 0) urqmdGen->SkipEvents(nStartEvent);
+
+	// if nEvents is equal 0 then all events (start with nStartEvent) of the given file should be processed
+	if (nEvents == 0)
+		nEvents = MpdGetNumEvents::GetNumURQMDEvents(const_cast<char *>(inFile.Data())) - nStartEvent;
+
+	return urqmdGen;
+}
+
+FairGenerator * MpdVHLLEGenCreator::create (TString & inFile, Int_t &, Int_t &)
+{
+	if (!CheckFileExist(inFile)) return NULL;
+	
+	MpdVHLLEGenerator* vhlleGen = new MpdVHLLEGenerator(inFile, kTRUE); // kTRUE corresponds to hydro + cascade, kFALSE -- hydro only
+	vhlleGen->SkipEvents(0);
+	return vhlleGen;
+}
+
+FairGenerator * MpdFLUIDGenCreator::create (TString & inFile, Int_t & nStartEvent, Int_t &)
+{
+	if (!CheckFileExist(inFile)) return NULL;
+	
+	Mpd3fdGenerator* fluidGen = new Mpd3fdGenerator(inFile);
+    if (nStartEvent > 0) fluidGen->SkipEvents(nStartEvent);
+    //fluidGen->SetPsiRP(0.); // set fixed Reaction Plane angle [rad] instead of random
+    //fluidGen->SetProtonNumberCorrection(79./197.); // Z/A Au for Theseus 2018-03-17-bc2a06d
+	return fluidGen;
+}
+
+/*FairGenerator * MpdPARTGenCreator::create (TString &, Int_t &, Int_t &)
+{
+	// ------- Particle Generator
+	FairParticleGenerator* partGen = new FairParticleGenerator(211, 10, 1, 0, 3, 1, 0, 0);
+	return partGen;
+}*/
+
+/*FairGenerator * MpdIONGenCreator::create (TString &, Int_t &, Int_t &)
+{
+	// ------- Ion Generator
+	FairIonGenerator *fIongen = new FairIonGenerator(79, 197, 79, 1, 0., 0., 25, 0., 0., -1.);
+	return fIongen;
+}*/
+
+/*FairGenerator * MpdBOXGenCreator::create (TString &, Int_t &, Int_t &)
+{
+	gRandom->SetSeed(0);
+	// ------- Box Generator
+	FairBoxGenerator* boxGen = new FairBoxGenerator(13, 100); // 13 = muon; 1 = multipl.
+	boxGen->SetPRange(0.25, 2.5); // GeV/c //setPRange vs setPtRange
+	boxGen->SetPhiRange(0, 360); // Azimuth angle range [degree]
+	boxGen->SetThetaRange(0, 180); // Polar angle in lab system range [degree]
+	boxGen->SetXYZ(0., 0., 0.); // mm o cm ??
+	return boxGen;
+}*/
+
+FairGenerator * MpdHSDGenCreator::create (TString & inFile, Int_t & nStartEvent, Int_t & nEvents)
+{
+// ------- HSD/PHSD Generator
+	if (!CheckFileExist(inFile)) return NULL;
+
+	MpdPHSDGenerator *hsdGen = new MpdPHSDGenerator(inFile.Data());
+	//hsdGen->SetPsiRP(0.); // set fixed Reaction Plane angle [rad] instead of random
+	if (nStartEvent > 0) hsdGen->SkipEvents(nStartEvent);
+
+	// if nEvents is equal 0 then all events (start with nStartEvent) of the given file should be processed
+	if (nEvents == 0)
+		nEvents = MpdGetNumEvents::GetNumPHSDEvents(const_cast<char *>(inFile.Data())) - nStartEvent;
+
+	return hsdGen;
+}
+
+FairGenerator * MpdLAQGSMGenCreator::create (TString & inFile, Int_t & nStartEvent, Int_t & nEvents)
+{
+// ------- LAQGSM Generator
+	if (!CheckFileExist(inFile)) return NULL;
+
+	MpdLAQGSMGenerator* guGen = new MpdLAQGSMGenerator(inFile.Data(),kTRUE,0, 1+nStartEvent+nEvents);
+	// kTRUE - for NICA/MPD, 1+nStartEvent+nEvents - search ions in selected part of file.
+	if (nStartEvent > 0) guGen->SkipEvents(nStartEvent);
+
+	// if nEvents is equal 0 then all events (start with nStartEvent) of the given file should be processed
+	if (nEvents == 0)
+		nEvents = MpdGetNumEvents::GetNumQGSMEvents(const_cast<char *>(inFile.Data())) - nStartEvent;
+
+	return guGen;
+}
+
+std::function<Int_t(FairRunSim * fRun)> MpdLAQGSMGenCreator::postActions(TString & inFile)
+{
+    return [&inFile](FairRunSim * fRun)
+    {
+        TString Pdg_table_name = TString::Format("%s%s%c%s", gSystem->BaseName(inFile.Data()), ".g", (fRun->GetName())[6], ".pdg_table.dat");
+        return (TDatabasePDG::Instance())->WritePDGTable(Pdg_table_name.Data());
+    };
+}
+
+/*FairGenerator * MpdHADGENGenCreator::create (TString &, Int_t &, Int_t &)
+{
+	THadgen* hadGen = new THadgen();
+	hadGen->SetRandomSeed(clock() + time(0));
+	hadGen->SetParticleFromPdgCode(0, 196.9665, 79);
+	hadGen->SetEnergy(6.5E3);
+	MpdGeneralGenerator* generalHad = new MpdGeneralGenerator(hadGen);
+	return generalHad;
+}*/

generators/genFactory/MpdGeneratorsFactory.h

+#ifndef MPDGENERATORFACTORY_H
+#define MPDGENERATORFACTORY_H
+
+#include <map>
+#include <memory>
+#include <functional>
+
+#include <TString.h>
+#include <TSystem.h>
+
+#include "MpdGeneratorType.h"
+
+#include "FairGenerator.h"
+#include "FairRunSim.h"
+
+class MpdGenericGenCreator;
+class MpdUrqmdGenCreator;
+class MpdVHLLEGenCreator;
+class MpdFLUIDGenCreator;
+//class MpdPARTGenCreator;
+//class MpdIONGenCreator;
+//class MpdBOXGenCreator;
+class MpdHSDGenCreator;
+class MpdLAQGSMGenCreator;
+//class MpdHADGENGenCreator;
+class MpdFactoryMadeGenerator;
+
+
+class MpdGeneratorsFactory
+{
+public:
+	MpdGeneratorsFactory();
+	virtual ~MpdGeneratorsFactory();
+	
+	std::shared_ptr<MpdFactoryMadeGenerator> create(MpdGenerator & Gen);
+	
+protected:
+	typedef std::map<MpdGeneratorType, MpdGenericGenCreator*> CreatorsMap;
+	CreatorsMap creators;
+	
+	template <class C> 
+	void add(const MpdGeneratorType & id);
+};
+
+class MpdFactoryMadeGenerator: public std::enable_shared_from_this<MpdFactoryMadeGenerator>
+{
+public:
+	MpdFactoryMadeGenerator(std::shared_ptr<FairGenerator> fg, std::function<Int_t(FairRunSim*)> pa, TString dataFile)
+		: f_gen(fg), f_pa(pa), f_dataFile(dataFile) {}
+	std::shared_ptr<MpdFactoryMadeGenerator> getptr() {return shared_from_this();}
+	
+	virtual const std::shared_ptr<FairGenerator> & GetGeneratorPtr() const {return f_gen;}
+	virtual Int_t PostActions(FairRunSim * fRun){return f_pa(fRun);}
+	//virtual void SetPostActions(std::function<Int_t(FairRunSim*)> pa){f_pa = pa;}
+	
+protected:
+	std::shared_ptr<FairGenerator> f_gen;
+	TString f_dataFile = "";
+	std::function<Int_t(FairRunSim*)> f_pa;
+};
+
+class MpdGenericGenCreator
+{
+public:
+	virtual FairGenerator * create (TString & inFile, Int_t & nStartEvent, Int_t & nEvents) = 0;
+	virtual std::function<Int_t(FairRunSim*)> postActions(TString &){return [](FairRunSim*){return 0;};};
+};
+
+class MpdUrqmdGenCreator : public MpdGenericGenCreator
+{
+public:
+	virtual FairGenerator * create (TString & inFile, Int_t & nStartEvent, Int_t & nEvents) override;
+};
+
+class MpdVHLLEGenCreator : public MpdGenericGenCreator
+{
+public:
+	virtual FairGenerator * create (TString & inFile, Int_t &, Int_t &) override;
+};
+
+class MpdFLUIDGenCreator : public MpdGenericGenCreator
+{
+public:
+	virtual FairGenerator * create (TString & inFile, Int_t & nStartEvent, Int_t &) override;
+};
+
+/*class MpdPARTGenCreator : public MpdGenericGenCreator
+{
+public:
+	virtual FairGenerator * create (TString &, Int_t &, Int_t &) override;
+};*/
+
+/*class MpdIONGenCreator : public MpdGenericGenCreator
+{
+public:
+	virtual FairGenerator * create (TString &, Int_t &, Int_t &) override;
+};*/
+
+/*class MpdBOXGenCreator : public MpdGenericGenCreator
+{
+public:
+	virtual FairGenerator * create (TString &, Int_t &, Int_t &) override;
+};*/
+
+class MpdHSDGenCreator : public MpdGenericGenCreator
+{
+public:
+	virtual FairGenerator * create (TString & inFile, Int_t & nStartEvent, Int_t & nEvents) override;
+};
+
+class MpdLAQGSMGenCreator : public MpdGenericGenCreator
+{
+public:
+	virtual FairGenerator * create (TString & inFile, Int_t & nStartEvent, Int_t & nEvents) override;
+	virtual std::function<Int_t(FairRunSim * fRun)> postActions(TString & inFile) override;
+};
+
+/*class MpdHADGENGenCreator : public MpdGenericGenCreator
+{
+public:
+	virtual FairGenerator * create (TString &, Int_t &, Int_t &) override;
+};*/
+
+
+#endif

macro/mpd/runMC_v2.C

+#if !defined(__CINT__) && !defined(__CLING__)
+#include "TString.h"
+#include "TStopwatch.h"
+#include "TROOT.h"
+#include "TSystem.h"
+#include "TMath.h"
+#include "TDatabasePDG.h"
+
+#include "FairRun.h"
+#include "FairRunAna.h"
+#include "FairRunSim.h"
+#include "FairRuntimeDb.h"
+#include "FairParRootFileIo.h"
+#include "FairTrajFilter.h"
+#include "FairUrqmdGenerator.h"
+#include "FairPrimaryGenerator.h"
+#include "MpdLAQGSMGenerator.h"
+#include "FairCave.h"
+#include "FairPipe.h"
+#include "FairMagnet.h"
+
+#include "TpcDetector.h"
+#include "MpdEmc.h"
+#include "MpdTof.h"
+#include "MpdZdc.h"
+#include "MpdFfd.h"
+#include "MpdMultiField.h"
+#include "MpdMultiFieldPar.h"
+#include "MpdConstField.h"
+#include "MpdFieldMap.h"
+#include "MpdGetNumEvents.h"
+
+#include <iostream>
+#include <vector>
+#endif
+
+R__ADD_INCLUDE_PATH($VMCWORKDIR)
+#include "macro/mpd/mpdloadlibs.C"
+#include "macro/mpd/geometry_stage1.C"
+//#include "macro/mpd/geometry_v2.C"
+R__LOAD_LIBRARY(libMpdGenFactory);
+#include "generators/genFactory/MpdGeneratorsFactory.h"
+
+// inFile - input file with generator data, default: auau.09gev.mbias.98k.ftn14
+// nStartEvent - for compatibility, any number
+// nEvents - number of events to transport, default: 1
+// outFile - output file with MC data, default: evetest.root
+// flag_store_FairRadLenPoint
+// FieldSwitcher: 0 - corresponds to the ConstantField (0, 0, 5) kG (It is used by default); 1 - corresponds to the FieldMap ($VMCWORKDIR/input/B-field_v2.dat)
+void runMC_v2(TString inFile = "auau.04gev.0_3fm.10k.f14.gz", TString outFile = "evetest.root", Int_t nStartEvent = 0, Int_t nEvents = 10,
+        Bool_t flag_store_FairRadLenPoint = kFALSE, Int_t FieldSwitcher = 0,
+        MpdGeneratorType GenType = MpdGeneratorType::HADGEN /*Choose generator: URQMD VHLLE FLUID PART ION BOX HSD LAQGSM HADGEN CUSTOM*/, Bool_t useGeant4 = false)
+{
+    TStopwatch timer;
+    timer.Start();
+    gDebug = 0;
+
+    FairRunSim* fRun = new FairRunSim();
+    // Choose the Geant Navigation System
+    fRun->SetName(useGeant4 ? "TGeant4" : "TGeant3");
+
+    geometry_stage1(fRun); // load mpd geometry
+    //geometry_v2(fRun); // load mpd geometry
+
+    // Use extended MC Event header instead of the default one.
+    //MpdMCEventHeader* mcHeader = new MpdMCEventHeader();
+    //fRun->SetMCEventHeader(mcHeader);
+
+    // Create and Set Event Generator
+    FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
+    fRun->SetGenerator(primGen);
+
+    // smearing of beam interaction point
+    primGen->SetBeam(0.0,0.0,0.1,0.1);
+    primGen->SetTarget(0.0,24.0);
+    primGen->SmearGausVertexZ(kTRUE);
+    primGen->SmearVertexXY(kTRUE);
+
+    // Use user defined decays https://fairroot.gsi.de/?q=node/57
+    fRun->SetUserDecay(kTRUE);
+
+    vector<MpdGenerator> Generators;
+    //struct MpdGenerator -> {GeneratorType, FilePath, StartEvent, nEvents}
+    if (GenType == MpdGeneratorType::CUSTOM)
+    {
+        Generators.push_back({MpdGeneratorType::PART, "", 0, 0});
+        Generators.push_back({MpdGeneratorType::ION, "", 0, 0});
+        //and so on... fill vector with all needed generators
+    }
+    else
+        Generators.push_back({GenType, inFile, nStartEvent, nEvents});
+
+    unique_ptr<MpdGeneratorsFactory> genFactory(new MpdGeneratorsFactory());
+    vector<shared_ptr<MpdFactoryMadeGenerator>> fmGens;
+    
+    for (auto Gen : Generators)
+    {
+        
+        switch (Gen.genType)
+        {
+            case MpdGeneratorType::URQMD:
+            case MpdGeneratorType::VHLLE:
+            case MpdGeneratorType::FLUID:
+            case MpdGeneratorType::HSD:
+            case MpdGeneratorType::LAQGSM:
+            {
+                shared_ptr<MpdFactoryMadeGenerator> fmgen = genFactory->create(Gen);
+                if (fmgen.get() == NULL)
+                {
+                    cout << "unable to create generator" << endl;
+                    exit(1);
+                }
+                fmGens.push_back(fmgen->getptr());
+                shared_ptr<FairGenerator> g = fmgen->GetGeneratorPtr();
+                primGen->AddGenerator(g.get());
+            }
+            break;
+            case MpdGeneratorType::PART:
+            {
+                // ------- Particle Generator
+                FairParticleGenerator* partGen = new FairParticleGenerator(211, 10, 1, 0, 3, 1, 0, 0);
+                primGen->AddGenerator(partGen);
+            }
+            break;
+            case MpdGeneratorType::ION:
+            {
+                // ------- Ion Generator
+                FairIonGenerator *fIongen = new FairIonGenerator(79, 197, 79, 1, 0., 0., 25, 0., 0., -1.);
+                primGen->AddGenerator(fIongen);
+            }
+            break;
+            case MpdGeneratorType::BOX:
+            {
+                gRandom->SetSeed(0);
+                // ------- Box Generator
+                FairBoxGenerator* boxGen = new FairBoxGenerator(13, 100); // 13 = muon; 1 = multipl.
+                boxGen->SetPRange(0.25, 2.5); // GeV/c //setPRange vs setPtRange
+                boxGen->SetPhiRange(0, 360); // Azimuth angle range [degree]
+                boxGen->SetThetaRange(0, 180); // Polar angle in lab system range [degree]
+                boxGen->SetXYZ(0., 0., 0.); // mm o cm ??
+                primGen->AddGenerator(boxGen);
+            }
+            break;
+            case MpdGeneratorType::HADGEN:
+            {
+                THadgen* hadGen = new THadgen();
+                hadGen->SetRandomSeed(clock() + time(0));
+                hadGen->SetParticleFromPdgCode(0, 196.9665, 79);
+                hadGen->SetEnergy(6.5E3);
+                MpdGeneralGenerator* generalHad = new MpdGeneralGenerator(hadGen);
+                primGen->AddGenerator(generalHad);
+            }
+            break;
+            case MpdGeneratorType::CUSTOM:
+            default:
+                cout << "unable to find generator type" << endl;
+                exit(1);
+            break;
+        }
+    }
+
+    fRun->SetOutputFile(outFile.Data());
+
+    // Magnetic Field Map - for proper use in the analysis MultiField is necessary here
+    MpdMultiField* fField = new MpdMultiField();
+
+    if (FieldSwitcher == 0) {
+        MpdConstField* fMagField = new MpdConstField();
+        fMagField->SetField(0., 0., 5.); // values are in kG:  1T = 10kG
+        fMagField->SetFieldRegion(-230, 230, -230, 230, -375, 375);
+        fField->AddField(fMagField);
+        fRun->SetField(fField);
+        cout << "FIELD at (0., 0., 0.) = (" <<
+                fMagField->GetBx(0., 0., 0.) << "; " << fMagField->GetBy(0., 0., 0.) << "; " << fMagField->GetBz(0., 0., 0.) << ")" << endl;
+    }
+    else if (FieldSwitcher == 1) {
+        MpdFieldMap* fMagField = new MpdFieldMap("B-field_v2", "A");
+        fMagField->Init();
+        fField->AddField(fMagField);
+        fRun->SetField(fField);
+        cout << "FIELD at (0., 0., 0.) = (" <<
+                fMagField->GetBx(0., 0., 0.) << "; " << fMagField->GetBy(0., 0., 0.) << "; " << fMagField->GetBz(0., 0., 0.) << ")" << endl;
+    }
+
+    fRun->SetStoreTraj(kTRUE);
+    fRun->SetRadLenRegister(flag_store_FairRadLenPoint); // radiation length manager
+
+    //  MpdTpcDigitizerTask* tpcDigitizer = new MpdTpcDigitizerTask();
+    //  tpcDigitizer->SetOnlyPrimary(kTRUE); /// Digitize only primary track
+    //  tpcDigitizer->SetMakeQA(kTRUE);  /// SetMakeQA(kTRUE) prepares Quality Assurance Histograms  
+    //  tpcDigitizer->SetDiffuse(kFALSE);
+    //  tpcDigitizer->SetDebug(kFALSE);
+    //  tpcDigitizer->SetDistort(kFALSE);
+    //  tpcDigitizer->SetResponse(kFALSE);
+    //  tpcDigitizer->SetDistribute(kFALSE);
+    //  fRun->AddTask(tpcDigitizer);
+
+    fRun->Init();
+
+    // -Trajectories Visualization (TGeoManager Only)
+    // Set cuts for storing the trajectories
+    FairTrajFilter* trajFilter = FairTrajFilter::Instance();
+    trajFilter->SetStepSizeCut(0.01); // 1 cm
+    //  trajFilter->SetVertexCut(-2000., -2000., 4., 2000., 2000., 100.);
+    trajFilter->SetMomentumCutP(.50); // p_lab > 500 MeV
+    //  trajFilter->SetEnergyCut(.2, 3.02); // 0 < Etot < 1.04 GeV
+
+    trajFilter->SetStorePrimaries(kTRUE);
+    trajFilter->SetStoreSecondaries(kFALSE);
+
+    // Fill the Parameter containers for this run
+    FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
+
+    Bool_t kParameterMerged = kTRUE;
+    FairParRootFileIo* output = new FairParRootFileIo(kParameterMerged);
+    //AZ output->open(parFile.Data());
+    output->open(gFile);
+    rtdb->setOutput(output);
+
+    MpdMultiFieldPar* Par = (MpdMultiFieldPar*) rtdb->getContainer("MpdMultiFieldPar");
+    if (fField)
+        Par->SetParameters(fField);
+    Par->setInputVersion(fRun->GetRunId(), 1);
+    Par->setChanged();
+    // Par->printParams();
+
+    rtdb->saveOutput();
+    rtdb->print();
+
+    // Transport nEvents
+    fRun->Run(nEvents);
+
+    for (auto fmgen : fmGens)
+    {
+        fmgen->PostActions(fRun);
+    }
+
+    timer.Stop();
+    Double_t rtime = timer.RealTime(), ctime = timer.CpuTime();
+    printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime);
+    cout << "Macro finished successfully." << endl;     // marker of successful execution for CDASH
+}