Couple STS reconstruction to digi data in CbmDigiEvent format

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";

macro/run/run_reco_digievent.C

+/* 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

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);

reco/tasks/CMakeLists.txt

@@ -10,6 +10,7 @@ Set(LIBRARY_NAME CbmRecoTasks)
 # -----  Compilation sources   ----------------------------
 set(SRCS
 CbmTaskBuildEvents.cxx
+CbmTaskMakeRecoEvents.cxx
 CbmTaskTriggerDigi.cxx
 )
 # ---------------------------------------------------------

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__ */

reco/tasks/CbmTaskMakeRecoEvents.cxx

+/* 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)

reco/tasks/CbmTaskMakeRecoEvents.h

+/* 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 */