run_ana_bmn.C

Go to the documentation of this file.

// Macro for physics analysis of simulated and experimental BM@N events
// dstFileName - input DST file with data (MC or exp. data), a ROOT file with a TChain object of DST files or
//               text file with a list of DST files; run info, mag.field and pars are taken from the first file
// anaFileName - output file with analysis results
// nStartEvent - number of first event to process (starts with zero), default: 0
// nEvents - number of events to process after nStartEvent, 0 - all events of the given file(s) after nStartEvent
void run_ana_bmn(TString dstFileName = "$VMCWORKDIR/macro/run/bmndst.root",
                 TString anaFileName = "$VMCWORKDIR/macro/run/bmnana.root",
                 Int_t nStartEvent = 0,
                 Int_t nEvents = 10)
{
    gDebug = 0;   // ROOT debug option
    // Verbosity level (0 = summary only, 1 = event-level, 2 = track-level, 3 = full debug)
    Int_t iVerbose = 0;
    // Log level for LOG(severity): "fatal", "error", "alarm", "important", "warn", "state", "info", "detail", "debug",
    // "debug[1-4]", "trace", "nolog"
    FairLogger::GetLogger()->SetLogScreenLevel("info");
 
    // -----   Timer   ---------------------------------------------------------
    TStopwatch timer;
    timer.Start();
 
    // get the list of input DST files (input file can also contain a TChain or a text list of DST files)
    vector<TString> dstFileList;
    if (BmnFunctionSet::GetInputFileList(dstFileName, dstFileList) < 0)
        exit(-1);
    // the first input file is used to get the run information, magnetic field and parameters
    TString firstFileName = dstFileList[0];
 
    // -----   Analysis run   --------------------------------------------
    FairRunAna* fRunAna = new FairRunAna();
    // fRunAna->SetEventHeader(new AnaEventHeader());
 
    // 'isExp' flag for FairRunAna tasks whether simulation or experimental data are used
    Bool_t isExp = (BmnFunctionSet::isSimulationFile(firstFileName) == 0);
 
    // get field value and run information (beam, target, energy) from the DST file
    Int_t fieldType = -1;
    Double_t fieldScale = 0., beam_energy = 0.;
    TString beam = "", target = "";
    FairField* magField =
        BmnFunctionAna::GetDSTField(firstFileName, fieldScale, fieldType, &beam, &target, &beam_energy);
    if (magField != nullptr)
        fRunAna->SetField(magField);
 
    // Declare input source as simulation file or experimental data
    Int_t run_period = -1, run_number = -1;
    BmnFileSource* fFileSource = new BmnFileSource(firstFileName, run_period, run_number);
    // the rest input files are chained to the first one (their branch structure must be the same)
    for (size_t i = 1; i < dstFileList.size(); i++)
        fFileSource->AddFile(dstFileList[i]);
    fRunAna->SetSource(fFileSource);
    // parameters are read only from the first file, so do not reinitialize them when the run ID changes in the chain
    if (dstFileList.size() > 1)
        fRunAna->SetContainerStatic(kTRUE);
 
    // print the run summary for both experimental and simulation data
    TString title = isExp ? "EXPERIMENTAL RUN SUMMARY" : "SIMULATION RUN SUMMARY";
    Int_t nSideBars = (58 - title.Length() - 2) / 2;   // the summary box is 58 symbols wide
    cout << "\n\n"
         << string(nSideBars, '|') << " " << title << " " << string(58 - nSideBars - title.Length() - 2, '|') << endl;
    cout << "||\t\t\t\t\t\t\t||" << endl;
    cout << left;
    cout << "||\t\tPeriod:\t\t" << setw(24) << run_period << "||" << endl;
    if (isExp)
        cout << "||\t\tNumber:\t\t" << setw(24) << run_number << "||" << endl;
    cout << "||\t\tBeam:\t\t" << setw(24) << (beam == "" ? "-" : beam.Data()) << "||" << endl;
    cout << "||\t\tEnergy:\t\t" << setw(24)
         << ((beam_energy > 0) ? TString::Format("%.2f AGeV", beam_energy).Data() : "-") << "||" << endl;
    cout << "||\t\tTarget:\t\t" << setw(24) << (target == "" ? "-" : target.Data()) << "||" << endl;
    cout << "||\t\tField scale:\t" << setw(24) << setprecision(4) << fieldScale << "||" << endl;
    if (dstFileList.size() > 1)
        cout << "||\t\tInput files:\t" << setw(24) << dstFileList.size() << "||" << endl;
    cout << "||\t\t\t\t\t\t\t||" << endl;
    cout << string(58, '|') << "\n\n" << endl;
 
    // if directory for the output file does not exist, then create
    if (BmnFunctionSet::CreateDirectoryTree(anaFileName, 1) < 0)
        exit(-2);
    fRunAna->SetSink(new FairRootFileSink(anaFileName));
    fRunAna->SetGenerateRunInfo(false);
 
    // if nEvents is equal 0 then all events of the given file(s) starting with "nStartEvent" should be processed
    if (nEvents == 0) {
        for (auto& fileName : dstFileList) {
            Int_t nFileEvents = MpdGetNumEvents::GetNumROOTEvents(fileName.Data());
            if (nFileEvents < 0)
                exit(-3);
            nEvents += nFileEvents;
        }
        nEvents -= nStartEvent;
    }
 
    // Add TObjectString file names to a TList which is passed as input to the FairParAsciiFileIo
    // The FairParAsciiFileIo creates a concatenated input parameter file all_<pid>.par, which is then can be used
    // TList* parFileNameList = new TList();
 
    // ====================================================================== //
    // ===                          Analysis Tasks                        === //
    // ====================================================================== //
    CbmKF* kalman_filter = new CbmKF();
    fRunAna->AddTask(kalman_filter);
 
    FairTask* ana_task = new BmnAnalysisTask();
    fRunAna->AddTask(ana_task);
 
    // -----   Parameter database (input)  -------------------------------------
    FairRuntimeDb* rtdb = fRunAna->GetRuntimeDb();
    FairParRootFileIo* parIo1 = new FairParRootFileIo();
    // FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
    parIo1->open(firstFileName.Data());
    // parIo2->open(parFileNameList, "in");
    rtdb->setFirstInput(parIo1);
    // rtdb->setSecondInput(parIo2);
    // -------------------------------------------------------------------------
 
    // -----   Initialize FaiRunAna --------------------------------------------
    fRunAna->GetMainTask()->SetVerbose(iVerbose);
    fRunAna->Init();   // geometry is loaded from DST FairGeoParSet container
    // -------------------------------------------------------------------------
 
    // -----   Parameter database (output)  ------------------------------------
    FairParRootFileIo* parOut = new FairParRootFileIo(kTRUE);
    parOut->open(gFile);
    rtdb->setOutput(parOut);
    rtdb->writeContainer(rtdb->getContainer("FairGeoParSet"), rtdb->getCurrentRun());
    rtdb->saveOutput();
    // rtdb->print();
    // -------------------------------------------------------------------------
 
    // -----   Run FairRunAna  -------------------------------------------------
    fRunAna->Run(nStartEvent, nStartEvent + nEvents);
    // -------------------------------------------------------------------------
 
    // -----   Finish   --------------------------------------------------------
    delete fRunAna;
    timer.Stop();
    Double_t rtime = timer.RealTime();
    Double_t ctime = timer.CpuTime();
    cout << endl << endl;
    cout << "Macro finished successfully." << endl;   // marker of successful execution
    if (dstFileList.size() > 1)
        cout << "Input  file list is " << dstFileName << " (" << dstFileList.size() << " files)" << endl;
    else
        cout << "Input  file is " + firstFileName << endl;
    cout << "Output file is " + anaFileName << endl;
    cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
    cout << endl;
    // ------------------------------------------------------------------------
}