From 915a7021c5867abd127afc46eda7de98b32c5127 Mon Sep 17 00:00:00 2001 From: Cornelius Feier-Riesen <cornelius.riesen@physik.uni-giessen.de> Date: Thu, 22 Jul 2021 17:00:12 +0200 Subject: [PATCH] Update run_reco.C macro modified: analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx modified: macro/analysis/dielectron/batch_job.py modified: macro/analysis/dielectron/batch_send_common.py modified: macro/analysis/dielectron/draw_all.py modified: macro/analysis/dielectron/hadd_many.py modified: macro/analysis/dielectron/run_analysis.C modified: macro/analysis/dielectron/run_litqa.C modified: macro/analysis/dielectron/run_reco.C --- .../lmvm/CbmAnaDielectronTaskDrawAll.cxx | 4 +- macro/analysis/dielectron/batch_job.py | 6 +- .../analysis/dielectron/batch_send_common.py | 2 - macro/analysis/dielectron/draw_all.py | 4 +- macro/analysis/dielectron/hadd_many.py | 10 +- macro/analysis/dielectron/run_analysis.C | 2 +- macro/analysis/dielectron/run_litqa.C | 4 +- macro/analysis/dielectron/run_reco.C | 470 ++++++++++++++++-- 8 files changed, 440 insertions(+), 62 deletions(-) mode change 100644 => 100755 macro/analysis/dielectron/run_reco.C diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx index 21b2d45e70..c7043d40e9 100644 --- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx @@ -1686,7 +1686,9 @@ void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos() // from 'Cocktail + BG' - TH1D* cock = NULL; + + TH1D* cock = nullptr; + if (step == kMc) cock = (TH1D*) GetCoctailMinv(kMc); else if (step == kAcc) cock = (TH1D*) GetCoctailMinv(kAcc); diff --git a/macro/analysis/dielectron/batch_job.py b/macro/analysis/dielectron/batch_job.py index b61b1c4caa..826203a3c5 100755 --- a/macro/analysis/dielectron/batch_job.py +++ b/macro/analysis/dielectron/batch_job.py @@ -41,11 +41,11 @@ def main(): litQaFile = dataDir + "/litqa." + taskId + ".root" analysisFile = dataDir + "/analysis." + taskId + ".root" geoSimFile = dataDir + "/geosim." + taskId + ".root" - print("Here: in batchJobPy") + sEvBuildRaw = "Ideal" #os.system(('root -l -b -q {}/run_sim.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, urqmdFile, plutoFile, mcFile, parFile, geoSimFile, geoSetup, nofEvents)) #os.system(('root -l -b -q {}/run_digi.C\(\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, mcFile, parFile, digiFile, nofEvents)) - #os.system(('root -l -b -q {}/run_reco.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, mcFile, parFile, digiFile, recoFile, geoSetup, nofEvents)) - #os.system(('root -l -b -q {}/run_litqa.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, mcFile, parFile, digiFile, recoFile, litQaFile, geoSetup, nofEvents)) + os.system(('root -l -b -q {}/run_reco.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, mcFile, parFile, digiFile, recoFile, geoSetup, sEvBuildRaw, nofEvents)) + os.system(('root -l -b -q {}/run_litqa.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, mcFile, parFile, digiFile, recoFile, litQaFile, geoSetup, nofEvents)) os.system(('root -l -b -q {}/run_analysis.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, mcFile, parFile, digiFile, recoFile, analysisFile, plutoParticle, colSystem, colEnergy, geoSetup, nofEvents)) def getPlutoPath(colSystem, colEnergy, plutoParticle, taskId): diff --git a/macro/analysis/dielectron/batch_send_common.py b/macro/analysis/dielectron/batch_send_common.py index a76bca4e44..10352ace63 100755 --- a/macro/analysis/dielectron/batch_send_common.py +++ b/macro/analysis/dielectron/batch_send_common.py @@ -1,7 +1,5 @@ #!/usr/bin/env python3 -## #changed on Feb 8th; first check if works or not ### - import os import shutil diff --git a/macro/analysis/dielectron/draw_all.py b/macro/analysis/dielectron/draw_all.py index aa6a559486..e5847052d1 100755 --- a/macro/analysis/dielectron/draw_all.py +++ b/macro/analysis/dielectron/draw_all.py @@ -19,8 +19,8 @@ def main(): drawSig = True os.system(('root -l -b -q {}/draw_analysis.C\(\\"{}\\",\\"{}\\"\)').format(macroDir, inFilesAna, outFile, uMvd, drawSig)) - #inFilesLitqa = dataDir + plutoParticle + "/litqa.all.root" - #os.system(('root -l -b -q {}/draw_litqa.C\(\\"{}\\",\\"{}\\"\)').format(macroDir, inFilesLitqa, outFile)) + inFilesLitqa = dataDir + plutoParticle + "/litqa.all.root" + os.system(('root -l -b -q {}/draw_litqa.C\(\\"{}\\",\\"{}\\"\)').format(macroDir, inFilesLitqa, outFile)) print("===== DRAW ALL =====") os.system(('root -l -b -q {}/draw_analysis_all.C').format(macroDir)) diff --git a/macro/analysis/dielectron/hadd_many.py b/macro/analysis/dielectron/hadd_many.py index 1ac5f91805..e65ba53bd1 100755 --- a/macro/analysis/dielectron/hadd_many.py +++ b/macro/analysis/dielectron/hadd_many.py @@ -8,11 +8,11 @@ def main(): for plutoParticle in plutoParticles: dataDirPluto = dataDir + "/" + plutoParticle - #inFilesLitQa = dataDirPluto + "/litqa.*.root" - #outFileLitQa = dataDirPluto + "/litqa.all.root" - #if os.path.exists(outFileLitQa): - # os.remove(outFileLitQa) - #os.system(('hadd -j -T -f {} {}').format(outFileLitQa, inFilesLitQa)) + inFilesLitQa = dataDirPluto + "/litqa.*.root" + outFileLitQa = dataDirPluto + "/litqa.all.root" + if os.path.exists(outFileLitQa): + os.remove(outFileLitQa) + os.system(('hadd -j -T -f {} {}').format(outFileLitQa, inFilesLitQa)) inFilesAna = dataDirPluto + "/analysis.*.root" outFileAna = dataDirPluto + "/analysis.all.root" diff --git a/macro/analysis/dielectron/run_analysis.C b/macro/analysis/dielectron/run_analysis.C index d5718321c7..5464c5a4a3 100755 --- a/macro/analysis/dielectron/run_analysis.C +++ b/macro/analysis/dielectron/run_analysis.C @@ -73,7 +73,7 @@ void run_analysis(const string& mcFile = "/lustre/nyx/cbm/users/criesen/c FairRuntimeDb* rtdb = run->GetRuntimeDb(); FairParRootFileIo* parIo1 = new FairParRootFileIo(); FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo(); - parIo1->open(parFile.c_str(), "READ"); + parIo1->open(parFile.c_str(), "UPDATE"); rtdb->setFirstInput(parIo1); if (!parFileList->IsEmpty()) { parIo2->open(parFileList, "in"); diff --git a/macro/analysis/dielectron/run_litqa.C b/macro/analysis/dielectron/run_litqa.C index 7a3020d00b..2c46dd2e4e 100644 --- a/macro/analysis/dielectron/run_litqa.C +++ b/macro/analysis/dielectron/run_litqa.C @@ -65,7 +65,7 @@ void run_litqa(const string& mcFile = "/lustre/nyx/cbm/users/criesen/cbm/data/ // RICH reco QA CbmRichRecoQa* richRecoQa = new CbmRichRecoQa(); richRecoQa->SetOutputDir(""); - // run->AddTask(richRecoQa); + //run->AddTask(richRecoQa); // Reconstruction Qa CbmLitTrackingQa* trackingQa = new CbmLitTrackingQa(); @@ -112,7 +112,7 @@ void run_litqa(const string& mcFile = "/lustre/nyx/cbm/users/criesen/cbm/data/ FairRuntimeDb* rtdb = run->GetRuntimeDb(); FairParRootFileIo* parIo1 = new FairParRootFileIo(); FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo(); - parIo1->open(parFile.c_str(), "READ"); + parIo1->open(parFile.c_str(), "UPDATE"); rtdb->setFirstInput(parIo1); if (!parFileList->IsEmpty()) { parIo2->open(parFileList, "in"); diff --git a/macro/analysis/dielectron/run_reco.C b/macro/analysis/dielectron/run_reco.C old mode 100644 new mode 100755 index 796fd207c6..09180a8480 --- a/macro/analysis/dielectron/run_reco.C +++ b/macro/analysis/dielectron/run_reco.C @@ -1,22 +1,119 @@ -/* Copyright (C) 2011-2020 Justus-Liebig-Universitaet Giessen, Giessen +/* Copyright (C) 2020-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt SPDX-License-Identifier: GPL-3.0-only - Authors: Andrey Lebedev, Semen Lebedev, Elena Lebedeva [committer] */ + Authors: Volker Friese [committer], Dominik Smith */ + +/** @file run_reco.C + ** @author Volker Friese <v.friese@gsi.de> + ** @since 14 November 2020 + **/ + + +// --- Includes needed for IDE +#include <RtypesCore.h> +#if !defined(__CLING__) +#include "CbmBuildEventsFromTracksReal.h" +#include "CbmBuildEventsIdeal.h" +#include "CbmBuildEventsQa.h" +#include "CbmDefs.h" +#include "CbmFindPrimaryVertex.h" +#include "CbmKF.h" +#include "CbmL1.h" +#include "CbmL1StsTrackFinder.h" +#include "CbmLitFindGlobalTracks.h" +#include "CbmMCDataManager.h" +#include "CbmMatchRecoToMC.h" +#include "CbmMuchFindHitsGem.h" +#include "CbmMvdClusterfinder.h" +#include "CbmMvdHitfinder.h" +#include "CbmPVFinderKF.h" +#include "CbmPrimaryVertexFinder.h" +#include "CbmPsdHitProducer.h" +#include "CbmRecoSts.h" +#include "CbmRichHitProducer.h" +#include "CbmRichReconstruction.h" +#include "CbmSetup.h" +#include "CbmStsFindTracks.h" +#include "CbmStsFindTracksEvents.h" +#include "CbmStsTrackFinder.h" +#include "CbmTaskBuildRawEvents.h" +#include "CbmTofSimpClusterizer.h" +#include "CbmTrdClusterFinder.h" +#include "CbmTrdHitProducer.h" + +#include <FairFileSource.h> +#include <FairMonitor.h> +#include <FairParAsciiFileIo.h> +#include <FairParRootFileIo.h> +#include <FairRunAna.h> +#include <FairRuntimeDb.h> +#include <FairSystemInfo.h> + +#include <TStopwatch.h> +#endif + +/*void run_reco(TString input = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/inmed/mc.1.root", + Int_t nTimeSlices = -1, Int_t firstTimeSlice = 0, TString output = "", + TString sEvBuildRaw = "", TString setup = "sis100_electron", TString paramFile = "", Bool_t useMC = false) +{*/ // original void run_reco(const string& mcFile = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/inmed/mc.1.root", const string& parFile = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/inmed/param.1.root", const string& digiFile = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/inmed/digi.1.root", const string& recoFile = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/inmed/reco.1.root", - const string& geoSetup = "sis100_electron", int nEvents = 1000) -{ - TTree::SetMaxTreeSize(90000000000); + const string& setup = "sis100_electron", const string& sEvBuildRaw = "", int nEvents = 1000) +{ // old + //TTree::SetMaxTreeSize(90000000000); + //remove(recoFile.c_str()); // was few lines below before; see commented line there + + + // --- Logger settings ---------------------------------------------------- + TString logLevel = "INFO"; + TString logVerbosity = "LOW"; + // ------------------------------------------------------------------------ + + + // ----- Environment -------------------------------------------------- TString myName = "run_reco"; TString srcDir = gSystem->Getenv("VMCWORKDIR"); + // ------------------------------------------------------------------------ + + + // ----- In- and output file names ------------------------------------ + /*if (input.IsNull()) input = "test"; + TString rawFile = input + ".raw.root"; + TString traFile = input + ".tra.root"; + if (output.IsNull()) output = input; + TString outFile = output + ".reco.root"; + TString monFile = output + ".moni_reco.root"; + if (paramFile.IsNull()) paramFile = input; + TString parFile = paramFile + ".par.root"; + std::cout << "Inputfile " << rawFile << std::endl; + std::cout << "Outfile " << outFile << std::endl; + std::cout << "Parfile " << parFile << std::endl;*/ - remove(recoFile.c_str()); - CbmSetup::Instance()->LoadSetup(geoSetup.c_str()); + // ----- Load the geometry setup ------------------------------------- + std::cout << std::endl << "-I- " << myName << ": Loading setup " << setup << std::endl; + CbmSetup* geo = CbmSetup::Instance(); + geo->LoadSetup(setup.c_str()); + // ------------------------------------------------------------------------ + + // ----- Some global switches ----------------------------------------- + Bool_t eventBased = !(sEvBuildRaw == ""); + Bool_t useMC = true; + Bool_t useMvd = geo->IsActive(ECbmModuleId::kMvd); + Bool_t useSts = geo->IsActive(ECbmModuleId::kSts); + Bool_t useRich = geo->IsActive(ECbmModuleId::kRich); + Bool_t useMuch = geo->IsActive(ECbmModuleId::kMuch); + Bool_t useTrd = geo->IsActive(ECbmModuleId::kTrd); + Bool_t useTof = geo->IsActive(ECbmModuleId::kTof); + Bool_t usePsd = geo->IsActive(ECbmModuleId::kPsd); + // ------------------------------------------------------------------------ + + + // ----- Parameter files as input to the runtime database ------------- std::cout << std::endl << "-I- " << myName << ": Defining parameter files " << std::endl; TList* parFileList = new TList(); TString geoTag; @@ -34,52 +131,295 @@ void run_reco(const string& mcFile = "/lustre/nyx/cbm/users/criesen/cbm/data/l // - TOF digitisation parameters if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTof, geoTag)) { - //TObjString* tofFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digi.par"); - //parFileList->Add(tofFile); - //std::cout << "-I- " << myName << ": Using parameter file " << tofFile->GetString() << std::endl; TObjString* tofBdfFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par"); parFileList->Add(tofBdfFile); std::cout << "-I- " << myName << ": Using parameter file " << tofBdfFile->GetString() << std::endl; } + // ----- Timer -------------------------------------------------------- TStopwatch timer; timer.Start(); - gDebug = 0; + // ------------------------------------------------------------------------ + + // ----- FairRunAna --------------------------------------------------- FairRunAna* run = new FairRunAna(); - FairFileSource* inputSource = new FairFileSource(digiFile.c_str()); - inputSource->AddFriend(mcFile.c_str()); + FairFileSource* inputSource = new FairFileSource(digiFile); // original: FairFileSource(rawFile); + if (useMC) { inputSource->AddFriend(mcFile); } // original: AddFriend(traFile); } run->SetSource(inputSource); run->SetOutputFile(recoFile.c_str()); - run->SetGenerateRunInfo(kFALSE); - - FairLogger::GetLogger()->SetLogScreenLevel("INFO"); - FairLogger::GetLogger()->SetLogVerbosityLevel("LOW"); - - CbmMCDataManager* mcManager = new CbmMCDataManager("MCManager", 1); - mcManager->AddFile(mcFile.c_str()); - run->AddTask(mcManager); - - // Reconstruction tasks - TString macroName = srcDir + "/macro/run/modules/reconstruct.C"; - std::cout << "Loading macro " << macroName << std::endl; - gROOT->LoadMacro(macroName); - Bool_t recoSuccess = gROOT->ProcessLine("reconstruct(true,true)"); - if (!recoSuccess) { - std::cerr << "-E-" << myName << ": error in executing " << macroName << std::endl; - return; + run->SetGenerateRunInfo(kTRUE); + //FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monFile); monFile not yet defined here + // ------------------------------------------------------------------------ + + // ----- MCDataManager ----------------------------------- + if (useMC) { + CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 0); // TODO: ("MCDataManager", 0)?? + mcManager->AddFile(mcFile.c_str()); + run->AddTask(mcManager); } - std::cout << "-I-" << myName << ": " << macroName << " excuted successfully" << std::endl; + // ------------------------------------------------------------------------ + + // ----- Logger settings ---------------------------------------------- + FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data()); + FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data()); + // ------------------------------------------------------------------------ + + + // ----- Raw event building from digis -------------------------------- + if (eventBased) { + if (sEvBuildRaw == "Ideal" || sEvBuildRaw == "ideal" || sEvBuildRaw == "IDEAL") { + FairTask* evBuildRaw = new CbmBuildEventsIdeal(); + run->AddTask(evBuildRaw); + std::cout << "-I- " << myName << ": Added task " << evBuildRaw->GetName() << std::endl; + eventBased = kTRUE; + } //? Ideal raw event building + else if (sEvBuildRaw == "Real" || sEvBuildRaw == "real" || sEvBuildRaw == "REAL") { + CbmTaskBuildRawEvents* evBuildRaw = new CbmTaskBuildRawEvents(); + + //Choose between NoOverlap, MergeOverlap, AllowOverlap + evBuildRaw->SetEventOverlapMode(EOverlapModeRaw::AllowOverlap); + + // Remove detectors where digis not found + if (!useRich) evBuildRaw->RemoveDetector(kRawEventBuilderDetRich); + if (!useMuch) evBuildRaw->RemoveDetector(kRawEventBuilderDetMuch); + if (!usePsd) evBuildRaw->RemoveDetector(kRawEventBuilderDetPsd); + if (!useTof) evBuildRaw->RemoveDetector(kRawEventBuilderDetTof); + if (!useTrd) evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd); + if (!useSts) { + std::cerr << "-E- " << myName << ": Sts must be present for raw event " + << "building using ''Real2019'' option. Terminating macro." << std::endl; + return; + } + // Set STS as reference detector + evBuildRaw->SetReferenceDetector(kRawEventBuilderDetSts); + evBuildRaw->SetTsParameters(0.0, 1.e7, 0.0); + + // Use CbmMuchDigi instead of CbmMuchBeamtimeDigi + evBuildRaw->ChangeMuchBeamtimeDigiFlag(kFALSE); + + evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kSts, 1000); + evBuildRaw->SetTriggerMaxNumber(ECbmModuleId::kSts, -1); + evBuildRaw->SetTriggerWindow(ECbmModuleId::kSts, -500, 500); + + run->AddTask(evBuildRaw); + std::cout << "-I- " << myName << ": Added task " << evBuildRaw->GetName() << std::endl; + eventBased = kTRUE; + } //? Real raw event building + else { + std::cerr << "-E- " << myName << ": Unknown option " << sEvBuildRaw + << " for raw event building! Terminating macro execution." << std::endl; + return; + } + } //? event-based reco + // ------------------------------------------------------------------------ + + + // ----------- QA for raw event builder ----------------------------------- + if (eventBased && useMC) { + CbmBuildEventsQA* evBuildQA = new CbmBuildEventsQA(); + run->AddTask(evBuildQA); + } + // ------------------------------------------------------------------------ + + + // ----- Local reconstruction in MVD ---------------------------------- + if (useMvd) { + + CbmMvdClusterfinder* mvdCluster = new CbmMvdClusterfinder("MVD Cluster Finder", 0, 0); + run->AddTask(mvdCluster); + std::cout << "-I- : Added task " << mvdCluster->GetName() << std::endl; + + CbmMvdHitfinder* mvdHit = new CbmMvdHitfinder("MVD Hit Finder", 0, 0); + mvdHit->UseClusterfinder(kTRUE); + run->AddTask(mvdHit); + std::cout << "-I- " << myName << ": Added task " << mvdHit->GetName() << std::endl; + } + // ------------------------------------------------------------------------ + + + // ----- Local reconstruction in STS ---------------------------------- + if (useSts) { + CbmRecoSts* stsReco = new CbmRecoSts(kCbmRecoTimeslice); + if (eventBased) stsReco->SetMode(kCbmRecoEvent); + run->AddTask(stsReco); + std::cout << "-I- " << myName << ": Added task " << stsReco->GetName() << std::endl; + } + // ------------------------------------------------------------------------ + + + // ----- Local reconstruction in RICH --------------------------------- + if (useRich) { + CbmRichHitProducer* richHitProd = new CbmRichHitProducer(); + run->AddTask(richHitProd); + std::cout << "-I- " << myName << ": Added task " << richHitProd->GetName() << std::endl; + } + // ------------------------------------------------------------------------ + + + // ----- Local reconstruction in MUCH --------------------------------- + if (useMuch) { + + // --- Parameter file name + TString geoTag; + geo->GetGeoTag(ECbmModuleId::kMuch, geoTag); + Int_t muchFlag = (geoTag.Contains("mcbm") ? 1 : 0); + TString parFile = gSystem->Getenv("VMCWORKDIR"); + parFile += "/parameters/much/much_" + geoTag(0, 4) + "_digi_sector.root"; + std::cout << "Using parameter file " << parFile << std::endl; + + // --- Hit finder for GEMs + FairTask* muchReco = new CbmMuchFindHitsGem(parFile.Data(), muchFlag); + run->AddTask(muchReco); + std::cout << "-I- " << myName << ": Added task " << muchReco->GetName() << " with parameter file " << parFile + << std::endl; + } + // ------------------------------------------------------------------------ + + + // ----- Local reconstruction in TRD ---------------------------------- + if (useTrd) { - CbmTrdSetTracksPidLike* trdLI = new CbmTrdSetTracksPidLike("TRDLikelihood", "TRDLikelihood"); - trdLI->SetUseMCInfo(kTRUE); - trdLI->SetUseMomDependence(kTRUE); - run->AddTask(trdLI); + Double_t triggerThreshold = 0.5e-6; // SIS100 + CbmTrdClusterFinder* trdCluster = new CbmTrdClusterFinder(); + if (eventBased) trdCluster->SetTimeBased(kFALSE); + else + trdCluster->SetTimeBased(kTRUE); + trdCluster->SetNeighbourEnable(true, false); + trdCluster->SetMinimumChargeTH(triggerThreshold); + trdCluster->SetRowMerger(true); - CbmMatchRecoToMC* matchRecoToMc = new CbmMatchRecoToMC(); - run->AddTask(matchRecoToMc); + // Uncomment if you want to use all available digis. + // In that case clusters hits will not be added to the CbmEvent + // trdCluster->SetUseOnlyEventDigis(kFALSE); - std::cout << std::endl << std::endl << "-I- " << myName << ": Set runtime DB" << std::endl; + run->AddTask(trdCluster); + std::cout << "-I- " << myName << ": Added task " << trdCluster->GetName() << std::endl; + + CbmTrdHitProducer* trdHit = new CbmTrdHitProducer(); + run->AddTask(trdHit); + std::cout << "-I- " << myName << ": Added task " << trdHit->GetName() << std::endl; + } + // ------------------------------------------------------------------------ + + + // ----- Local reconstruction in TOF ---------------------------------- + if (useTof) { + CbmTofSimpClusterizer* tofCluster = new CbmTofSimpClusterizer("TofSimpClusterizer", 0); + tofCluster->SetOutputBranchPersistent("TofHit", kTRUE); + tofCluster->SetOutputBranchPersistent("TofDigiMatch", kTRUE); + run->AddTask(tofCluster); + std::cout << "-I- " << myName << ": Added task " << tofCluster->GetName() << std::endl; + } + // ------------------------------------------------------------------------ + + + // ----- Local reconstruction in PSD ---------------------------------- + if (usePsd) { + CbmPsdHitProducer* psdHit = new CbmPsdHitProducer(); + run->AddTask(psdHit); + std::cout << "-I- " << myName << ": Added task " << psdHit->GetName() << std::endl; + } + + + // ----- Track finding in STS (+ MVD) ---------------------------------- + if (useMvd || useSts) { + CbmKF* kalman = new CbmKF(); + run->AddTask(kalman); + CbmL1* l1 = 0; + if (useMC) { l1 = new CbmL1("L1", 2, 3); } + else { + l1 = new CbmL1("L1", 0); + } + l1->SetDataMode(!eventBased); + + // --- Material budget file names + TString mvdGeoTag; + if (geo->GetGeoTag(ECbmModuleId::kMvd, mvdGeoTag)) { + TString parFile = gSystem->Getenv("VMCWORKDIR"); + parFile += "/parameters/mvd/mvd_matbudget_" + mvdGeoTag + ".root"; + std::cout << "Using material budget file " << parFile << std::endl; + l1->SetMvdMaterialBudgetFileName(parFile.Data()); + } + TString stsGeoTag; + if (geo->GetGeoTag(ECbmModuleId::kSts, stsGeoTag)) { + TString parFile = gSystem->Getenv("VMCWORKDIR"); + parFile += "/parameters/sts/sts_matbudget_" + stsGeoTag + ".root"; + std::cout << "Using material budget file " << parFile << std::endl; + l1->SetStsMaterialBudgetFileName(parFile.Data()); + } + run->AddTask(l1); + std::cout << "-I- " << myName << ": Added task " << l1->GetName() << std::endl; + + CbmStsTrackFinder* stsTrackFinder = new CbmL1StsTrackFinder(); + if (eventBased) { + FairTask* stsFindTracks = new CbmStsFindTracksEvents(stsTrackFinder, useMvd); + run->AddTask(stsFindTracks); + std::cout << "-I- " << myName << ": Added task " << stsFindTracks->GetName() << std::endl; + } + else { + FairTask* stsFindTracks = new CbmStsFindTracks(0, stsTrackFinder, useMvd); + run->AddTask(stsFindTracks); + std::cout << "-I- " << myName << ": Added task " << stsFindTracks->GetName() << std::endl; + } + } + // ------------------------------------------------------------------------ + + + // ==== From here on, the time-based and the event-based reconstruction + // ==== chains differ, since time-based version of primary vertex finding + // ==== and global tracking are not yet available. For time-based + // ==== reconstruction, a track-based event finder is used; no global + // ==== tracks are produced. + + if (eventBased) { + + // ----- Primary vertex finding ------------------------------------- + CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF(); + CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder); + run->AddTask(findVertex); + std::cout << "-I- " << myName << ": Added task " << findVertex->GetName() << std::endl; + // ---------------------------------------------------------------------- + + + // --- Global track finding ----------------------------------------- + CbmLitFindGlobalTracks* finder = new CbmLitFindGlobalTracks(); + finder->SetTrackingType("branch"); + finder->SetMergerType("nearest_hit"); + run->AddTask(finder); + std::cout << "-I- : Added task " << finder->GetName() << std::endl; + // ---------------------------------------------------------------------- + + + // ----- RICH reconstruction ---------------------------------------- + if (useRich) { + CbmRichReconstruction* richReco = new CbmRichReconstruction(); + run->AddTask(richReco); + std::cout << "-I- : Added task " << richReco->GetName() << std::endl; + } + // ---------------------------------------------------------------------- + + } //? event-based reco + + else { + + // -----Â Event building from STS tracks ----------------------------- + run->AddTask(new CbmBuildEventsFromTracksReal()); + // ---------------------------------------------------------------------- + + } //? time-based reco + + + // ----- Match reco to MC ------------------------------------------------ + if (useMC) { + CbmMatchRecoToMC* match1 = new CbmMatchRecoToMC(); + run->AddTask(match1); + } + + + // ----- Parameter database -------------------------------------------- + std::cout << std::endl << std::endl; + std::cout << "-I- " << myName << ": Set runtime DB" << std::endl; FairRuntimeDb* rtdb = run->GetRuntimeDb(); FairParRootFileIo* parIo1 = new FairParRootFileIo(); FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo(); @@ -89,22 +429,60 @@ void run_reco(const string& mcFile = "/lustre/nyx/cbm/users/criesen/cbm/data/l parIo2->open(parFileList, "in"); rtdb->setSecondInput(parIo2); } + // ------------------------------------------------------------------------ - std::cout << std::endl << "-I- " << myName << ": Initialise run" << std::endl; - run->Init(); + // ----- Run initialisation ------------------------------------------- + std::cout << std::endl; + std::cout << "-I- " << myName << ": Initialise run" << std::endl; + run->Init(); rtdb->setOutput(parIo1); rtdb->saveOutput(); rtdb->print(); + // ------------------------------------------------------------------------ + + + // ----- Register light ions (d, t, He3, He4) ------------------------- + std::cout << std::endl; + TString registerLightIonsMacro = gSystem->Getenv("VMCWORKDIR"); + registerLightIonsMacro += "/macro/KF/registerLightIons.C"; + std::cout << "Loading macro " << registerLightIonsMacro << std::endl; + gROOT->LoadMacro(registerLightIonsMacro); + gROOT->ProcessLine("registerLightIons()"); + // ------------------------------------------------------------------------ + // ----- Start run ---------------------------------------------------- + std::cout << std::endl << std::endl; std::cout << "-I- " << myName << ": Starting run" << std::endl; - run->Run(0, nEvents); + run->Run(0, nEvents); // original: run->Run(firstTimeSlice, nTimeSlices); + // ------------------------------------------------------------------------ + + // ----- Finish ------------------------------------------------------- timer.Stop(); + FairMonitor::GetMonitor()->Print(); + Double_t rtime = timer.RealTime(); + Double_t ctime = timer.CpuTime(); std::cout << std::endl << std::endl; - std::cout << "Macro finished succesfully." << std::endl; - std::cout << "Reco file is " << recoFile << std::endl; + std::cout << "Macro finished successfully." << std::endl; + std::cout << "Output file is " << recoFile << std::endl; std::cout << "Parameter file is " << parFile << std::endl; - std::cout << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s" << std::endl; - std::cout << std::endl << "Test passed" << std::endl << "All ok" << std::endl; -} + std::cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << std::endl; + FairSystemInfo sysInfo; + Float_t maxMemory = sysInfo.GetMaxMemory(); + std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">"; + std::cout << maxMemory; + std::cout << "</DartMeasurement>" << std::endl; + Float_t cpuUsage = ctime / rtime; + std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">"; + std::cout << cpuUsage; + std::cout << "</DartMeasurement>" << std::endl; + // ------------------------------------------------------------------------ + + + // ----- This is to prevent a malloc error when exiting ROOT ---------- + // The source of the error is unknown. Related to TGeoManager. + //RemoveGeoManager(); // TODO: comment this out? + // ------------------------------------------------------------------------ + +} // End of main macro function -- GitLab