#if !defined(__CINT__) && !defined(__CLING__) #include "TString.h" #include "TStopwatch.h" #include "TROOT.h" #include "TSystem.h" #include "TMath.h" #include "TDatabasePDG.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 using namespace std; #endif R__ADD_INCLUDE_PATH($VMCWORKDIR) #include "macro/mpd/mpdloadlibs.C" #include "macro/mpd/geometry_stage1.C" //#include "macro/mpd/geometry_v2.C" #define HSD // Choose generator: URQMD VHLLE FLUID PART ION BOX HSD LAQGSM HADGEN #define GEANT3 // Choose: GEANT3 GEANT4 // 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(TString inFile = "$VMCWORKDIR/macro/mpd/MST.dat.gz", TString outFile = "$VMCWORKDIR/macro/mpd/evetest.root", Int_t nStartEvent = 0, Int_t nEvents = 500, Bool_t flag_store_FairRadLenPoint = kFALSE, Int_t FieldSwitcher = 0) { TStopwatch timer; timer.Start(); gDebug = 0; FairRunSim* fRun = new FairRunSim(); // Choose the Geant Navigation System #ifdef GEANT3 fRun->SetName("TGeant3"); #else fRun->SetName("TGeant4"); #endif 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); // Use external decayer //fRun->SetPythiaDecayer(TString("$VMCWORKDIR/gconfig/LambdaDecayConfig.C")); #ifdef MCDST // McDst generator if (!CheckFileExist(inFile)) return; MpdMcDstGenerator* mcDstGen = new MpdMcDstGenerator(inFile); primGen->AddGenerator(mcDstGen); #else #ifdef URQMD // <---- Urqmd Generator if (!CheckFileExist(inFile)) return; MpdUrqmdGenerator* urqmdGen = new MpdUrqmdGenerator(inFile); // Event plane angle (in degrees) will be generated by uniform distribution from min to max Float_t min = 0.0, max = 30.0; urqmdGen->SetEventPlane(min * TMath::DegToRad(), max * TMath::DegToRad()); primGen->AddGenerator(urqmdGen); 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(inFile.Data()) - nStartEvent; #else #ifdef VHLLE if (!CheckFileExist(inFile)) return; MpdVHLLEGenerator* vhlleGen = new MpdVHLLEGenerator(inFile, kTRUE); // kTRUE corresponds to hydro + cascade, kFALSE -- hydro only vhlleGen->SkipEvents(0); primGen->AddGenerator(vhlleGen); #else #ifdef FLUID if (!CheckFileExist(inFile)) return; 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 primGen->AddGenerator(fluidGen); #else #ifdef PART // <---- Particle Generator FairParticleGenerator* partGen = new FairParticleGenerator(211, 10, 1, 0, 3, 1, 0, 0); primGen->AddGenerator(partGen); #else #ifdef ION // <---- Ion Generator FairIonGenerator *fIongen = new FairIonGenerator(79, 197, 79, 1, 0., 0., 25, 0., 0., -1.); primGen->AddGenerator(fIongen); #else #ifdef BOX // <---- Box Generator gRandom->SetSeed(0); 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); #else #ifdef HSD // <---- HSD/PHSD Generator if (!CheckFileExist(inFile)) return; MpdPHSDGenerator *hsdGen = new MpdPHSDGenerator(inFile.Data()); //hsdGen->SetPsiRP(0.); // set fixed Reaction Plane angle [rad] instead of random primGen->AddGenerator(hsdGen); 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(inFile.Data()) - nStartEvent; #else #ifdef LAQGSM // <---- LAQGSM Generator if (!CheckFileExist(inFile)) return; MpdLAQGSMGenerator* guGen = new MpdLAQGSMGenerator(inFile.Data(), kTRUE, 4, 1+nStartEvent+nEvents); // kTRUE - for NICA/MPD, 1+nStartEvent+nEvents - search ions in selected part of file. primGen->AddGenerator(guGen); 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(inFile.Data()) - nStartEvent; #else #ifdef 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); #endif #endif #endif #endif #endif #endif #endif #endif #endif #endif 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(); //AZ - Enable decays of unstable particles /* MpdStack *stack = NULL; if (TString(fRun->GetName()).Contains("TGeant3")) stack = (MpdStack*) ((TGeant3*)gMC)->GetStack(); else stack = (MpdStack*) ((TGeant4*)gMC)->GetStack(); stack->SetDecayUnstable(); */ // -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); #ifdef LAQGSM TString Pdg_table_name = TString::Format("%s%s%c%s", gSystem->BaseName(inFile.Data()), ".g", (fRun->GetName())[6], ".pdg_table.dat"); (TDatabasePDG::Instance())->WritePDGTable(Pdg_table_name.Data()); #endif 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 }