From 5ee591f92be9f6af540ffe4e2de8fb53736e2cf7 Mon Sep 17 00:00:00 2001 From: Volker Friese <v.friese@gsi.de> Date: Thu, 9 Dec 2021 09:21:06 +0000 Subject: [PATCH] Couple STS reconstruction to digi data in CbmDigiEvent format --- macro/reco/reco_digi.C | 2 +- macro/run/run_reco_digievent.C | 258 +++++++++++++++++++++++++++ reco/detectors/sts/CbmRecoSts.cxx | 2 +- reco/tasks/CMakeLists.txt | 1 + reco/tasks/CbmRecoTasksLinkDef.h | 1 + reco/tasks/CbmTaskMakeRecoEvents.cxx | 150 ++++++++++++++++ reco/tasks/CbmTaskMakeRecoEvents.h | 82 +++++++++ 7 files changed, 494 insertions(+), 2 deletions(-) create mode 100644 macro/run/run_reco_digievent.C create mode 100644 reco/tasks/CbmTaskMakeRecoEvents.cxx create mode 100644 reco/tasks/CbmTaskMakeRecoEvents.h diff --git a/macro/reco/reco_digi.C b/macro/reco/reco_digi.C index 8b77b4967d..849ac2be70 100644 --- a/macro/reco/reco_digi.C +++ b/macro/reco/reco_digi.C @@ -75,7 +75,7 @@ void reco_digi(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice TString rawFile = input + ".raw.root"; TString traFile = input + ".tra.root"; if (output.IsNull()) output = input; - TString outFile = output + ".reco.root"; + TString outFile = output + ".events.root"; TString monFile = output + ".moni_reco.root"; if (paramFile.IsNull()) paramFile = input; TString parFile = paramFile + ".par.root"; diff --git a/macro/run/run_reco_digievent.C b/macro/run/run_reco_digievent.C new file mode 100644 index 0000000000..f1963615b2 --- /dev/null +++ b/macro/run/run_reco_digievent.C @@ -0,0 +1,258 @@ +/* Copyright (C) 2020-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + 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 + + +/** @brief Macro for CBM reconstruction + ** @author Volker Friese <v.friese@gsi.de> + ** @since 14 November 2020 + ** @param input Name of input file (w/o extension .raw.root) + ** @param nTimeSlices Number of time-slices to process + ** @param firstTimeSlice First time-slice (entry) to be processed + ** @param output Name of output file (w/o extension .rec.root) + ** @param sEvBuildRaw Option for raw event building + ** @param setup Name of predefined geometry setup + ** @param paramFile Parameter ROOT file (w/o extension .par.root) + ** @param debugWithMC Option to provide the trackfinder with MC information + ** + ** This macro performs event-by-event reconstruction from from the digis in DigiEvents. + ** It can be used for real data after unpacking, triggering and event building or + ** for simulated data after triggering and event building with macro/reco/reco_digi.C. + ** + ** The file names must be specified without extensions. The convention is + ** that the raw (input) file is [input].raw.root. The output file + ** will be [input].rec.root if not specified by the user. The parameter file + ** has the extension .par.root. It is assumed to be [input].par.root if + ** not specified by the user. + ** + ** If no argument is specified, the input will be set to "test". This allows + ** to execute the macro chain (run_tra_file.C, run_digi.C and run_reco.C) + ** from the ROOT prompt without user intervention. + ** + **/ +void run_reco_digievent(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice = 0, TString output = "", + TString setup = "sis100_electron", TString paramFile = "") +{ + + // ======================================================================== + // Adjust this part according to your requirements + + // --- Logger settings ---------------------------------------------------- + TString logLevel = "INFO"; + TString logVerbosity = "LOW"; + // ------------------------------------------------------------------------ + + // ----- Environment -------------------------------------------------- + TString myName = "run_reco_digievent"; // this macro's name for screen output + TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory + // ------------------------------------------------------------------------ + + + // ----- In- and output file names ------------------------------------ + if (input.IsNull()) input = "test"; + TString rawFile = input + ".events.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; + + // ----- Load the geometry setup ------------------------------------- + std::cout << std::endl; + std::cout << "-I- " << myName << ": Loading setup " << setup << std::endl; + CbmSetup* geo = CbmSetup::Instance(); + geo->LoadSetup(setup); + // ------------------------------------------------------------------------ + + + // ----- Some global switches ----------------------------------------- + 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; + std::cout << "-I- " << myName << ": Defining parameter files " << std::endl; + TList* parFileList = new TList(); + TString geoTag; + + // - TRD digitisation parameters + if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTrd, geoTag)) { + const Char_t* npar[4] = {"asic", "digi", "gas", "gain"}; + TObjString* trdParFile(NULL); + for (Int_t i(0); i < 4; i++) { + trdParFile = new TObjString(srcDir + "/parameters/trd/trd_" + geoTag + "." + npar[i] + ".par"); + parFileList->Add(trdParFile); + std::cout << "-I- " << myName << ": Using parameter file " << trdParFile->GetString() << std::endl; + } + } + + // - TOF digitisation parameters + if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTof, geoTag)) { + 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; + } + // ------------------------------------------------------------------------ + + // In general, the following parts need not be touched + // ======================================================================== + + + // ----- Timer -------------------------------------------------------- + TStopwatch timer; + timer.Start(); + // ------------------------------------------------------------------------ + + + // ----- FairRunAna --------------------------------------------------- + FairRunAna* run = new FairRunAna(); + FairFileSource* inputSource = new FairFileSource(rawFile); + run->SetSource(inputSource); + run->SetOutputFile(outFile); + run->SetGenerateRunInfo(kTRUE); + FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monFile); + // ------------------------------------------------------------------------ + + + // ----- Logger settings ---------------------------------------------- + FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data()); + FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data()); + // ------------------------------------------------------------------------ + + + // ---- Make Reco Events ---------------------------------------------- + // ---- This is required if the input is in DigiEvent format + auto makeEvents = std::make_unique<CbmTaskMakeRecoEvents>(); + LOG(info) << "-I- " << myName << ": Adding task " << makeEvents->GetName(); + run->AddTask(makeEvents.release()); + // ------------------------------------------------------------------------ + + + // ----- Local reconstruction in STS ---------------------------------- + if (useSts) { + auto recoSts = std::make_unique<CbmRecoSts>(kCbmRecoEvent); + std::cout << "-I- " << myName << ": Adding task " << recoSts->GetName(); + run->AddTask(recoSts.release()); + } + // ------------------------------------------------------------------------ + + + // ----- 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(); + parIo1->open(parFile.Data(), "UPDATE"); + rtdb->setFirstInput(parIo1); + if (!parFileList->IsEmpty()) { + parIo2->open(parFileList, "in"); + rtdb->setSecondInput(parIo2); + } + // ------------------------------------------------------------------------ + + + // ----- Run initialisation ------------------------------------------- + std::cout << std::endl; + std::cout << "-I- " << myName << ": Initialise run" << std::endl; + run->Init(); + rtdb->setOutput(parIo1); + rtdb->saveOutput(); + rtdb->print(); + // ------------------------------------------------------------------------ + + + // ----- Start run ---------------------------------------------------- + std::cout << std::endl << std::endl; + std::cout << "-I- " << myName << ": Starting run" << std::endl; + 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 successfully." << std::endl; + std::cout << "Output file is " << outFile << std::endl; + std::cout << "Parameter file is " << parFile << 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(); + // ------------------------------------------------------------------------ + +} // End of main macro function diff --git a/reco/detectors/sts/CbmRecoSts.cxx b/reco/detectors/sts/CbmRecoSts.cxx index b3eecca84b..7ac75b6d28 100644 --- a/reco/detectors/sts/CbmRecoSts.cxx +++ b/reco/detectors/sts/CbmRecoSts.cxx @@ -226,7 +226,7 @@ void CbmRecoSts::Finish() LOG(info) << "Hits / TSlice : " << fixed << setprecision(2) << fNofHits / Double_t(fNofTs); LOG(info) << "Digis per cluster : " << fixed << setprecision(2) << digiCluster; LOG(info) << "Clusters per hit : " << fixed << setprecision(2) << clusterHit; - LOG(info) << "Time per TSlice : " << fixed << setprecision(2) << 1000. * fTimeRun / Double_t(fNofTs) << " s "; + LOG(info) << "Time per TSlice : " << fixed << setprecision(2) << 1000. * fTimeRun / Double_t(fNofTs) << " ms "; fTimeRun /= Double_t(fNofEvents); fTime1Run /= Double_t(fNofEvents); diff --git a/reco/tasks/CMakeLists.txt b/reco/tasks/CMakeLists.txt index e6e7d6dc31..ef3960862a 100644 --- a/reco/tasks/CMakeLists.txt +++ b/reco/tasks/CMakeLists.txt @@ -10,6 +10,7 @@ Set(LIBRARY_NAME CbmRecoTasks) # ----- Compilation sources ---------------------------- set(SRCS CbmTaskBuildEvents.cxx +CbmTaskMakeRecoEvents.cxx CbmTaskTriggerDigi.cxx ) # --------------------------------------------------------- diff --git a/reco/tasks/CbmRecoTasksLinkDef.h b/reco/tasks/CbmRecoTasksLinkDef.h index 45e9ff31f2..0e54302a7c 100644 --- a/reco/tasks/CbmRecoTasksLinkDef.h +++ b/reco/tasks/CbmRecoTasksLinkDef.h @@ -12,6 +12,7 @@ // --- Classes #pragma link C++ class CbmTaskBuildEvents + ; +#pragma link C++ class CbmTaskMakeRecoEvents + ; #pragma link C++ class CbmTaskTriggerDigi + ; #endif /* __CINT__ */ diff --git a/reco/tasks/CbmTaskMakeRecoEvents.cxx b/reco/tasks/CbmTaskMakeRecoEvents.cxx new file mode 100644 index 0000000000..7094ba466f --- /dev/null +++ b/reco/tasks/CbmTaskMakeRecoEvents.cxx @@ -0,0 +1,150 @@ +/* Copyright (C) 2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Volker Friese [committer] */ + + +#include "CbmTaskMakeRecoEvents.h" + +#include "CbmDefs.h" +#include "CbmEvent.h" + +#include <FairLogger.h> + +#include <TClonesArray.h> +#include <TStopwatch.h> + +#include <algorithm> +#include <cassert> +#include <iomanip> +#include <iostream> +#include <vector> + + +using namespace std; + + +// ----- Constructor ----------------------------------------------------- +CbmTaskMakeRecoEvents::CbmTaskMakeRecoEvents() : FairTask("MakeRecoEvents") {} +// --------------------------------------------------------------------------- + + +// ----- Destructor ------------------------------------------------------ +CbmTaskMakeRecoEvents::~CbmTaskMakeRecoEvents() {} +// --------------------------------------------------------------------------- + + +// ----- Execution ------------------------------------------------------- +void CbmTaskMakeRecoEvents::Exec(Option_t*) +{ + + // --- Timer and counters + TStopwatch timer; + timer.Start(); + + // --- No action if no DigiEvent branch is present + if (!fDigiEvents) return; + + // --- Clear output arrays + fStsDigis->clear(); + fRecoEvents->Clear(); + + + // --- Event loop + Int_t eventNr = 0; + for (auto& digiEvent : *fDigiEvents) { + + // --- Copy StsDigis + size_t startIndex = fStsDigis->size(); + fStsDigis->insert(fStsDigis->end(), digiEvent.fData.fSts.fDigis.begin(), digiEvent.fData.fSts.fDigis.end()); + size_t stopIndex = startIndex + digiEvent.fData.fSts.fDigis.size(); + + // --- Create and fill CbmEvent object + CbmEvent* recoEvent = new ((*fRecoEvents)[eventNr]) CbmEvent(eventNr); + for (size_t index = startIndex; index < stopIndex; index++) { + recoEvent->AddData(ECbmDataType::kStsDigi, index); + } + + eventNr++; + } //# events + + // --- Timeslice log + timer.Stop(); + stringstream logOut; + logOut << setw(20) << left << GetName() << " ["; + logOut << fixed << setw(8) << setprecision(1) << right << timer.RealTime() * 1000. << " ms] "; + logOut << "TS " << fNumTs << ", events " << fDigiEvents->size() << ", STS digis " << fStsDigis->size(); + LOG(info) << logOut.str(); + assert(fDigiEvents->size() == static_cast<size_t>(fRecoEvents->GetEntriesFast())); + + // --- Run statistics + fNumTs++; + fTimeTot += timer.RealTime(); + assert(fDigiEvents->size() == static_cast<size_t>(fRecoEvents->GetEntriesFast())); + fNumEvents += fDigiEvents->size(); +} +// ---------------------------------------------------------------------------- + + +// ----- End-of-timeslice action ------------------------------------------ +void CbmTaskMakeRecoEvents::Finish() +{ + + std::cout << std::endl; + LOG(info) << "====================================="; + LOG(info) << GetName() << ": Run summary"; + LOG(info) << "Timeslices : " << fNumTs; + LOG(info) << "Events : " << fNumEvents; + LOG(info) << "Time / TS : " << fixed << setprecision(2) << 1000. * fTimeTot / double(fNumTs) << " ms"; + LOG(info) << "====================================="; +} +// ---------------------------------------------------------------------------- + + +// ----- Initialisation --------------------------------------------------- +InitStatus CbmTaskMakeRecoEvents::Init() +{ + + LOG(info) << "=================================================="; + LOG(info) << GetName() << ": Initialising "; + + // --- Get FairRootManager instance + FairRootManager* frm = FairRootManager::Instance(); + assert(frm); + + // --- Try to get input vector (CbmDigiEvent) + fDigiEvents = frm->InitObjectAs<const std::vector<CbmDigiEvent>*>("DigiEvent"); + + // --- If DigiEvents are present, create StsDigi and CbmEvent branches + if (fDigiEvents) { + LOG(info) << GetName() << ": Found branch DigiEvent"; + if (frm->GetObject("CbmEvent")) { + LOG(error) << GetName() << ": Found branch CbmEvent! Aborting..."; + return kFATAL; + } + fRecoEvents = new TClonesArray("CbmEvent", 1); + frm->Register("CbmEvent", "Reco events", fRecoEvents, IsOutputBranchPersistent("CbmEvent")); + if (frm->GetObject("StsDigi")) { + LOG(error) << GetName() << ": Found branch StsDigi! Aborting..."; + return kFATAL; + } + fStsDigis = new std::vector<CbmStsDigi>; + frm->RegisterAny("StsDigi", fStsDigis, kFALSE); + } + + // --- If no DigiEvent branch is present, there must be a CbmEvent branch + else { + fRecoEvents = dynamic_cast<TClonesArray*>(frm->GetObject("CbmEvent")); + if (fRecoEvents == nullptr) { + LOG(error) << GetName() << ": Neither DigiEvent nor CbmEvent branch found! Aborting..."; + return kFATAL; + } + } + + LOG(info) << "=================================================="; + std::cout << std::endl; + + return kSUCCESS; +} +// ---------------------------------------------------------------------------- + +ClassImp(CbmTaskMakeRecoEvents) diff --git a/reco/tasks/CbmTaskMakeRecoEvents.h b/reco/tasks/CbmTaskMakeRecoEvents.h new file mode 100644 index 0000000000..b5edfdb9e0 --- /dev/null +++ b/reco/tasks/CbmTaskMakeRecoEvents.h @@ -0,0 +1,82 @@ +/* Copyright (C) 2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Volker Friese [committer] */ + + +#ifndef CBMTASKMAKERECOEVENTS_H +#define CBMTASKMAKERECOEVENTS_H 1 + + +#include "CbmDefs.h" +#include "CbmDigiEvent.h" + +#include <FairTask.h> + +#include <vector> + +class TClonesArray; + + +/** @class CbmTaskMakeRecoEvents + ** @brief Task class for interfacing storable raw events in the CbmDigiEvent format to + ** the current offline reconstruction chain. + ** + ** @author Volker Friese <v.friese@gsi.de> + ** @since 08.12.2021 + ** + ** This tasks creates the established data interfaces (digi branches, CbmEvent) as input + ** to the reconstruction tasks from digis stored in the new event data format CbmDigiEvent + ** as created by trigger and event builder from experiment or simulated data. It is to be + ** understood as intermediate solution until the reconstruction routines will be properly + ** interfaced to the new format. The expense is a duplication of digis in memory. + ** + ** The task has to be run prior to any reconstruction task making use of digis (cluster and + ** hit finder). + ** + ** TOFO: The current implementation is for STS only and shall be expanded to all digi types. + **/ +class CbmTaskMakeRecoEvents : public FairTask { + + +public: + /** @brief Constructor **/ + CbmTaskMakeRecoEvents(); + + + /** @brief Copy constructor (disabled) **/ + CbmTaskMakeRecoEvents(const CbmTaskMakeRecoEvents&) = delete; + + + /** @brief Destructor **/ + virtual ~CbmTaskMakeRecoEvents(); + + + /** @brief Task execution **/ + virtual void Exec(Option_t* opt); + + + /** @brief Finish timeslice **/ + virtual void Finish(); + + + /** @brief Assignment operator (disabled) **/ + CbmTaskMakeRecoEvents& operator=(const CbmTaskMakeRecoEvents&) = delete; + + +private: // methods + /** @brief Task initialisation **/ + virtual InitStatus Init(); + + +private: // members + const std::vector<CbmDigiEvent>* fDigiEvents = nullptr; + TClonesArray* fRecoEvents = nullptr; + std::vector<CbmStsDigi>* fStsDigis = nullptr; + double fTimeTot = 0.; ///< Execution time + size_t fNumTs = 0; ///< Number of processed timeslices + size_t fNumEvents = 0; ///< Number of events + + ClassDef(CbmTaskMakeRecoEvents, 1); +}; + +#endif /* CBMTASKMAKERECOEVENTS_H */ -- GitLab