run_sim_bmn.C

Go to the documentation of this file.

// BM@N simulation: primary particles from an event generator are transported through the setup, and
// then the MC points are converted to detector digits. Output: MC tracks, points and digits in simFile
// genFile - input file with generator data, if needed
// simFile - output file with MC data, default: bmnsim.root
// nStartEvent - start event in the input generator file to begin transporting, default: 0
// nEvents - number of events to transport
// generatorName - generator name for the input file (enumeration above)
// geometryTag - geometry tag from Geometry Database, "run9:" means Run 9 latest geometry (by creation date).
//               The tag can be followed by ';' and a path to an XML file with a local customization of the setup
//               e.g. "run9:;$VMCWORKDIR/my_setup.xml" (see database/geo_db/examples/README.md)
// saveGeoInfo - whether save 'GeoTracks' for Event Display and create "full_geometry.root" with setup geometry
// useRealEffects - whether we use realistic effects at simulation (Lorentz, misalignment)
enum enumGenerators   // enumeration of generator names corresponding input files
{
    UNIGEN,
    URQMD,
    QGSM,
    HSD,
    BOX,
    PART,
    ION,
    DCMQGSM,
    DCMSMM
};
 
void run_sim_bmn(TString genFile = "",
                 TString simFile = "$VMCWORKDIR/macro/run/bmnsim.root",
                 Int_t nStartEvent = 0,
                 Int_t nEvents = 10,
                 enumGenerators generatorName = BOX,
                 TString geometryTag = "run9:",
                 Bool_t saveGeoInfo = kFALSE,
                 Bool_t useRealEffects = kTRUE)
{
    gDebug = 0;   // ROOT debug option
    // if (gEnv) gEnv->SetValue("Root.Stacktrace", 0);
    //  Verbosity level (0 = summary only, 1 = event-level, 2 = track-level, 3 = full debug)
    Int_t iVerbose = 0;
    // LLog level for LOG(severity): "fatal", "error", "alarm", "important", "warn", "state", "info", "detail",
    // "debug", "debug[1-4]", "trace", "nolog"
    FairLogger::GetLogger()->SetLogScreenLevel("info");
    gRandom->SetSeed(0);   // set seed for PRNG
 
    TStopwatch timer;
    timer.Start();
 
    // -----   Create simulation run   ----------------------------------------
    FairRunSim* fRun = new FairRunSim();
    fRun->SetUseFairLinks(kTRUE);
 
    // Choose the Geant Navigation System
    fRun->SetName("TGeant4");   // or "TGeant3"
 
    // The geometry tag can carry a path to an XML file with a local customization of the setup
    // (see database/geo_db/examples/README.md) after the ';' separator, e.g.:
    //   "run9:;my_setup.xml"        - the latest Run 9 geometry customized by the file
    TString geoLocalSettings = "";
    if (BmnFunctionSet::ParseGeometryTag(geometryTag, geoLocalSettings) != 0)
        exit(-1);
 
    // Get period from the given geometry tag and set FaiRun::RunId
    UInt_t run_period = 9;
    TPRegexp re_run_number("(run(\\d+))");
    TString run_tag = geometryTag;
    if (re_run_number.MatchB(run_tag)) {
        re_run_number.Substitute(run_tag, "$2");
        run_period = run_tag.Atoi();
        fRun->SetRunId(run_period);
    }
    // Load the BM@N setup geometry from the Geometry Database
    GeoSetup gSetup;
    if (TString(geometryTag.Strip(TString::kBoth)).EndsWith(":")) {
        geometryTag = TString(geometryTag.Strip(TString::kBoth)).Chop();   // remove trailing ":" if latest geometry
        geometryTag = gSetup.GetLastApprovedSetupTag(geometryTag.Data()).c_str();
    }
    Int_t errStatus =
        gSetup.LoadSetupToFairRunSim(geometryTag.Data(), geoLocalSettings.IsNull() ? nullptr : geoLocalSettings.Data());
    if (errStatus != 0) {
        cout << "ERROR: Geometry not loaded. Error status " << errStatus << endl;
        exit(-1);
    }
 
    // Use the experiment specific MC Event header instead of the default one
    // This one stores additional information about the reaction plane
    // MpdMCEventHeader* mcHeader = new MpdMCEventHeader();
    // fRun->SetMCEventHeader(mcHeader);
 
    // Create and Set Event Generator
    FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
    fRun->SetGenerator(primGen);
 
    // Smearing of beam interaction point, if needed, and primary vertex position
    // DO NOT do it in corresponding gen. sections to avoid incorrect summation!!!
    // parameters for period 8
    primGen->SetBeam(0.4376, 0.1239 - 0.3, 0.4537, 0.6925);   // (beamX0, beamY0, beamSigmaX, beamSigmaY)
    primGen->SetBeamAngle(0.011, 0.0, 0.0, 0.0);   // (beamAngleX0, beamAngleY0, beamAngleSigmaX, beamAngleSigmaY)
    primGen->SetTarget(-0.01942, 0.175);           // (targetZ, targetDz)
    primGen->SmearVertexZ(kTRUE);
    primGen->SmearGausVertexXY(kTRUE);
 
    switch (generatorName) {
        case UNIGEN: {
            if (!BmnFunctionSet::CheckFileExist(genFile, 1))
                exit(-1);
 
            auto* unigen = new MpdUnigenGenerator(genFile);
            unigen->RegisterIons();
            primGen->AddGenerator(unigen);
            unigen->SetEventPlane(-M_PI, M_PI);
            auto n_events_in_tree = unigen->GetNEntries();
            if (nStartEvent >= n_events_in_tree) {
                cout << "You are trying to specify the event out reach of current tree" << generatorName << endl;
                exit(-3);
            }
            auto n_avail_events = n_events_in_tree - nStartEvent - 1;
            if (nEvents <= 0)
                nEvents = n_avail_events;
            else
                nEvents = std::min<int>(n_avail_events, nEvents);
            cout << nEvents << " are availiable and will be processed" << endl;
            if (nStartEvent > 0)
                unigen->SkipEvents(nStartEvent);
            break;
        }
            // ------- UrQMD Generator
        case URQMD: {
            if (!BmnFunctionSet::CheckFileExist(genFile, 1))
                exit(-1);
 
            MpdUrqmdGenerator* urqmdGen = new MpdUrqmdGenerator(genFile);
            // urqmdGen->SetEventPlane(0., 360.);
            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(genFile.Data()) - nStartEvent;
            break;
        }
 
            // ------- Particle Generator
        case PART: {
            // FairParticleGenerator generates a single particle type (PDG, mult, [GeV] px, py, pz, [cm] vx, vy, vz)
            FairParticleGenerator* partGen = new FairParticleGenerator(211, 10, 1, 0, 3, 1, 0, 0);
            primGen->AddGenerator(partGen);
            break;
        }
 
            // ------- Ion Generator
        case ION: {
            // Start beam from a far point to check mom. reconstruction procedure (Z, A, q, mult, [GeV] px, py, pz, [cm]
            // vx, vy, vz)
            FairIonGenerator* fIongen = new FairIonGenerator(54, 124, 54, 1, 0., 0., 4.65, 0., 0., 0.);   // Xe
            primGen->AddGenerator(fIongen);
            break;
        }
 
            // ------- Box Generator
        case BOX: {
            FairBoxGenerator* boxGen = new FairBoxGenerator(2212, 1);   //(PDG code, multiplicity)
            boxGen->SetPRange(0.2, 5.0);                                // GeV/c, setPRange vs setPtRange
            boxGen->SetPhiRange(0, 360);                                // Azimuth angle range [degree]
            boxGen->SetThetaRange(0, 40.0);                             // Polar angle in lab system range [degree]
            primGen->AddGenerator(boxGen);
            break;
        }
 
            // ------- HSD/PHSD Generator
        case HSD: {
            if (!BmnFunctionSet::CheckFileExist(genFile, 1))
                exit(-1);
 
            MpdPHSDGenerator* hsdGen = new MpdPHSDGenerator(genFile.Data());
            // hsdGen->SetPsiRP(0.); // set fixed Reaction Plane angle 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(genFile.Data()) - nStartEvent;
            break;
        }
 
            // ------- LAQGSM/DCM-QGSM Generator
        case QGSM:
        case DCMQGSM: {
 
            if (!BmnFunctionSet::CheckFileExist(genFile, 1))
                exit(-1);
 
            MpdLAQGSMGenerator* guGen = new MpdLAQGSMGenerator(genFile.Data(), kFALSE);
            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(genFile.Data()) - nStartEvent;
            break;
        }
        case DCMSMM: {
 
            if (!BmnFunctionSet::CheckFileExist(genFile, 1))
                exit(-1);
 
            MpdDCMSMMGenerator* smmGen = new MpdDCMSMMGenerator(genFile.Data());
            primGen->AddGenerator(smmGen);
            if (nStartEvent > 0)
                smmGen->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::GetNumDCMSMMEvents(genFile.Data()) - nStartEvent;
            break;
        }
 
        default: {
            cout << "ERROR: Generator name was not pre-defined: " << generatorName << endl;
            exit(-3);
        }
    }   // end of switch (generatorName)
 
    if (simFile == "")
        simFile = "$VMCWORKDIR/macro/run/bmnsim.root";
    // if directory for the output file does not exist, then create
    if (BmnFunctionSet::CreateDirectoryTree(simFile, 1) < 0)
        exit(-2);
    fRun->SetSink(new FairRootFileSink(simFile.Data()));
    fRun->SetIsMT(false);
 
    // -----   Define beam, target and energy by the generator filename   -----------
    TString beamParticle = "", targetParticle = "";
    Double_t beamEnergy = -1;   // beam energy (GeV)
    BmnFunctionSet::ParseFileName(gSystem->BaseName(genFile), beamParticle, targetParticle, beamEnergy);
 
    // the parsed values are written to the FairBaseParSet container of the simulation file: the beam energy - to the
    // beam momentum field, the beam and target names - to the container description as 'beam;target'.
    // At the reconstruction stage BmnFillDstTask takes them from there to fill the run header of the DST file
    if (beamEnergy > 0)
        fRun->SetBeamMom(beamEnergy);   // FairRunSim writes it to FairBaseParSet during the initialization
    fRun->GetRuntimeDb()
        ->getContainer("FairBaseParSet")
        ->setDescription(TString::Format("%s;%s", beamParticle.Data(), targetParticle.Data()));
 
    // -----   Create magnetic field   ----------------------------------------
    Double_t fieldScale = 0.929;   // default for 3.8 GeV: 0.929 for 1900 A map (112.0 mV);
    // trying to define magnetic field by the beam energy
    if (beamEnergy > 0) {
        double fieldValue = BmnFunctionSet::GetMagneticFieldByEnergy(beamParticle, beamEnergy);
        if (fieldValue > 0)
            fieldScale = fieldValue / 111.75;   //? 112.0 gives 0.926 for 3.8 GeV, not 0.929
    }
    BmnFieldMap* magField = nullptr;
    // if GetMagneticField(-1) -> get default value from the Geometry Database
    magField = dynamic_cast<BmnFieldMap*>(gSetup.GetMagneticField(fieldScale));
    // magField = new BmnNewFieldMap("FieldMap_1900_extrap_noPed.root");
    // magField->SetScale(fieldScale);
    fRun->SetField(magField);
 
    fRun->SetStoreTraj(kTRUE);
    fRun->SetRadLenRegister(kFALSE);   // radiation length manager
 
    // -----   Digitizers: converting MC points to detector digits   -----------
 
    // SiBT-Digitizer
    BmnSiBTDigitizer* sibtDigit = new BmnSiBTDigitizer(run_period);
    fRun->AddTask(sibtDigit);
 
    // SiProf-Digitizer
    BmnSiProfDigitizer* siprofDigit = new BmnSiProfDigitizer(run_period);
    fRun->AddTask(siprofDigit);
 
    // BD-Digitizer
    BmnBdDigitizer* bdDigit = new BmnBdDigitizer();
    fRun->AddTask(bdDigit);
 
    // SiMD-Digitizer
    BmnSiMDDigitizer* simdDigit = new BmnSiMDDigitizer();
    fRun->AddTask(simdDigit);
 
    // VSP-Digitizer
    BmnVSPDigitizer* vspDigit = new BmnVSPDigitizer(run_period);
    vspDigit->SetUseRealEffects(useRealEffects);
    vspDigit->SetNoise(45.0, 150.0);
    fRun->AddTask(vspDigit);
 
    // SI-Digitizer
    BmnSiliconDigitizer* siliconDigit = new BmnSiliconDigitizer(run_period);
    siliconDigit->SetUseRealEffects(useRealEffects);
    siliconDigit->SetNoise(45.0, 150.0);   // AZ-150224
    fRun->AddTask(siliconDigit);
 
    // GEM-Digitizer
    BmnGemStripDigitizer* gemDigit = new BmnGemStripDigitizer(run_period);
    if (useRealEffects) {
        BmnGemStripMedium::GetInstance().SetCurrentConfiguration(
            // BmnGemStripMediumConfiguration::ARC4H10_80_20_E_1720_2240_3230_3730_B_0_7T);   // 3.8 AGeV
            // BmnGemStripMediumConfiguration::ARC4H10_80_20_E_1720_2240_3230_3730_B_0_3T);   // 1.6 AGeV
            // BmnGemStripMediumConfiguration::ARC4H10_80_20_E_1720_2240_3230_3730_B_0_4T);   // 2.2 AGeV
            BmnGemStripMediumConfiguration::ARC4H10_80_20_E_1720_2240_3230_3730_B_0_6T);   // 3.0 AGeV
        gemDigit->SetUseRealEffects(useRealEffects);
    }
    gemDigit->SetNoise(10.0, 25.0);   // AZ-150224
    BmnGemStripModule::SetUseANNGeneration(kFALSE);
    // BmnGemStripModule::SetANNModelPath("$VMCWORKDIR/parameters/gem/ANNModels/model_21.bmnmlp");
    fRun->AddTask(gemDigit);
 
    // CSC-Digitizer
    BmnCSCDigitizer* cscDigit = new BmnCSCDigitizer(run_period);
    fRun->AddTask(cscDigit);
 
    // // ZDC-Digitizer (for period < 8)
    // BmnZdcDigitizer* zdcDigit = new BmnZdcDigitizer();
    // zdcDigit->SetScale(39e3);
    // zdcDigit->SetThreshold(500.);
    // fRun->AddTask(zdcDigit);
 
    // FHCal-Digitizer (for period >= 8)
    BmnFHCalDigitizer* fhcalDigit = new BmnFHCalDigitizer(run_period);
    fRun->AddTask(fhcalDigit);
 
    // ScWall-Digitizer
    BmnScWallDigitizer* fScWallDigit = new BmnScWallDigitizer(run_period);
    fRun->AddTask(fScWallDigit);
 
    // Hodo-Digitizer
    BmnHodoDigitizer* fHodoDigit = new BmnHodoDigitizer(run_period);
    fRun->AddTask(fHodoDigit);
 
    // Ndet-Digitizer
    BmnNdetDigitizer* fNdetDigit = new BmnNdetDigitizer(run_period);
    fRun->AddTask(fNdetDigit);
 
    // // ECAL-Digitizer (for period < 8)
    // BmnEcalDigitizer* ecalDigit = new BmnEcalDigitizer();
    // fRun->AddTask(ecalDigit);
 
    fRun->SetStoreTraj(saveGeoInfo);
    fRun->Init();
    magField->Print("");
 
    // Cuts for Trajectories Visualization (TGeoManager only)
    if (saveGeoInfo) {
        FairTrajFilter* trajFilter = FairTrajFilter::Instance();
        // Set cuts for storing the trajectories
        trajFilter->SetStepSizeCut(0.1);   // [cm]
        trajFilter->SetVertexCut(-200., -200., -200., 200., 200.,
                                 1100.);      // (vxMin, vyMin, vzMin, vxMax, vyMax, vzMax)
        trajFilter->SetMomentumCutP(50e-3);   // p_lab > 50 MeV
        // trajFilter->SetEnergyCut(0., 60.);  // 0 < Etot < 60 GeV
        trajFilter->SetStorePrimaries(kTRUE);
        trajFilter->SetStoreSecondaries(kTRUE);   // 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);
 
    // the description keeps the full geometry tag: the setup tag and the XML customization file (if any) after ';'
    TString geoDescription = geometryTag;
    if (!geoLocalSettings.IsNull())
        geoDescription += ";" + geoLocalSettings;
    rtdb->getContainer("FairGeoParSet")->setDescription(geoDescription.Data());
    rtdb->saveOutput();
    rtdb->print();
 
    // Transport nEvents
    fRun->Run(nEvents);
 
    // gGeoManager->CheckOverlaps(0.0001);
    // gGeoManager->PrintOverlaps();
    if (saveGeoInfo) {
        fRun->CreateGeometryFile("full_geometry.root");   // save the full setup geometry to the additional file
    }
 
    if ((generatorName == QGSM) || (generatorName == DCMQGSM)) {
        TString Pdg_table_name = TString::Format("%s%s%c%s", gSystem->BaseName(genFile.Data()), ".g",
                                                 (fRun->GetName())[6], ".pdg_table.dat");
        (TDatabasePDG::Instance())->WritePDGTable(Pdg_table_name.Data());
    }
    delete fRun;
    delete magField;
 
    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 successfully execution for software testing systems
}