diff --git a/core/data/CMakeLists.txt b/core/data/CMakeLists.txt
index caf693ea6cb476b825219a39f23841b47bdda259..6121c9f1f474d641dd85ebb593e00e0399d6044a 100644
--- a/core/data/CMakeLists.txt
+++ b/core/data/CMakeLists.txt
@@ -42,6 +42,7 @@ set(SRCS
CbmCluster.cxx
CbmModuleList.cxx
CbmErrorMessage.cxx
+ CbmRawEvent.cxx
CbmMCTrack.cxx
CbmMCEventInfo.cxx
diff --git a/core/data/CbmRawEvent.cxx b/core/data/CbmRawEvent.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..93a6371a2f626c7bc2ee620b8e09c7c70a0ab90c
--- /dev/null
+++ b/core/data/CbmRawEvent.cxx
@@ -0,0 +1,51 @@
+/*
+ * CbmRawEvent.cxx
+ *
+ * Created on: 27.03.2021
+ * Author: vfriese
+ */
+
+
+#include "CbmRawEvent.h"
+
+
+// ----- Number of digis per system
+size_t CbmRawEvent::GetNofDigis(ECbmModuleId system) const
+{
+ switch (system) {
+ case ECbmModuleId::kMvd: return fDigisMvd.size(); break;
+ case ECbmModuleId::kSts: return fDigisSts.size(); break;
+ case ECbmModuleId::kRich: return fDigisRich.size(); break;
+ case ECbmModuleId::kMuch: return fDigisMuch.size(); break;
+ case ECbmModuleId::kTrd: return fDigisTrd.size(); break;
+ case ECbmModuleId::kTof: return fDigisTof.size(); break;
+ case ECbmModuleId::kPsd: return fDigisPsd.size(); break;
+ default: return 0; break;
+ }
+}
+
+
+// ----- Total number of digis
+size_t CbmRawEvent::GetNofDigis() const
+{
+ size_t result = fDigisMvd.size();
+ result += fDigisSts.size();
+ result += fDigisRich.size();
+ result += fDigisMuch.size();
+ result += fDigisTrd.size();
+ result += fDigisTof.size();
+ result += fDigisPsd.size();
+ return result;
+}
+
+
+// ----- String output
+std::string CbmRawEvent::ToString() const
+{
+ std::stringstream ss;
+ ss << "Event " << fNumber << " at t = " << fTime << " ns. Digis: " << GetNofDigis() << " (MVD " << fDigisMvd.size()
+ << " STS " << fDigisSts.size() << " RICH " << fDigisRich.size() << " MUCH " << fDigisMuch.size() << " TRD "
+ << fDigisTrd.size() << " TOF " << fDigisTof.size() << " PSD " << fDigisPsd.size() << ")";
+ return ss.str();
+}
+// -------------------------------------------------------------------------
diff --git a/core/data/CbmRawEvent.h b/core/data/CbmRawEvent.h
new file mode 100644
index 0000000000000000000000000000000000000000..c902606942179556577b582d2ef111016b976224
--- /dev/null
+++ b/core/data/CbmRawEvent.h
@@ -0,0 +1,192 @@
+/** @file CbmRawEvent.h
+ ** @author Volker Friese <v.friese@gsi.de>
+ ** @date 25.03.2021
+ **/
+
+
+#ifndef CBMRAWEVENT_H
+#define CBMRAWEVENT_H 1
+
+#include "CbmDefs.h" // for ECbmDataType, ECbmModuleId::kStsTrack
+//#include "CbmMatch.h" // for CbmMatch
+#include "CbmMuchDigi.h"
+#include "CbmMvdDigi.h"
+#include "CbmPsdDigi.h"
+#include "CbmRichDigi.h"
+#include "CbmStsDigi.h"
+#include "CbmTofDigi.h"
+#include "CbmTrdDigi.h"
+//#include "CbmVertex.h" // for CbmVertex, found in core/data/global
+
+//#include <Rtypes.h> // for THashConsistencyHolder, ClassDef
+//#include <RtypesCore.h> // for Double_t, UInt_t, Int_t
+//#include <TMatrixFSymfwd.h> // for TMatrixFSym
+//#include <TObject.h> // for TObject
+
+
+//#include <map> // for map, map<>::mapped_type
+//#include <span>
+//#include <string> // for string
+//#include <utility>
+#include <vector> // for vector
+
+
+/** @class CbmRawEvent
+ ** @brief Data class as collections of digis for each detector system
+ ** @author V.Friese <v.friese@gsi.de>
+ ** @version 1.0
+ **
+ **/
+class CbmRawEvent : public TObject {
+
+public:
+ /** @brief Default constructor **/
+ CbmRawEvent() {};
+
+
+ /** @brief Constructor with event number and time
+ ** @param number Event number
+ ** @param time Event trigger time [ns]
+ **/
+ CbmRawEvent(uint64_t number, Double_t time = 0) : fNumber(number), fTime(time) {};
+
+
+ /** @brief Copy constructor **/
+ CbmRawEvent(const CbmRawEvent&) = default;
+
+
+ /** @brief Move constructor **/
+ CbmRawEvent(CbmRawEvent&&) = default;
+
+
+ /** @brief Assignment operator **/
+ CbmRawEvent& operator=(const CbmRawEvent&) = default;
+
+
+ /** @brief Destructor **/
+ virtual ~CbmRawEvent() {};
+
+
+ /** @brief Add a digi to the event
+ ** @param data Digi
+ ** TODO: Somehow, this does not work, so dedicated methods were implemented.
+ **/
+ template<class Digi>
+ void AddDigi(ECbmModuleId system, Digi digi)
+ {
+ switch (system) {
+ case ECbmModuleId::kMvd: fDigisMvd.push_back(digi); break;
+ case ECbmModuleId::kSts: fDigisSts.push_back(digi); break;
+ case ECbmModuleId::kRich: fDigisRich.push_back(digi); break;
+ case ECbmModuleId::kMuch: fDigisMuch.push_back(digi); break;
+ case ECbmModuleId::kTrd: fDigisTrd.push_back(digi); break;
+ case ECbmModuleId::kTof: fDigisTof.push_back(digi); break;
+ case ECbmModuleId::kPsd: fDigisPsd.push_back(digi); break;
+ default: break;
+ }
+ }
+
+
+ void AddDigiMvd(CbmMvdDigi digi) { fDigisMvd.push_back(digi); }
+ void AddDigiSts(CbmStsDigi digi) { fDigisSts.push_back(digi); }
+ void AddDigiRich(CbmRichDigi digi) { fDigisRich.push_back(digi); }
+ void AddDigiMuch(CbmMuchDigi digi) { fDigisMuch.push_back(digi); }
+ void AddDigiTrd(CbmTrdDigi digi) { fDigisTrd.push_back(digi); }
+ void AddDigiTof(CbmTofDigi digi) { fDigisTof.push_back(digi); }
+ void AddDigiPsd(CbmPsdDigi digi) { fDigisPsd.push_back(digi); }
+
+
+ /** @brief Read-only access to a single digi
+ ** @param system System identifier (ECbmModuleId)
+ ** @param index Index of digi in its vector
+ ** @return Reference to digi object
+ **/
+ template<class Digi>
+ const Digi& GetDigi(ECbmModuleId system, size_t index) const
+ {
+ switch (system) {
+ case ECbmModuleId::kMvd: return fDigisMvd.at(index); break;
+ case ECbmModuleId::kSts: return fDigisSts.at(index); break;
+ case ECbmModuleId::kRich: return fDigisRich.at(index); break;
+ case ECbmModuleId::kMuch: return fDigisMuch.at(index); break;
+ case ECbmModuleId::kTrd: return fDigisTrd.at(index); break;
+ case ECbmModuleId::kTof: return fDigisTof.at(index); break;
+ case ECbmModuleId::kPsd: return fDigisPsd.at(index); break;
+ default: break;
+ }
+ }
+
+
+ /** @brief Read-only access to digi vectors
+ ** @param system System identifier (ECbmModuleId)
+ ** @return Reference to digi vector for the specified system
+ **/
+ template<class Digi>
+ const std::vector<Digi>& GetDigis(ECbmModuleId system) const
+ {
+ switch (system) {
+ case ECbmModuleId::kMvd: return (&fDigisMvd); break;
+ case ECbmModuleId::kSts: return (&fDigisSts); break;
+ case ECbmModuleId::kRich: return (&fDigisRich); break;
+ case ECbmModuleId::kMuch: return (&fDigisMuch); break;
+ case ECbmModuleId::kTrd: return (&fDigisTrd); break;
+ case ECbmModuleId::kTof: return (&fDigisTof); break;
+ case ECbmModuleId::kPsd: return (&fDigisPsd); break;
+ default: break;
+ }
+ }
+
+
+ /** @brief Number of digis for a given system
+ ** @param systemId System identifier (ECbmModuleId)
+ ** @return Number of digis for the specified system
+ **/
+ size_t GetNofDigis(ECbmModuleId) const;
+
+
+ /** @brief Total number of digis (all systems)
+ ** @return Number of digis in event
+ **/
+ size_t GetNofDigis() const;
+
+
+ /** @brief Event number
+ ** @return Event number
+ **/
+ uint64_t GetNumber() const { return fNumber; }
+
+
+ /** @brief Event time
+ ** @return Event trigger time [ns]
+ **/
+ Double_t GetTime() const { return fTime; }
+
+
+ /** @brief Set event time
+ ** @param time Event trigger time [ns]
+ **/
+ void SetTime(Double_t time) { fTime = time; }
+
+
+ /** @brief String output **/
+ std::string ToString() const;
+
+
+private:
+ /** Event meta data **/
+ uint64_t fNumber; ///< Event number
+ double fTime; ///< Trigger time for this event [ns]
+
+ std::vector<CbmMvdDigi> fDigisMvd = {};
+ std::vector<CbmStsDigi> fDigisSts = {};
+ std::vector<CbmRichDigi> fDigisRich = {};
+ std::vector<CbmMuchDigi> fDigisMuch = {};
+ std::vector<CbmTrdDigi> fDigisTrd = {};
+ std::vector<CbmTofDigi> fDigisTof = {};
+ std::vector<CbmPsdDigi> fDigisPsd = {};
+
+
+ ClassDef(CbmRawEvent, 1);
+};
+
+#endif /* CBMRAWEVENT_H_ */
diff --git a/core/data/DataLinkDef.h b/core/data/DataLinkDef.h
index 8c8156ecb4c39248df45cba7de7dcb933541639d..be8b44cb3c204983c3ebaf83cadfc6654a4f20cb 100644
--- a/core/data/DataLinkDef.h
+++ b/core/data/DataLinkDef.h
@@ -26,6 +26,7 @@
#pragma link C++ class CbmLink + ;
#pragma link C++ class CbmModuleList;
#pragma link C++ class CbmErrorMessage + ;
+#pragma link C++ class CbmRawEvent + ;
// ---In base
#pragma link C++ class CbmDigiBranchBase + ;
@@ -127,10 +128,9 @@
#pragma link C++ class CbmDigiVector < CbmPsdDsp> + ;
#pragma link C++ class vector < CbmEventStore> + ;
-#pragma read sourceClass="CbmTofDigi" version="[1-2]" targetClass="CbmTofDigi" \
- source="UInt_t fuAddress" target="fuAddress" \
- include="Logger.h" \
- code="{ UInt_t system = (onfile.fuAddress >> 0) & ((1 << 4) - 1); \
+#pragma read sourceClass = "CbmTofDigi" version = "[1-2]" targetClass = "CbmTofDigi" source = \
+ "UInt_t fuAddress" target = "fuAddress" include = "Logger.h" code = \
+ "{ UInt_t system = (onfile.fuAddress >> 0) & ((1 << 4) - 1); \
UInt_t smId = (onfile.fuAddress >> 4) & ((1 << 8) - 1); \
UInt_t smType = (onfile.fuAddress >> 12) & ((1 << 4) - 1); \
UInt_t rpcId = (onfile.fuAddress >> 16) & ((1 << 7) - 1); \
@@ -152,10 +152,8 @@
+ ((rpcType & ((1 << 4) - 1)) << 28); \
}"
-#pragma read sourceClass="CbmTofHit" version="[1-4]" targetClass="CbmTofHit" \
- source="" target="" \
- include="Logger.h" \
- code="{ \
+#pragma read sourceClass = "CbmTofHit" version = "[1-4]" targetClass = "CbmTofHit" source = "" target = "" include = \
+ "Logger.h" code = "{ \
LOG(error); \
LOG(error) << \"You are trying to read an outdated version of CbmTofHit\"; \
LOG(error) << \"where the unique tof address can't be converted\"; \
@@ -164,10 +162,8 @@
LOG(fatal) << \"Stop execution.\"; \
}"
-#pragma read sourceClass="CbmTofPoint" version="[1-3]" targetClass="CbmTofPoint" \
- source="" target="" \
- include="Logger.h" \
- code="{ \
+#pragma read sourceClass = "CbmTofPoint" version = "[1-3]" targetClass = "CbmTofPoint" source = "" target = \
+ "" include = "Logger.h" code = "{ \
LOG(error); \
LOG(error) << \"You are trying to read an outdated version of CbmTofPoint\"; \
LOG(error) << \"where the unique tof address can't be converted\"; \
diff --git a/macro/run/run_reco_mcbm.C b/macro/run/run_reco_mcbm.C
new file mode 100644
index 0000000000000000000000000000000000000000..24d42280f788de947f2e0764a852b2ee1b54ec66
--- /dev/null
+++ b/macro/run/run_reco_mcbm.C
@@ -0,0 +1,311 @@
+/* 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 "CbmEvbuildRawTask.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 useMC Option to provide the trackfinder with MC information
+ **
+ ** This macro performs from the digis in a time-slice. It can be used
+ ** for simulated data (result of run_digi.C) or real data after unpacking.
+ **
+ ** The macro covers both time-based reconstruction and event-based
+ ** reconstruction using raw events build from digis. This can be selected
+ ** by the forth argument. If left empty, no raw event builder will be
+ ** employed and reconstruction will be time-based. The option "Ideal"
+ ** selects the ideal raw event builder, which associates digis to events
+ ** based on the MC truth. The option "Real" selects a real raw event builder
+ ** (latest version, for older versions use "Real2018" or "Real2019").
+ **
+ **
+ ** 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_mcbm(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice = 0, TString output = "",
+ TString sEvBuildRaw = "", TString setup = "sis100_electron", TString paramFile = "",
+ Bool_t useMC = false)
+{
+
+ // ========================================================================
+ // Adjust this part according to your requirements
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "debug";
+ TString logVerbosity = "LOW";
+ // ------------------------------------------------------------------------
+
+ // ----- Environment --------------------------------------------------
+ TString myName = "run_reco"; // 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 + ".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;
+
+ // ----- 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 eventBased = !sEvBuildRaw.IsNull();
+ 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);
+ if (useMC) { inputSource->AddFriend(traFile); }
+ run->SetSource(inputSource);
+ run->SetOutputFile(outFile);
+ run->SetGenerateRunInfo(kTRUE);
+ FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monFile);
+ // ------------------------------------------------------------------------
+
+ // ----- MCDataManager -----------------------------------
+ if (useMC) {
+ CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 0);
+ mcManager->AddFile(traFile);
+ run->AddTask(mcManager);
+ }
+ // ------------------------------------------------------------------------
+
+ // ----- Logger settings ----------------------------------------------
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ // ------------------------------------------------------------------------
+
+
+ // ----- Raw event building from digis --------------------------------
+ if (eventBased) {
+ if (sEvBuildRaw.EqualTo("Ideal", TString::ECaseCompare::kIgnoreCase)) {
+ 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.EqualTo("Real", TString::ECaseCompare::kIgnoreCase)) {
+ CbmEvbuildRawTask* evBuildRaw = new CbmEvbuildRawTask();
+
+ //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(kRawEventBuilderDetTof);
+ evBuildRaw->SetTsParameters(0.0, 1.e7, 0.0);
+
+ // Use CbmMuchDigi instead of CbmMuchBeamtimeDigi
+ evBuildRaw->ChangeMuchBeamtimeDigiFlag(kFALSE);
+
+ //evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kSts, 1000);
+ evBuildRaw->SetTriggerMinNumber(ECbmModuleId::kTof, 5);
+ //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
+ // ------------------------------------------------------------------------
+
+
+ // ----- 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/eventbuilder/CMakeLists.txt b/reco/eventbuilder/CMakeLists.txt
index 17f0aaa1af6f661c70b2b23a2d3d5318f393d875..9473646a6df9dbec5e2f4b65e55152414569cb68 100644
--- a/reco/eventbuilder/CMakeLists.txt
+++ b/reco/eventbuilder/CMakeLists.txt
@@ -67,6 +67,8 @@ link_directories( ${LINK_DIRECTORIES})
set(SRCS
+digis/CbmEvbuildRawAlgo.cxx
+digis/CbmEvbuildRawTask.cxx
digis/CbmAlgoBuildRawEvents.cxx
digis/CbmTaskBuildRawEvents.cxx
digis/CbmBuildEventsIdeal.cxx
diff --git a/reco/eventbuilder/CbmEventBuilderLinkDef.h b/reco/eventbuilder/CbmEventBuilderLinkDef.h
index f6743998e726cae0cb26d36259d71097d861b518..0837321fcb068b1e50c1b634377452b0c42d32e4 100644
--- a/reco/eventbuilder/CbmEventBuilderLinkDef.h
+++ b/reco/eventbuilder/CbmEventBuilderLinkDef.h
@@ -8,6 +8,8 @@
#pragma link off all classes;
#pragma link off all functions;
+#pragma link C++ class CbmEvbuildRawAlgo + ;
+#pragma link C++ class CbmEvbuildRawTask + ;
#pragma link C++ class CbmAlgoBuildRawEvents + ;
#pragma link C++ class CbmTaskBuildRawEvents + ;
#pragma link C++ class CbmBuildEventsFromTracksIdeal + ;
diff --git a/reco/eventbuilder/digis/CbmEvbuildRawAlgo.cxx b/reco/eventbuilder/digis/CbmEvbuildRawAlgo.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0610a9e18bdd0dfe42598da67c20a46d51ed439a
--- /dev/null
+++ b/reco/eventbuilder/digis/CbmEvbuildRawAlgo.cxx
@@ -0,0 +1,990 @@
+/********************************************************************************
+ * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
+ * *
+ * This software is distributed under the terms of the *
+ * GNU Lesser General Public Licence (LGPL) version 3, *
+ * copied verbatim in the file "LICENSE" *
+ ********************************************************************************/
+#include "CbmEvbuildRawAlgo.h"
+
+/// CBM headers
+#include "CbmMuchBeamTimeDigi.h"
+#include "CbmMuchDigi.h"
+#include "CbmPsdDigi.h"
+#include "CbmRawEvent.h"
+#include "CbmRichDigi.h"
+#include "CbmStsDigi.h"
+#include "CbmTofDigi.h"
+#include "CbmTrdDigi.h"
+
+#include "TimesliceMetaData.h"
+
+/// FAIRROOT headers
+#include <FairLogger.h>
+#include <FairRootManager.h>
+#include <FairRunOnline.h>
+
+/// FAIRSOFT headers (geant, boost, ...)
+#include <TCanvas.h>
+#include <TClonesArray.h>
+#include <TFolder.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <THttpServer.h>
+#include <TStopwatch.h>
+
+template<>
+void CbmEvbuildRawAlgo::LoopOnSeeds<Double_t>();
+
+Bool_t CbmEvbuildRawAlgo::InitAlgo()
+{
+ LOG(info) << "CbmEvbuildRawAlgo::InitAlgo => Starting sequence";
+
+ if (fbGetTimings) {
+ fTimer = new TStopwatch;
+ fTimer->Start();
+ }
+
+ /// Check if reference detector is set and seed data are available,
+ /// otherwise look for explicit seed times
+ if (fRefDet.detId == ECbmModuleId::kNotExist) {
+ if (fSeedTimes == nullptr) {
+ LOG(fatal) << "No reference detector set and no seed times supplied, stopping there!";
+ }
+ }
+ else {
+ if (fSeedTimes != nullptr) {
+ LOG(fatal) << "Cannot have explicit seed times and reference detector, stopping there!";
+ }
+ if (kFALSE == CheckDataAvailable(fRefDet)) {
+ LOG(fatal) << "Reference detector set but no digi input found, stopping there!";
+ }
+ }
+
+ /// Check if data for detectors in selection list are available
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ if (kFALSE == CheckDataAvailable(*det)) {
+ LOG(fatal) << "No digi input for one of selection detector, stopping there!";
+ }
+ }
+
+ /// Access the TS metadata to know TS start time if needed
+ if (fbUseTsMetaData) {
+ if (!fTimeSliceMetaDataArray) {
+ LOG(fatal) << "No TS metadata input found"
+ << " => Please check in the unpacking macro if the following line was "
+ "present!"
+ << std::endl
+ << "source->SetWriteOutputFlag(kTRUE); // For writing TS metadata";
+ }
+ }
+ if (fbFillHistos) { CreateHistograms(); }
+ if (fTimer != nullptr) {
+ fTimer->Stop();
+ Double_t rtime = fTimer->RealTime();
+ Double_t ctime = fTimer->CpuTime();
+ LOG(info) << "CbmEvbuildRawAlgo::Init(): Real time " << rtime << " s, CPU time " << ctime << " s";
+ }
+
+ LOG(info) << "CbmEvbuildRawAlgo::InitAlgo => Done";
+ return kTRUE;
+}
+
+void CbmEvbuildRawAlgo::Finish()
+{
+ if (fbGetTimings) { PrintTimings(); }
+}
+
+void CbmEvbuildRawAlgo::PrintTimings()
+{
+ if (fTimer == nullptr) { LOG(fatal) << "Trying to print timings but timer not set"; }
+ else {
+ Double_t rtime = fTimer->RealTime();
+ Double_t ctime = fTimer->CpuTime();
+ LOG(info) << "CbmEvbuildRawAlgo: Real time " << rtime << " s, CPU time " << ctime << " s";
+ }
+}
+
+void CbmEvbuildRawAlgo::ClearEventVector()
+{
+ /// Need to delete the object the pointer points to first
+ int counter = 0;
+ for (CbmRawEvent* event : fEventVector) {
+ LOG(debug) << "Event " << counter << " has " << event->GetNofDigis() << " digis";
+ delete event;
+ counter++;
+ }
+ fEventVector.clear();
+}
+
+void CbmEvbuildRawAlgo::ProcessTs()
+{
+ LOG_IF(info, fuNrTs % 1000 == 0) << "Begin of TS " << fuNrTs;
+ if (fTimer != nullptr) { fTimer->Start(kFALSE); }
+ LOG(debug) << "Init timeslice";
+ InitTs();
+ LOG(debug) << "Init seed window";
+ InitSeedWindow();
+ LOG(debug) << "Build events";
+ BuildEvents();
+
+ /// Store last event with trigger if not done by other seed
+ if (nullptr != fCurrentEvent) {
+ /// TODO: store start time of current event ?
+ // fCurrentEvent->SetStartTime( fPrevTime ); // Replace Seed time with time of first digi in event?
+ //fCurrentEvent->SetEndTime(fdPrevEvtEndTime);
+ fEventVector.push_back(fCurrentEvent);
+ fuCurEv++;
+
+ /// Prevent building over TS edge
+ fCurrentEvent = nullptr;
+ }
+
+ LOG(debug) << "Found " << fEventVector.size() << " triggered events";
+ if (fbFillHistos) { FillHistos(); }
+ if (fTimer != nullptr) { fTimer->Stop(); }
+
+ fuNrTs++;
+}
+
+void CbmEvbuildRawAlgo::InitTs()
+{
+ /// Reset TS based variables (analysis per TS = no building over the border)
+
+ /// Reference detector
+ if (fRefDet.detId != ECbmModuleId::kNotExist) {
+ fRefDet.fuStartIndex = 0;
+ fRefDet.fuEndIndex = 0;
+ }
+ /// Loop on detectors in selection list
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ (*det).fuStartIndex = 0;
+ (*det).fuEndIndex = 0;
+ }
+}
+
+void CbmEvbuildRawAlgo::InitSeedWindow()
+{
+ /// Access the TS metadata if needed to know TS start time and overlap size
+ Double_t dTsStartTime = fdTsStartTime;
+ Double_t dOverlapStart = fdTsStartTime + fdTsLength;
+ Double_t dOverlapSize = fdTsOverLength;
+
+ if (fbUseTsMetaData) {
+ const TimesliceMetaData* pTsMetaData = dynamic_cast<TimesliceMetaData*>(fTimeSliceMetaDataArray->At(0));
+ if (nullptr == pTsMetaData)
+ LOG(fatal) << Form("CbmEvbuildRawAlgo::LoopOnSeeds => "
+ "No TS metadata found for TS %6u.",
+ fuNrTs);
+ dTsStartTime = pTsMetaData->GetStartTime();
+ dOverlapStart = pTsMetaData->GetOverlapStartTime();
+ dOverlapSize = pTsMetaData->GetOverlapDuration();
+ }
+
+ /// Print warning in first TS if time window borders out of potential overlap
+ if ((0.0 < fdEarliestTimeWinBeg && dOverlapSize < fdLatestTimeWinEnd) || (dOverlapSize < fdWidestTimeWinRange)) {
+ LOG(warning) << "CbmEvbuildRawAlgo::LoopOnSeeds => "
+ << Form("Event window not fitting in TS overlap, risk of "
+ "incomplete events: %f %f %f %f",
+ fdEarliestTimeWinBeg, fdLatestTimeWinEnd, fdWidestTimeWinRange, dOverlapSize);
+ } // if end of event window does not fit in overlap for a seed at edge of TS core
+
+ /// Define an acceptance window for the seeds in order to use the overlap
+ /// part of the TS to avoid incomplete events
+ if (fbIgnoreTsOverlap) {
+ fdSeedWindowBeg = dTsStartTime;
+ fdSeedWindowEnd = dOverlapStart;
+ }
+ else {
+ fdSeedWindowBeg = dTsStartTime + (0.0 < fdEarliestTimeWinBeg ? 0.0 : -fdEarliestTimeWinBeg);
+ fdSeedWindowEnd = dOverlapStart + (0.0 < fdEarliestTimeWinBeg ? 0.0 : -fdEarliestTimeWinBeg);
+ }
+}
+
+void CbmEvbuildRawAlgo::BuildEvents()
+{
+ /// Call LoopOnSeed with proper template argument
+ switch (fRefDet.detId) {
+ case ECbmModuleId::kSts: {
+ LoopOnSeeds<CbmStsDigi>();
+ break;
+ }
+ case ECbmModuleId::kMuch: {
+ if (fbUseMuchBeamtimeDigi) { LoopOnSeeds<CbmMuchBeamTimeDigi>(); }
+ else {
+ LoopOnSeeds<CbmMuchDigi>();
+ }
+ break;
+ }
+ case ECbmModuleId::kTrd: {
+ LoopOnSeeds<CbmTrdDigi>();
+ break;
+ }
+ case ECbmModuleId::kTof: {
+ LoopOnSeeds<CbmTofDigi>();
+ break;
+ }
+ case ECbmModuleId::kRich: {
+ LoopOnSeeds<CbmRichDigi>();
+ break;
+ }
+ case ECbmModuleId::kPsd: {
+ LoopOnSeeds<CbmPsdDigi>();
+ break;
+ }
+ case ECbmModuleId::kT0: {
+ LoopOnSeeds<CbmTofDigi>();
+ break;
+ }
+ case ECbmModuleId::kNotExist: { //explicit seed times
+ LoopOnSeeds<Double_t>();
+ break;
+ }
+ default: {
+ LOG(fatal) << "CbmEvbuildRawAlgo::BuildEvents => "
+ << "Trying to search event seeds with unsupported det: " << fRefDet.sName;
+ break;
+ }
+ }
+}
+
+template<>
+void CbmEvbuildRawAlgo::LoopOnSeeds<Double_t>()
+{
+ if (ECbmModuleId::kNotExist == fRefDet.detId) {
+ const UInt_t uNbSeeds = fSeedTimes->size();
+ /// Loop on size of vector
+ for (UInt_t uSeed = 0; uSeed < uNbSeeds; ++uSeed) {
+ LOG(debug) << Form("Checking seed %6u / %6u", uSeed, uNbSeeds);
+ Double_t dTime = fSeedTimes->at(uSeed);
+
+ /// Check if seed in acceptance window (is this needed here?)
+ if (dTime < fdSeedWindowBeg) { continue; }
+ else if (fdSeedWindowEnd < dTime) {
+ break;
+ }
+ /// Check Seed and build event if needed
+ CheckSeed(dTime, uSeed);
+ }
+ }
+ else {
+ LOG(fatal) << "Trying to read explicit seeds while reference detector is set.";
+ }
+}
+
+template<class DigiSeed>
+void CbmEvbuildRawAlgo::LoopOnSeeds()
+{
+ if (ECbmModuleId::kT0 == fRefDet.detId) {
+ const UInt_t uNbRefDigis = fT0DigiVec->size();
+ /// Loop on size of vector
+ for (UInt_t uDigi = 0; uDigi < uNbRefDigis; ++uDigi) {
+ LOG(debug) << Form("Checking seed %6u / %6u", uDigi, uNbRefDigis);
+ Double_t dTime = fT0DigiVec->at(uDigi).GetTime();
+ /// Check Seed and build event if needed
+ CheckSeed(dTime, uDigi);
+ }
+ }
+ else {
+ const UInt_t uNbRefDigis = GetNofDigis(fRefDet.detId);
+ /// Loop on size of vector
+ for (UInt_t uDigi = 0; uDigi < uNbRefDigis; ++uDigi) {
+ LOG(debug) << Form("Checking seed %6u / %6u", uDigi, uNbRefDigis);
+ const DigiSeed* pDigi = GetDigi<DigiSeed>(uDigi);
+ const Double_t dTime = pDigi->GetTime();
+
+ /// Check if seed in acceptance window
+ if (dTime < fdSeedWindowBeg) { continue; }
+ else if (fdSeedWindowEnd < dTime) {
+ break;
+ }
+ /// Check Seed and build event if needed
+ CheckSeed(dTime, uDigi);
+ }
+ }
+}
+
+Double_t CbmEvbuildRawAlgo::GetSeedTimeWinRange()
+{
+ if (ECbmModuleId::kNotExist != fRefDet.detId) { return fRefDet.GetTimeWinRange(); }
+ else {
+ return fdSeedTimeWinEnd - fdSeedTimeWinBeg;
+ }
+}
+
+void CbmEvbuildRawAlgo::CheckSeed(Double_t dSeedTime, UInt_t uSeedDigiIdx)
+{
+ /// If previous event valid and event overlap not allowed, check if we are in overlap
+ /// and react accordingly
+ if (nullptr != fCurrentEvent
+ && (EOverlapModeRaw::AllowOverlap != fOverMode || dSeedTime - fdPrevEvtTime < GetSeedTimeWinRange())
+ && dSeedTime - fdPrevEvtTime < fdWidestTimeWinRange) {
+ /// Within overlap range
+ switch (fOverMode) {
+ case EOverlapModeRaw::NoOverlap: {
+ /// No overlap allowed => reject
+ LOG(debug1) << "Reject seed due to overlap";
+ return;
+ break;
+ }
+ case EOverlapModeRaw::MergeOverlap: {
+ /// Merge overlap mode => do nothing and go on filling current event
+ break;
+ }
+ case EOverlapModeRaw::AllowOverlap: {
+ /// In allow overlap mode => reject only if reference det is in overlap
+ /// to avoid cloning events due to single seed cluster
+ LOG(debug1) << "Reject seed because part of cluster of previous one";
+ return;
+ break;
+ }
+ }
+ } // if( prev Event exists and mode forbiden overlap present )
+ else {
+ /// Out of overlap range or in overlap allowed mode
+ /// => store previous event if not empty and create new one
+ if (nullptr != fCurrentEvent) {
+ /// TODO: store start time of current event ?
+ // fCurrentEvent->SetStartTime( fPrevTime ); // Replace Seed time with time of first digi in event?
+ //fCurrentEvent->SetEndTime(fdPrevEvtEndTime);
+ fEventVector.push_back(fCurrentEvent);
+ fuCurEv++;
+ }
+ fCurrentEvent = new CbmRawEvent(fuCurEv, dSeedTime);
+ } // else of if( prev Event exists and mode forbiden overlap present )
+
+ if (fRefDet.detId != ECbmModuleId::kNotExist) {
+ /// If window open for reference detector, search for other reference Digis matching it
+ /// Otherwise only add the current seed
+ if (fRefDet.fdTimeWinBeg < fRefDet.fdTimeWinEnd) {
+ SearchMatches(dSeedTime, fRefDet);
+ /// Also add the seed if the window starts after the seed
+ if (0 < fRefDet.fdTimeWinBeg) AddDigiToEvent(fRefDet, uSeedDigiIdx);
+ }
+ else {
+ AddDigiToEvent(fRefDet, uSeedDigiIdx);
+ }
+ }
+
+ /// Search for matches for each detector in selection list
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ SearchMatches(dSeedTime, *det);
+ }
+
+ CheckTriggerCondition(dSeedTime);
+}
+
+void CbmEvbuildRawAlgo::SearchMatches(Double_t dSeedTime, RawEventBuilderDetector& detMatch)
+{
+ switch (detMatch.detId) {
+ case ECbmModuleId::kSts: {
+ SearchMatches<CbmStsDigi>(dSeedTime, detMatch);
+ break;
+ }
+ case ECbmModuleId::kMuch: {
+ if (fbUseMuchBeamtimeDigi) { SearchMatches<CbmMuchBeamTimeDigi>(dSeedTime, detMatch); }
+ else {
+ SearchMatches<CbmMuchDigi>(dSeedTime, detMatch);
+ }
+ break;
+ }
+ case ECbmModuleId::kTrd: {
+ SearchMatches<CbmTrdDigi>(dSeedTime, detMatch);
+ break;
+ }
+ case ECbmModuleId::kTof: {
+ SearchMatches<CbmTofDigi>(dSeedTime, detMatch);
+ break;
+ }
+ case ECbmModuleId::kRich: {
+ SearchMatches<CbmRichDigi>(dSeedTime, detMatch);
+ break;
+ }
+ case ECbmModuleId::kPsd: {
+ SearchMatches<CbmPsdDigi>(dSeedTime, detMatch);
+ break;
+ }
+ case ECbmModuleId::kT0: {
+ SearchMatches<CbmTofDigi>(dSeedTime, detMatch);
+ break;
+ }
+ default: {
+ LOG(fatal) << "CbmEvbuildRawAlgo::LoopOnSeeds => "
+ << "Trying to search matches with unsupported det: " << detMatch.sName << std::endl
+ << "You may want to add support for it in the method.";
+ break;
+ }
+ }
+}
+
+template<class DigiCheck>
+void CbmEvbuildRawAlgo::SearchMatches(Double_t dSeedTime, RawEventBuilderDetector& detMatch)
+{
+ /// This algo relies on time sorted vectors for the selected detectors
+ UInt_t uLocalIndexStart = detMatch.fuStartIndex;
+ UInt_t uLocalIndexEnd = detMatch.fuStartIndex;
+
+ /// Check the Digis until out of window
+ if (ECbmModuleId::kT0 == detMatch.detId) {
+
+ /// Loop on size of vector
+ const UInt_t uNbSelDigis = fT0DigiVec->size();
+ for (UInt_t uDigi = detMatch.fuStartIndex; uDigi < uNbSelDigis; ++uDigi) {
+ const Double_t dTime = fT0DigiVec->at(uDigi).GetTime();
+ const Double_t dTimeDiff = dTime - dSeedTime;
+
+ /// Check if within time window, update start/stop indices if needed
+ if (dTimeDiff < detMatch.fdTimeWinBeg) {
+ ++uLocalIndexStart;
+ continue;
+ }
+ else if (detMatch.fdTimeWinEnd < dTimeDiff) {
+ /// Store as end the first digi out of window to avoid double counting in case of
+ /// merged overlap event mode
+ uLocalIndexEnd = uDigi;
+ break;
+ }
+ AddDigiToEvent(detMatch, uDigi);
+ if (fdPrevEvtEndTime < dTime) fdPrevEvtEndTime = dTime;
+ }
+
+ /// catch the case where we reach the end of the vector before being out of the time window
+ if (uLocalIndexEnd < uLocalIndexStart) uLocalIndexEnd = uNbSelDigis;
+ }
+ else {
+ const UInt_t uNbSelDigis = GetNofDigis(detMatch.detId);
+ /// Loop on size of vector
+ for (UInt_t uDigi = detMatch.fuStartIndex; uDigi < uNbSelDigis; ++uDigi) {
+ const DigiCheck* pDigi = GetDigi<DigiCheck>(uDigi);
+ const Double_t dTime = pDigi->GetTime();
+ const Double_t dTimeDiff = dTime - dSeedTime;
+ LOG(debug4) << detMatch.sName << Form(" => Checking match %6u / %6u, dt %f", uDigi, uNbSelDigis, dTimeDiff);
+
+ /// Check if within time window, update start/stop indices if needed
+ if (dTimeDiff < detMatch.fdTimeWinBeg) {
+ ++uLocalIndexStart;
+ continue;
+ }
+ else if (detMatch.fdTimeWinEnd < dTimeDiff) {
+ /// Store as end the first digi out of window to avoid double counting in case of
+ /// merged overlap event mode
+ uLocalIndexEnd = uDigi;
+ break;
+ }
+ AddDigiToEvent(detMatch, uDigi);
+ if (fdPrevEvtEndTime < dTime) fdPrevEvtEndTime = dTime;
+ }
+ /// catch the case where we reach the end of the vector before being out of the time window
+ if (uLocalIndexEnd < uLocalIndexStart) uLocalIndexEnd = uNbSelDigis;
+ }
+
+ /// Update the StartIndex and EndIndex for the next event seed
+ detMatch.fuStartIndex = uLocalIndexStart;
+ detMatch.fuEndIndex = uLocalIndexEnd;
+}
+
+void CbmEvbuildRawAlgo::AddDigiToEvent(RawEventBuilderDetector& det, Int_t _entry)
+{
+ switch (det.detId) {
+ case ECbmModuleId::kSts: fCurrentEvent->AddDigiSts(fStsDigis->at(_entry)); break;
+ case ECbmModuleId::kRich: fCurrentEvent->AddDigiRich(fRichDigis->at(_entry)); break;
+ case ECbmModuleId::kMuch: fCurrentEvent->AddDigiMuch(fMuchDigis->at(_entry)); break;
+ case ECbmModuleId::kTrd: fCurrentEvent->AddDigiTrd(fTrdDigis->at(_entry)); break;
+ case ECbmModuleId::kTof: fCurrentEvent->AddDigiTof(fTofDigis->at(_entry)); break;
+ case ECbmModuleId::kPsd: fCurrentEvent->AddDigiPsd(fPsdDigis->at(_entry)); break;
+ default: break;
+ }
+}
+
+void CbmEvbuildRawAlgo::CheckTriggerCondition(Double_t dSeedTime)
+{
+ /// Check if event is filling trigger conditions and clear it if not
+ if (HasTrigger(fCurrentEvent)) {
+ fdPrevEvtTime = dSeedTime;
+
+ /// In case of NoOverlap or MergeOverlap, we can and should start checking the next window
+ /// from end of current window in order to save CPU and avoid duplicating
+ if (EOverlapModeRaw::NoOverlap == fOverMode || EOverlapModeRaw::MergeOverlap == fOverMode) {
+ /// Update reference detector
+ if (fRefDet.detId != ECbmModuleId::kNotExist) { fRefDet.fuStartIndex = fRefDet.fuEndIndex; }
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ (*det).fuStartIndex = (*det).fuEndIndex;
+ }
+ }
+ }
+ else {
+ LOG(debug1) << "Reject seed due to Trigger requirements";
+ delete fCurrentEvent;
+ fCurrentEvent = nullptr; /// delete does NOT set a pointer to nullptr...
+ }
+}
+
+Bool_t CbmEvbuildRawAlgo::HasTrigger(CbmRawEvent* event)
+{
+ /// Check first reference detector
+ if (fRefDet.detId != ECbmModuleId::kNotExist) {
+ if (kFALSE == CheckTriggerConditions(event, fRefDet)) { return kFALSE; }
+ }
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ if (kFALSE == CheckTriggerConditions(event, *det)) return kFALSE;
+ }
+ /// All Ok, trigger is there
+ return kTRUE;
+}
+
+Bool_t CbmEvbuildRawAlgo::CheckTriggerConditions(CbmRawEvent* event, RawEventBuilderDetector& det)
+{
+ /// Check if both Trigger conditions disabled for this detector
+ if (0 == det.fuTriggerMinDigis && det.fiTriggerMaxDigis < 0) { return kTRUE; }
+
+ /// Check if detector present
+ if (!CheckDataAvailable(det.detId)) {
+ LOG(warning) << "Event does not have digis storage for " << det.sName
+ << " while the following trigger min/max are defined: " << det.fuTriggerMinDigis << " "
+ << det.fiTriggerMaxDigis;
+ return kFALSE;
+ }
+
+ /// Check trigger rejection by minimal number or absence
+ const Int_t iNbDigis = event->GetNofDigis(det.detId);
+ if ((-1 == iNbDigis) || (static_cast<UInt_t>(iNbDigis) < det.fuTriggerMinDigis)) {
+ LOG(debug2) << "Event does not have enough digis: " << iNbDigis << " vs " << det.fuTriggerMinDigis << " for "
+ << det.sName;
+ return kFALSE;
+ }
+
+ /// Check trigger rejection by maximal number
+ else if (0 < det.fiTriggerMaxDigis && det.fiTriggerMaxDigis < iNbDigis) {
+ LOG(debug2) << "Event Has too many digis: " << iNbDigis << " vs " << det.fiTriggerMaxDigis << " for " << det.sName;
+ return kFALSE;
+ }
+ else {
+ return kTRUE;
+ }
+}
+
+//----------------------------------------------------------------------
+
+Bool_t CbmEvbuildRawAlgo::CheckDataAvailable(RawEventBuilderDetector& det)
+{
+ if (!CheckDataAvailable(det.detId)) {
+ LOG(info) << "No " << det.sName << " digi input found.";
+ return kFALSE;
+ }
+ return kTRUE;
+}
+
+bool CbmEvbuildRawAlgo::CheckDataAvailable(ECbmModuleId detId)
+{
+ switch (detId) {
+ case ECbmModuleId::kSts: {
+ return fStsDigis != nullptr;
+ }
+ case ECbmModuleId::kMuch: {
+ if (fbUseMuchBeamtimeDigi) { return fMuchBeamTimeDigis != nullptr; }
+ else {
+ return fMuchDigis != nullptr;
+ }
+ }
+ case ECbmModuleId::kTrd: {
+ return fTrdDigis != nullptr;
+ }
+ case ECbmModuleId::kTof: {
+ return fTofDigis != nullptr;
+ }
+ case ECbmModuleId::kRich: {
+ return fRichDigis != nullptr;
+ }
+ case ECbmModuleId::kPsd: {
+ return fPsdDigis != nullptr;
+ }
+ case ECbmModuleId::kT0: {
+ return fT0DigiVec != nullptr;
+ }
+ default: {
+ LOG(fatal) << "CbmEvbuildRawAlgo::CheckDataAvailable => "
+ << "Unsupported detector.";
+ return -1;
+ }
+ }
+}
+
+UInt_t CbmEvbuildRawAlgo::GetNofDigis(ECbmModuleId detId)
+{
+ switch (detId) {
+ case ECbmModuleId::kSts: {
+ return fStsDigis->size();
+ }
+ case ECbmModuleId::kMuch: {
+ if (fbUseMuchBeamtimeDigi) { return fMuchBeamTimeDigis->size(); }
+ else {
+ return fMuchDigis->size();
+ }
+ }
+ case ECbmModuleId::kTrd: {
+ return fTrdDigis->size();
+ }
+ case ECbmModuleId::kTof: {
+ return fTofDigis->size();
+ }
+ case ECbmModuleId::kRich: {
+ return fRichDigis->size();
+ }
+ case ECbmModuleId::kPsd: {
+ return fPsdDigis->size();
+ }
+ case ECbmModuleId::kT0: {
+ return fT0DigiVec->size(); //what to do here? Not in digi manager.
+ }
+ default: {
+ LOG(fatal) << "CbmEvbuildRawAlgo::GetNofDigis => "
+ << "Trying to get digi number with unsupported detector.";
+ return -1;
+ }
+ }
+}
+
+template<>
+const CbmStsDigi* CbmEvbuildRawAlgo::GetDigi(UInt_t uDigi)
+{
+ return &((*fStsDigis)[uDigi]);
+}
+template<>
+const CbmMuchBeamTimeDigi* CbmEvbuildRawAlgo::GetDigi(UInt_t uDigi)
+{
+ return &((*fMuchBeamTimeDigis)[uDigi]);
+}
+template<>
+const CbmMuchDigi* CbmEvbuildRawAlgo::GetDigi(UInt_t uDigi)
+{
+ return &((*fMuchDigis)[uDigi]);
+}
+template<>
+const CbmTrdDigi* CbmEvbuildRawAlgo::GetDigi(UInt_t uDigi)
+{
+ return &((*fTrdDigis)[uDigi]);
+}
+template<>
+const CbmTofDigi* CbmEvbuildRawAlgo::GetDigi(UInt_t uDigi)
+{
+ return &((*fTofDigis)[uDigi]);
+}
+template<>
+const CbmRichDigi* CbmEvbuildRawAlgo::GetDigi(UInt_t uDigi)
+{
+ return &((*fRichDigis)[uDigi]);
+}
+template<>
+const CbmPsdDigi* CbmEvbuildRawAlgo::GetDigi(UInt_t uDigi)
+{
+ return &((*fPsdDigis)[uDigi]);
+}
+
+//----------------------------------------------------------------------
+void CbmEvbuildRawAlgo::CreateHistograms()
+{
+ outFolder = new TFolder("AlgoBuildRawEventsHist", " AlgoBuildRawEvents Histos");
+ outFolder->Clear();
+
+ fhEventTime = new TH1F("hEventTime", "seed time of the events; Seed time [s]; Events", 1000, 0, 0.001);
+ fhEventTime->SetCanExtend(TH1::kAllAxes);
+
+ fhEventDt =
+ new TH1F("fhEventDt", "interval in seed time of consecutive events; Seed time [s]; Events", 1000, 0, 0.0001);
+ fhEventDt->SetCanExtend(TH1::kAllAxes);
+
+ fhEventSize = new TH1F("hEventSize", "nb of all digis in the event; Nb Digis []; Events []", 1000, 0, 100);
+ fhEventSize->SetCanExtend(TH1::kAllAxes);
+
+ fhNbDigiPerEvtTime = new TH2I("hNbDigiPerEvtTime",
+ "nb of all digis per event vs seed time of the events; Seed time "
+ "[s]; Nb Digis []; Events []",
+ 1000, 0, 0.001, 1000, 0, 100);
+ fhNbDigiPerEvtTime->SetCanExtend(TH2::kAllAxes);
+
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ /// In case name not provided, do not create the histo to avoid name conflicts!
+ if ("Invalid" == (*det).sName) {
+ fvhNbDigiPerEvtTimeDet.push_back(nullptr);
+ continue;
+ }
+ TH2I* hNbDigiPerEvtTimeDet = new TH2I(Form("hNbDigiPerEvtTime%s", (*det).sName.data()),
+ Form("nb of %s digis per event vs seed time of the events; Seed time "
+ "[s]; Nb Digis []; Events []",
+ (*det).sName.data()),
+ 1000, 0, 0.001, 1000, 0, 100);
+ hNbDigiPerEvtTimeDet->SetCanExtend(TH2::kAllAxes);
+ fvhNbDigiPerEvtTimeDet.push_back(hNbDigiPerEvtTimeDet);
+ }
+
+ AddHistoToVector(fhEventTime, "evtbuild");
+ AddHistoToVector(fhEventDt, "evtbuild");
+ AddHistoToVector(fhEventSize, "evtbuild");
+ AddHistoToVector(fhNbDigiPerEvtTime, "evtbuild");
+ outFolder->Add(fhEventTime);
+ outFolder->Add(fhEventDt);
+ outFolder->Add(fhEventSize);
+ outFolder->Add(fhNbDigiPerEvtTime);
+
+ for (std::vector<TH2*>::iterator itHist = fvhNbDigiPerEvtTimeDet.begin(); itHist != fvhNbDigiPerEvtTimeDet.end();
+ ++itHist) {
+ if (nullptr != (*itHist)) {
+ AddHistoToVector((*itHist), "evtbuild");
+ outFolder->Add((*itHist));
+ }
+ }
+}
+
+void CbmEvbuildRawAlgo::FillHistos()
+{
+ Double_t dPreEvtTime = -1.0;
+ for (CbmRawEvent* evt : fEventVector) {
+ fhEventTime->Fill(evt->GetTime() * 1e-9);
+ if (0.0 <= dPreEvtTime) { fhEventDt->Fill((evt->GetTime() - dPreEvtTime) * 1e-9); }
+
+ fhEventSize->Fill(evt->GetNofDigis());
+ fhNbDigiPerEvtTime->Fill(evt->GetTime() * 1e-9, evt->GetNofDigis());
+
+ /// Loop on selection detectors
+ for (UInt_t uDetIdx = 0; uDetIdx < fvDets.size(); ++uDetIdx) {
+ if (nullptr == fvhNbDigiPerEvtTimeDet[uDetIdx]) continue;
+ fvhNbDigiPerEvtTimeDet[uDetIdx]->Fill(evt->GetTime() * 1e-9, evt->GetNofDigis(fvDets[uDetIdx].detId));
+ }
+ dPreEvtTime = evt->GetTime();
+ }
+}
+
+void CbmEvbuildRawAlgo::ResetHistograms(Bool_t /*bResetTime*/)
+{
+ fhEventTime->Reset();
+ fhEventDt->Reset();
+ fhEventSize->Reset();
+
+ fhNbDigiPerEvtTime->Reset();
+ /// Loop on histograms
+ for (std::vector<TH2*>::iterator itHist = fvhNbDigiPerEvtTimeDet.begin(); itHist != fvhNbDigiPerEvtTimeDet.end();
+ ++itHist) {
+ (*itHist)->Reset();
+ }
+ /*
+ if( kTRUE == bResetTime )
+ {
+ /// Also reset the Start time for the evolution plots!
+ fdStartTime = -1.0;
+ }
+ */
+}
+
+void CbmEvbuildRawAlgo::SetReferenceDetector(ECbmModuleId refDet, ECbmDataType dataTypeIn, std::string sNameIn,
+ UInt_t uTriggerMinDigisIn, Int_t iTriggerMaxDigisIn,
+ Double_t fdTimeWinBegIn, Double_t fdTimeWinEndIn)
+{
+ /// FIXME: Deprecated method to be removed later. For now create temp object.
+ SetReferenceDetector(RawEventBuilderDetector(refDet, dataTypeIn, sNameIn, uTriggerMinDigisIn, iTriggerMaxDigisIn,
+ fdTimeWinBegIn, fdTimeWinEndIn));
+}
+void CbmEvbuildRawAlgo::AddDetector(ECbmModuleId selDet, ECbmDataType dataTypeIn, std::string sNameIn,
+ UInt_t uTriggerMinDigisIn, Int_t iTriggerMaxDigisIn, Double_t fdTimeWinBegIn,
+ Double_t fdTimeWinEndIn)
+{
+ /// FIXME: Deprecated method to be removed later. For now create temp object.
+ AddDetector(RawEventBuilderDetector(selDet, dataTypeIn, sNameIn, uTriggerMinDigisIn, iTriggerMaxDigisIn,
+ fdTimeWinBegIn, fdTimeWinEndIn));
+}
+
+void CbmEvbuildRawAlgo::SetReferenceDetector(RawEventBuilderDetector refDetIn)
+{
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ if ((*det) == refDetIn) {
+ LOG(warning) << "CbmEvbuildRawAlgo::SetReferenceDetector => "
+ "Reference detector already in selection detector list!"
+ << refDetIn.sName;
+ LOG(warning) << " => "
+ "It will be automatically removed from selection detector list!";
+ LOG(warning) << " => "
+ "Please also remember to update the selection windows to store "
+ "clusters!";
+ RemoveDetector(refDetIn);
+ }
+ }
+
+ if (fRefDet == refDetIn) {
+ LOG(warning) << "CbmEvbuildRawAlgo::SetReferenceDetector => "
+ "Doing nothing, identical reference detector already in use";
+ }
+ else {
+ LOG(info) << "CbmEvbuildRawAlgo::SetReferenceDetector => "
+ << "Replacing " << fRefDet.sName << " with " << refDetIn.sName << " as reference detector";
+ LOG(warning) << " => "
+ "You may want to use AddDetector after this command to add in "
+ "selection "
+ << refDetIn.sName;
+ LOG(warning) << " => "
+ "Please also remember to update the selection windows!";
+ }
+ fRefDet = refDetIn;
+
+ /// Update the variables storing the earliest and latest time window boundaries
+ UpdateTimeWinBoundariesExtrema();
+ /// Update the variable storing the size if widest time window for overlap detection
+ UpdateWidestTimeWinRange();
+}
+
+void CbmEvbuildRawAlgo::AddDetector(RawEventBuilderDetector selDet)
+{
+ if (fRefDet == selDet) {
+ LOG(fatal) << "CbmEvbuildRawAlgo::AddDetector => Cannot "
+ "add the reference detector as selection detector!"
+ << std::endl
+ << "=> Maybe first change the reference detector with "
+ "SetReferenceDetector?";
+ }
+
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ if ((*det) == selDet) {
+ LOG(warning) << "CbmEvbuildRawAlgo::AddDetector => "
+ "Doing nothing, selection detector already in list!"
+ << selDet.sName;
+ return;
+ }
+ }
+ fvDets.push_back(selDet);
+
+ /// Update the variables storing the earliest and latest time window boundaries
+ UpdateTimeWinBoundariesExtrema();
+ /// Update the variable storing the size if widest time window for overlap detection
+ UpdateWidestTimeWinRange();
+}
+
+void CbmEvbuildRawAlgo::RemoveDetector(RawEventBuilderDetector selDet)
+{
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ if ((*det) == selDet) {
+ fvDets.erase(det);
+ return;
+ }
+ }
+ LOG(warning) << "CbmEvbuildRawAlgo::RemoveDetector => Doing "
+ "nothing, selection detector not in list!"
+ << selDet.sName;
+}
+
+void CbmEvbuildRawAlgo::SetTriggerMinNumber(ECbmModuleId selDet, UInt_t uVal)
+{
+ /// Check first if reference detector
+ if (fRefDet.detId == selDet) {
+ fRefDet.fuTriggerMinDigis = uVal;
+ LOG(debug) << "Set Trigger min limit for " << fRefDet.sName << " to " << uVal;
+ return;
+ }
+
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ if ((*det).detId == selDet) {
+ (*det).fuTriggerMinDigis = uVal;
+ LOG(debug) << "Set Trigger min limit for " << (*det).sName << " to " << uVal;
+ return;
+ }
+ }
+ LOG(warning) << "CbmEvbuildRawAlgo::SetTriggerMinNumber => "
+ "Doing nothing, detector neither reference nor in selection list!"
+ << selDet;
+}
+
+void CbmEvbuildRawAlgo::SetTriggerMaxNumber(ECbmModuleId selDet, Int_t iVal)
+{
+ /// Check first if reference detector
+ if (fRefDet.detId == selDet) {
+ fRefDet.fiTriggerMaxDigis = iVal;
+ LOG(debug) << "Set Trigger min limit for " << fRefDet.sName << " to " << iVal;
+ return;
+ }
+
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ if ((*det).detId == selDet) {
+ (*det).fiTriggerMaxDigis = iVal;
+ LOG(debug) << "Set Trigger min limit for " << (*det).sName << " to " << iVal;
+ return;
+ }
+ }
+ LOG(warning) << "CbmEvbuildRawAlgo::SetTriggerMaxNumber => "
+ "Doing nothing, detector neither reference nor in selection list!"
+ << selDet;
+}
+
+void CbmEvbuildRawAlgo::SetTriggerWindow(ECbmModuleId selDet, Double_t dWinBeg, Double_t dWinEnd)
+{
+ /// Check if valid time window: end strictly after beginning
+ if (dWinEnd <= dWinBeg) {
+ LOG(fatal) << "CbmEvbuildRawAlgo::SetTriggerWindow => "
+ "Invalid time window: [ "
+ << dWinBeg << ", " << dWinEnd << " ]";
+ }
+
+ Bool_t bFound = kFALSE;
+ /// Check first if reference detector
+ if (fRefDet.detId == selDet) {
+ fRefDet.fdTimeWinBeg = dWinBeg;
+ fRefDet.fdTimeWinEnd = dWinEnd;
+ bFound = kTRUE;
+ }
+
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ if ((*det).detId == selDet) {
+ (*det).fdTimeWinBeg = dWinBeg;
+ (*det).fdTimeWinEnd = dWinEnd;
+ bFound = kTRUE;
+ }
+ }
+
+ if (kFALSE == bFound) {
+ LOG(warning) << "CbmEvbuildRawAlgo::SetTriggerWindow => "
+ "Doing nothing, detector neither reference nor in selection list!"
+ << selDet;
+ }
+
+ /// Update the variables storing the earliest and latest time window boundaries
+ UpdateTimeWinBoundariesExtrema();
+ /// Update the variable storing the size if widest time window for overlap detection
+ UpdateWidestTimeWinRange();
+}
+
+void CbmEvbuildRawAlgo::UpdateTimeWinBoundariesExtrema()
+{
+ /// Initialize with reference detector or explicit times
+ if (fRefDet.detId != ECbmModuleId::kNotExist) {
+ fdEarliestTimeWinBeg = fRefDet.fdTimeWinBeg;
+ fdLatestTimeWinEnd = fRefDet.fdTimeWinEnd;
+ }
+ else {
+ fdEarliestTimeWinBeg = fdSeedTimeWinBeg;
+ fdLatestTimeWinEnd = fdSeedTimeWinEnd;
+ }
+
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ fdEarliestTimeWinBeg = std::min(fdEarliestTimeWinBeg, (*det).fdTimeWinBeg);
+ fdLatestTimeWinEnd = std::max(fdLatestTimeWinEnd, (*det).fdTimeWinEnd);
+ }
+}
+
+void CbmEvbuildRawAlgo::UpdateWidestTimeWinRange()
+{
+ /// Initialize with reference detector
+ fdWidestTimeWinRange = GetSeedTimeWinRange();
+
+ /// Loop on selection detectors
+ for (std::vector<RawEventBuilderDetector>::iterator det = fvDets.begin(); det != fvDets.end(); ++det) {
+ fdWidestTimeWinRange = std::max(fdWidestTimeWinRange, (*det).fdTimeWinEnd - (*det).fdTimeWinBeg);
+ }
+}
+
+ClassImp(CbmEvbuildRawAlgo)
diff --git a/reco/eventbuilder/digis/CbmEvbuildRawAlgo.h b/reco/eventbuilder/digis/CbmEvbuildRawAlgo.h
new file mode 100644
index 0000000000000000000000000000000000000000..cdf52673488994db4d237717c81e02f77a65b969
--- /dev/null
+++ b/reco/eventbuilder/digis/CbmEvbuildRawAlgo.h
@@ -0,0 +1,343 @@
+/********************************************************************************
+ * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
+ * *
+ * This software is distributed under the terms of the *
+ * GNU Lesser General Public Licence (LGPL) version 3, *
+ * copied verbatim in the file "LICENSE" *
+ ********************************************************************************/
+#ifndef CBMEVBUILDRAWALGO_H
+#define CBMEVBUILDRAWALGO_H
+
+/// CBM headers
+#include "CbmAlgoBuildRawEvents.h"
+#include "CbmMuchBeamTimeDigi.h"
+#include "CbmMuchDigi.h"
+#include "CbmPsdDigi.h"
+#include "CbmRichDigi.h"
+#include "CbmStsDigi.h"
+#include "CbmTofDigi.h"
+#include "CbmTrdDigi.h"
+
+/// FAIRROOT headers
+#include "FairTask.h"
+
+/// FAIRSOFT headers (geant, boost, ...)
+
+/// C/C++ headers
+#include <tuple>
+
+#include <boost/any.hpp>
+
+#include <array>
+#include <map>
+#include <set>
+#include <vector>
+
+class TimesliceMetaData;
+class CbmRawEvent;
+class TClonesArray;
+class TH1;
+class TH2;
+class TNamed;
+class TStopwatch;
+class TCanvas;
+
+/*
+enum class EOverlapModeRaw
+{
+ NoOverlap,
+ MergeOverlap,
+ AllowOverlap
+};
+*/
+
+/*
+class RawEventBuilderDetector {
+public:
+ RawEventBuilderDetector() = default;
+
+ RawEventBuilderDetector(ECbmModuleId detIdIn, ECbmDataType dataTypeIn, std::string sNameIn)
+ : detId {detIdIn}
+ , dataType {dataTypeIn}
+ , sName {sNameIn}
+ {
+ }
+
+ RawEventBuilderDetector(ECbmModuleId detIdIn, ECbmDataType dataTypeIn, std::string sNameIn, UInt_t uTriggerMinDigisIn,
+ Int_t iTriggerMaxDigisIn, Double_t fdTimeWinBegIn, Double_t fdTimeWinEndIn)
+ : RawEventBuilderDetector(detIdIn, dataTypeIn, sNameIn)
+ {
+ fuTriggerMinDigis = uTriggerMinDigisIn;
+ fiTriggerMaxDigis = iTriggerMaxDigisIn;
+ fdTimeWinBeg = fdTimeWinBegIn;
+ fdTimeWinEnd = fdTimeWinEndIn;
+ }
+
+ bool operator==(const RawEventBuilderDetector& other) const { return (other.detId == this->detId); }
+ bool operator!=(const RawEventBuilderDetector& other) const { return (other.detId != this->detId); }
+
+ Double_t GetTimeWinRange() { return fdTimeWinEnd - fdTimeWinBeg; }
+
+ /// Settings
+ ECbmModuleId detId = ECbmModuleId::kNotExist;
+ ECbmDataType dataType = ECbmDataType::kUnknown;
+ std::string sName = "Invalid";
+ /// Minimum number of T0 digis needed to generate a trigger, 0 means don't use for trigger generation
+ UInt_t fuTriggerMinDigis = 0;
+ /// Maximum number of digis per detector to generate an event, -1 means no cut, 0 means anti-coinc trigger
+ Int_t fiTriggerMaxDigis = -1;
+ /// Selection Window
+ Double_t fdTimeWinBeg = -100;
+ Double_t fdTimeWinEnd = 100;
+ /// Book-keeping variables
+ UInt_t fuStartIndex = 0;
+ UInt_t fuEndIndex = 0;
+};
+
+/// Pre-defined detector types
+static const RawEventBuilderDetector kRawEventBuilderDetSts =
+ RawEventBuilderDetector(ECbmModuleId::kSts, ECbmDataType::kStsDigi, "Sts");
+static const RawEventBuilderDetector kRawEventBuilderDetMuch =
+ RawEventBuilderDetector(ECbmModuleId::kMuch, ECbmDataType::kMuchDigi, "Much");
+static const RawEventBuilderDetector kRawEventBuilderDetTrd =
+ RawEventBuilderDetector(ECbmModuleId::kTrd, ECbmDataType::kTrdDigi, "Trd");
+static const RawEventBuilderDetector kRawEventBuilderDetTof =
+ RawEventBuilderDetector(ECbmModuleId::kTof, ECbmDataType::kTofDigi, "Tof");
+static const RawEventBuilderDetector kRawEventBuilderDetRich =
+ RawEventBuilderDetector(ECbmModuleId::kRich, ECbmDataType::kRichDigi, "Rich");
+static const RawEventBuilderDetector kRawEventBuilderDetPsd =
+ RawEventBuilderDetector(ECbmModuleId::kPsd, ECbmDataType::kPsdDigi, "Psd");
+static const RawEventBuilderDetector kRawEventBuilderDetT0 =
+ RawEventBuilderDetector(ECbmModuleId::kT0, ECbmDataType::kT0Digi, "T0");
+static const RawEventBuilderDetector kRawEventBuilderDetUndef = RawEventBuilderDetector();
+
+*/
+
+class CbmEvbuildRawAlgo {
+public:
+ /** Default constructor **/
+ CbmEvbuildRawAlgo() = default;
+
+ CbmEvbuildRawAlgo(const CbmEvbuildRawAlgo&) = delete;
+ CbmEvbuildRawAlgo operator=(const CbmEvbuildRawAlgo&) = delete;
+
+ /** Destructor **/
+ ~CbmEvbuildRawAlgo() {};
+
+ /** Initiliazation at the beginning of a run **/
+ Bool_t InitAlgo();
+
+ /** Executed for each TS. **/
+ void ProcessTs();
+
+ /** Finish called at the end of the run **/
+ void Finish();
+
+ void SetFillHistos(Bool_t var) { fbFillHistos = var; }
+ void ResetHistograms(Bool_t bResetTime = kTRUE);
+
+ /** stopwatch timing **/
+ void SetTimings(Bool_t var) { fbGetTimings = var; }
+ void PrintTimings();
+
+ void SetReferenceDetector(ECbmModuleId refDet, ECbmDataType dataTypeIn, std::string sNameIn,
+ UInt_t uTriggerMinDigisIn = 0, Int_t iTriggerMaxDigisIn = -1,
+ Double_t fdTimeWinBegIn = -100, Double_t fdTimeWinEndIn = 100);
+ void AddDetector(ECbmModuleId selDet, ECbmDataType dataTypeIn, std::string sNameIn, UInt_t uTriggerMinDigisIn = 0,
+ Int_t iTriggerMaxDigisIn = -1, Double_t fdTimeWinBegIn = -100, Double_t fdTimeWinEndIn = 100);
+
+ void SetReferenceDetector(RawEventBuilderDetector refDetIn);
+ void AddDetector(RawEventBuilderDetector selDet);
+ void RemoveDetector(RawEventBuilderDetector selDet);
+
+ void SetTriggerMinNumber(ECbmModuleId selDet, UInt_t uVal);
+ void SetTriggerMaxNumber(ECbmModuleId selDet, Int_t iVal);
+ void SetTriggerWindow(ECbmModuleId selDet, Double_t dWinBeg, Double_t dWinEnd);
+
+ void SetTsParameters(Double_t dTsStartTime, Double_t dTsLength, Double_t dTsOverLength)
+ {
+ fdTsStartTime = dTsStartTime;
+ fdTsLength = dTsLength;
+ fdTsOverLength = dTsOverLength;
+ fbUseTsMetaData = kFALSE;
+ }
+
+ void SetSeedTimeWindow(Double_t timeWinBeg, Double_t timeWinEnd)
+ {
+ fdSeedTimeWinBeg = timeWinBeg;
+ fdSeedTimeWinEnd = timeWinEnd;
+ UpdateTimeWinBoundariesExtrema();
+ UpdateWidestTimeWinRange();
+ }
+
+ /// Control flags
+ void SetEventOverlapMode(EOverlapModeRaw mode) { fOverMode = mode; }
+ void SetIgnoreTsOverlap(Bool_t bFlagIn = kTRUE) { fbIgnoreTsOverlap = bFlagIn; }
+
+ void ChangeMuchBeamtimeDigiFlag(Bool_t bFlagIn = kFALSE) { fbUseMuchBeamtimeDigi = bFlagIn; }
+
+ /// For monitor algos
+ void AddHistoToVector(TNamed* pointer, std::string sFolder = "")
+ {
+ fvpAllHistoPointers.push_back(std::pair<TNamed*, std::string>(pointer, sFolder));
+ }
+ std::vector<std::pair<TNamed*, std::string>> GetHistoVector() { return fvpAllHistoPointers; }
+ void AddCanvasToVector(TCanvas* pointer, std::string sFolder = "")
+ {
+ fvpAllCanvasPointers.push_back(std::pair<TCanvas*, std::string>(pointer, sFolder));
+ }
+ std::vector<std::pair<TCanvas*, std::string>> GetCanvasVector() { return fvpAllCanvasPointers; }
+
+ /// Set digi containers
+ void SetT0Digis(const std::vector<CbmTofDigi>* T0DigiVec) { fT0DigiVec = T0DigiVec; }
+ void SetStsDigis(std::vector<CbmStsDigi>* StsDigis) { fStsDigis = StsDigis; }
+ void SetMuchDigis(std::vector<CbmMuchDigi>* MuchDigis) { fMuchDigis = MuchDigis; }
+ void SetTrdDigis(std::vector<CbmTrdDigi>* TrdDigis) { fTrdDigis = TrdDigis; }
+ void SetTofDigis(std::vector<CbmTofDigi>* TofDigis) { fTofDigis = TofDigis; }
+ void SetRichDigis(std::vector<CbmRichDigi>* RichDigis) { fRichDigis = RichDigis; }
+ void SetPsdDigis(std::vector<CbmPsdDigi>* PsdDigis) { fPsdDigis = PsdDigis; }
+ void SetMuchBeamTimeDigis(std::vector<CbmMuchBeamTimeDigi>* MuchBeamTimeDigis)
+ {
+ fMuchBeamTimeDigis = MuchBeamTimeDigis;
+ }
+
+ void SetSeedTimes(std::vector<Double_t>* SeedTimes) { fSeedTimes = SeedTimes; }
+
+ // TS metadata
+ void SetTimeSliceMetaDataArray(TClonesArray* TimeSliceMetaDataArray)
+ {
+ fTimeSliceMetaDataArray = TimeSliceMetaDataArray;
+ }
+
+ // Output folder for histograms
+ TFolder* GetOutFolder() { return outFolder; }
+
+ /// Data output access
+ std::vector<CbmRawEvent*>& GetEventVector() { return fEventVector; }
+ void ClearEventVector();
+
+private:
+ /// Internal methods
+ Bool_t CheckDataAvailable(RawEventBuilderDetector& det);
+ void InitTs();
+ void InitSeedWindow();
+ void BuildEvents();
+
+ void CreateHistograms();
+ void FillHistos();
+
+ template<class DigiSeed>
+ void LoopOnSeeds();
+
+ void CheckSeed(Double_t dSeedTime, UInt_t uSeedDigiIdx);
+ void CheckTriggerCondition(Double_t dSeedTime);
+
+ template<class DigiCheck>
+ void SearchMatches(Double_t dSeedTime, RawEventBuilderDetector& detMatch);
+ void SearchMatches(Double_t dSeedTime, RawEventBuilderDetector& detMatch);
+ void AddDigiToEvent(RawEventBuilderDetector& det, Int_t uIdx);
+ Bool_t HasTrigger(CbmRawEvent*);
+ Bool_t CheckTriggerConditions(CbmRawEvent* event, RawEventBuilderDetector& det);
+
+ void UpdateTimeWinBoundariesExtrema();
+ void UpdateWidestTimeWinRange();
+
+ TFolder* outFolder; // oputput folder to store histograms
+
+ /// Constants
+ static constexpr Double_t kdDefaultTimeWinBeg = -100.0;
+ static constexpr Double_t kdDefaultTimeWinEnd = 100.0;
+
+ /// User parameters
+ /// Control flags
+ Bool_t fbIgnoreTsOverlap = kFALSE; //! Ignore data in Overlap part of the TS
+ Bool_t fbFillHistos {kTRUE}; //! Switch ON/OFF filling of histograms
+ Bool_t fbUseMuchBeamtimeDigi = kTRUE; //! Switch between MUCH digi classes
+ Bool_t fbGetTimings = kFALSE; //! Measure CPU time using stopwatch
+ Bool_t fbUseTsMetaData = kTRUE; //! Read Ts Parameters from input tree
+ /// Event building mode and detectors selection
+ EOverlapModeRaw fOverMode {EOverlapModeRaw::AllowOverlap};
+
+ TStopwatch* fTimer = nullptr; //! is create when fbGetTimings is set before init
+
+ RawEventBuilderDetector fRefDet = RawEventBuilderDetector(ECbmModuleId::kT0, ECbmDataType::kT0Digi, "T0");
+ std::vector<RawEventBuilderDetector> fvDets = {
+ RawEventBuilderDetector(ECbmModuleId::kSts, ECbmDataType::kStsDigi, "kSts"),
+ RawEventBuilderDetector(ECbmModuleId::kMuch, ECbmDataType::kMuchDigi, "kMuch"),
+ RawEventBuilderDetector(ECbmModuleId::kTrd, ECbmDataType::kTrdDigi, "kTrd"),
+ RawEventBuilderDetector(ECbmModuleId::kTof, ECbmDataType::kTofDigi, "kTof"),
+ RawEventBuilderDetector(ECbmModuleId::kRich, ECbmDataType::kRichDigi, "kRich"),
+ RawEventBuilderDetector(ECbmModuleId::kPsd, ECbmDataType::kPsdDigi, "kPsd")};
+
+ Double_t fdEarliestTimeWinBeg = kdDefaultTimeWinBeg;
+ Double_t fdLatestTimeWinEnd = kdDefaultTimeWinEnd;
+ Double_t fdWidestTimeWinRange = kdDefaultTimeWinEnd - kdDefaultTimeWinBeg;
+ ///Seed window
+ Double_t fdSeedWindowBeg = 0;
+ Double_t fdSeedWindowEnd = 0;
+
+ Double_t fdTsStartTime = -1;
+ Double_t fdTsLength = -1;
+ Double_t fdTsOverLength = -1;
+
+ /// Data input
+ TClonesArray* fTimeSliceMetaDataArray = nullptr; //!
+
+ const std::vector<CbmTofDigi>* fT0DigiVec = nullptr;
+ const std::vector<CbmMuchDigi>* fMuchDigis = nullptr;
+ const std::vector<CbmMuchBeamTimeDigi>* fMuchBeamTimeDigis = nullptr;
+ const std::vector<CbmStsDigi>* fStsDigis = nullptr;
+ const std::vector<CbmTrdDigi>* fTrdDigis = nullptr;
+ const std::vector<CbmTofDigi>* fTofDigis = nullptr;
+ const std::vector<CbmRichDigi>* fRichDigis = nullptr;
+ const std::vector<CbmPsdDigi>* fPsdDigis = nullptr;
+
+ // If explicit seed times are supplied.
+ const std::vector<Double_t>* fSeedTimes = nullptr;
+ Double_t fdSeedTimeWinBeg = -100.0;
+ Double_t fdSeedTimeWinEnd = 100.0;
+
+ bool CheckDataAvailable(ECbmModuleId detId);
+ UInt_t GetNofDigis(ECbmModuleId detId);
+ template<class Digi>
+ const Digi* GetDigi(UInt_t uDigi);
+
+ Double_t GetSeedTimeWinRange();
+
+ /// Data ouptut
+ CbmRawEvent* fCurrentEvent = nullptr; //! pointer to the event which is currently build
+ std::vector<CbmRawEvent*> fEventVector = {}; //! vector with all created events
+
+ /// Monitoring histograms
+ /// => Pointers should be filled with TH1*, TH2*, TProfile*, ...
+ /// ==> To check if object N is of type T, use "T ObjectPointer = dynamic_cast<T>( fvpAllHistoPointers[N].first );" and check for nullptr
+ /// ==> To get back the original class name use "fvpAllHistoPointers[N].first->ClassName()" which returns a const char * (e.g. "TH1I")
+ /// ===> Usage example with feeding a THttpServer:
+ /// ===> #include "TH2.h"
+ /// ===> std::string sClassName = vHistos[ uHisto ].first.ClassName();
+ /// ===> if( !strncmp( sClassName, "TH1", 3 ) )
+ /// ===> server->Register( vHistos[ uHisto ].second.data(), dynamic_cast< TH1 * >(vHistos[ uHisto ].first) );
+ /// ===> else if( !strncmp( sClassName, "TH2", 3 ) )
+ /// ===> server->Register( vHistos[ uHisto ].second.data(), dynamic_cast< TH2 * >(vHistos[ uHisto ].first) );
+ std::vector<std::pair<TNamed*, std::string>>
+ fvpAllHistoPointers; //! Vector of pointers to histograms + optional folder name
+ std::vector<std::pair<TCanvas*, std::string>>
+ fvpAllCanvasPointers; //! Vector of pointers to canvases + optional folder name
+
+ TH1* fhEventTime = nullptr; //! histogram with the seed time of the events
+ TH1* fhEventDt = nullptr; //! histogram with the interval in seed time of consecutive events
+ TH1* fhEventSize = nullptr; //! histogram with the nb of all digis in the event
+ TH2* fhNbDigiPerEvtTime = nullptr; //! histogram with the nb of all digis per event vs seed time of the events
+ std::vector<TH2*> fvhNbDigiPerEvtTimeDet =
+ {}; //! histograms with the nb of digis in each detector per event vs seed time of the events
+
+ /// Internal state variables
+ UInt_t fuCurEv = 0; //! Event Counter
+ UInt_t fuNrTs = 0; //! Timeslice Counter
+ Double_t fdPrevEvtTime = 0.; //! Save previous time information
+ Double_t fdPrevEvtEndTime = 0.; //! Save previous event last digi time information
+
+ ClassDefNV(CbmEvbuildRawAlgo, 1);
+};
+
+#endif // CBMEVBUILDRAWALGO_H
diff --git a/reco/eventbuilder/digis/CbmEvbuildRawTask.cxx b/reco/eventbuilder/digis/CbmEvbuildRawTask.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..729a2c9522157e4fc4a98320e3664d7f5e89a529
--- /dev/null
+++ b/reco/eventbuilder/digis/CbmEvbuildRawTask.cxx
@@ -0,0 +1,500 @@
+/********************************************************************************
+ * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
+ * *
+ * This software is distributed under the terms of the *
+ * GNU Lesser General Public Licence (LGPL) version 3, *
+ * copied verbatim in the file "LICENSE" *
+ ********************************************************************************/
+#include "CbmEvbuildRawTask.h"
+
+#include "CbmDigiManager.h"
+#include "CbmRawEvent.h"
+
+#include <FairLogger.h>
+#include <FairRootManager.h>
+#include <FairRunOnline.h>
+
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TFolder.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <THttpServer.h>
+#include <TStopwatch.h>
+
+#include <iomanip>
+
+using std::cout;
+using std::endl;
+using std::fixed;
+using std::left;
+using std::pair;
+using std::right;
+using std::setprecision;
+using std::setw;
+using std::stringstream;
+
+
+CbmEvbuildRawTask::CbmEvbuildRawTask() : FairTask("BuildRawEvents")
+{
+ /// Create Algo. To be made generic/switchable when more event building algo are available!
+ fpAlgo = new CbmEvbuildRawAlgo();
+}
+
+CbmEvbuildRawTask::~CbmEvbuildRawTask() {}
+
+void CbmEvbuildRawTask::SetParContainers() {}
+
+void CbmEvbuildRawTask::SetSeedTimeFiller(RawEventBuilderDetector seedDet)
+{
+ if (fSeedTimeDetList.size() > 0) { LOG(fatal) << "Cannot use seed detector list and single instance together."; }
+
+ fSeedTimeDet = seedDet;
+ if (fSeedTimeDet != kRawEventBuilderDetUndef) {
+ if (fSeedTimes == nullptr) { fSeedTimes = new std::vector<Double_t>; }
+ }
+ else {
+ if (fSeedTimes != nullptr) {
+ fSeedTimes->clear();
+ delete fSeedTimes;
+ fSeedTimes = nullptr;
+ }
+ }
+ fpAlgo->SetSeedTimes(fSeedTimes);
+}
+
+void CbmEvbuildRawTask::AddSeedTimeFillerToList(RawEventBuilderDetector seedDet)
+{
+ if (fSeedTimeDet != kRawEventBuilderDetUndef) {
+ LOG(fatal) << "Cannot use seed detector list and single instance together.";
+ }
+ if (fSeedTimes == nullptr) { fSeedTimes = new std::vector<Double_t>; }
+
+ fSeedTimeDetList.push_back(seedDet);
+ fpAlgo->SetSeedTimes(fSeedTimes);
+}
+
+InitStatus CbmEvbuildRawTask::Init()
+{
+ if (fbGetTimings) {
+ fTimer = new TStopwatch;
+ fTimer->Start();
+ fCopyTimer = new TStopwatch;
+ fCopyTimer->Reset();
+ }
+
+ /// Get a handle from the IO manager
+ FairRootManager* ioman = FairRootManager::Instance();
+
+ //T0 not included in digi manager.
+ fT0Digis = ioman->InitObjectAs<std::vector<CbmTofDigi> const*>("T0Digi");
+ if (!fT0Digis) { LOG(info) << "No T0 digi input."; }
+ else {
+ LOG(info) << "T0 digi input.";
+ fpAlgo->SetT0Digis(fT0Digis);
+ }
+
+ // Get a pointer to the previous already existing data level
+ fDigiMan = CbmDigiManager::Instance();
+ if (fbUseMuchBeamtimeDigi) { fDigiMan->UseMuchBeamTimeDigi(); }
+ fDigiMan->Init();
+
+ //Init STS digis
+ if (!fDigiMan->IsPresent(ECbmModuleId::kSts)) { LOG(info) << "No STS digi input."; }
+ else {
+ LOG(info) << "STS digi input.";
+ fStsDigis = new std::vector<CbmStsDigi>;
+ fpAlgo->SetStsDigis(fStsDigis);
+ }
+
+ //Init MUCH digis
+ if (!fDigiMan->IsPresent(ECbmModuleId::kMuch)) { LOG(info) << "No MUCH digi input."; }
+ else {
+ LOG(info) << "MUCH digi input.";
+ if (fbUseMuchBeamtimeDigi) {
+ fMuchBeamTimeDigis = new std::vector<CbmMuchBeamTimeDigi>;
+ fpAlgo->SetMuchBeamTimeDigis(fMuchBeamTimeDigis);
+ }
+ else {
+ fMuchDigis = new std::vector<CbmMuchDigi>;
+ fpAlgo->SetMuchDigis(fMuchDigis);
+ }
+ }
+
+ //Init TRD digis
+ if (!fDigiMan->IsPresent(ECbmModuleId::kTrd)) { LOG(info) << "No TRD digi input."; }
+ else {
+ LOG(info) << "TRD digi input.";
+ fTrdDigis = new std::vector<CbmTrdDigi>;
+ fpAlgo->SetTrdDigis(fTrdDigis);
+ }
+
+ //Init TOF digis
+ if (!fDigiMan->IsPresent(ECbmModuleId::kTof)) { LOG(info) << "No TOF digi input."; }
+ else {
+ LOG(info) << "TOF digi input.";
+ fTofDigis = new std::vector<CbmTofDigi>;
+ fpAlgo->SetTofDigis(fTofDigis);
+ }
+
+ //Init RICH digis
+ if (!fDigiMan->IsPresent(ECbmModuleId::kRich)) { LOG(info) << "No RICH digi input."; }
+ else {
+ LOG(info) << "RICH digi input.";
+ fRichDigis = new std::vector<CbmRichDigi>;
+ fpAlgo->SetRichDigis(fRichDigis);
+ }
+
+ //Init PSD digis
+ if (!fDigiMan->IsPresent(ECbmModuleId::kPsd)) { LOG(info) << "No PSD digi input."; }
+ else {
+ LOG(info) << "PSD digi input.";
+ fPsdDigis = new std::vector<CbmPsdDigi>;
+ fpAlgo->SetPsdDigis(fPsdDigis);
+ }
+
+ /// Register output array (CbmEvent)
+ fEvents = new TClonesArray("CbmRawEvent", 100);
+ ioman->Register("RawEvent", "CBM Raw Event", fEvents, IsOutputBranchPersistent("RawEvent"));
+ if (!fEvents) LOG(fatal) << "Output branch was not created";
+
+ // Set timeslice meta data
+ fpAlgo->SetTimeSliceMetaDataArray(dynamic_cast<TClonesArray*>(ioman->GetObject("TimesliceMetaData")));
+
+ if (fTimer != nullptr) {
+ fTimer->Stop();
+ Double_t rtime = fTimer->RealTime();
+ Double_t ctime = fTimer->CpuTime();
+ LOG(info) << "CbmEvbuildRawTask::Init(): Real time " << rtime << " s, CPU time " << ctime << " s";
+ }
+
+ /// Call Algo Init method
+ if (kTRUE == fpAlgo->InitAlgo()) return kSUCCESS;
+ else
+ return kFATAL;
+}
+
+
+InitStatus CbmEvbuildRawTask::ReInit() { return kSUCCESS; }
+
+void CbmEvbuildRawTask::Exec(Option_t* /*option*/)
+{
+ if (fTimer != nullptr) { fTimer->Start(kFALSE); }
+ TStopwatch timer;
+ timer.Start();
+
+ LOG(debug2) << "CbmEvbuildRawTask::Exec => Starting sequence";
+ //Warning: Int_t must be used for the loop counters instead of UInt_t,
+ //as the digi manager can return -1, which would be casted to +1
+ //during comparison, leading to an error.
+
+ if (fSeedTimeDet != kRawEventBuilderDetUndef && fSeedTimeDetList.size() > 0) {
+ LOG(fatal) << "Cannot use seed detector list and single instance together.";
+ }
+
+ //Reset explicit seed times if set
+ if (fSeedTimeDet != kRawEventBuilderDetUndef || fSeedTimeDetList.size() > 0) { fSeedTimes->clear(); }
+
+ if (fCopyTimer != nullptr) { fCopyTimer->Start(kFALSE); }
+
+ //Read STS digis
+ if (fDigiMan->IsPresent(ECbmModuleId::kSts)) {
+ fStsDigis->clear();
+ for (Int_t i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kSts); i++) {
+ const CbmStsDigi* Digi = fDigiMan->Get<CbmStsDigi>(i);
+ fStsDigis->push_back(*Digi);
+ if (fSeedTimeDet.detId == ECbmModuleId::kSts) { fSeedTimes->push_back(Digi->GetTime()); }
+ }
+ LOG(debug) << "Read: " << fStsDigis->size() << " STS digis.";
+ LOG(debug) << "In DigiManager: " << fDigiMan->GetNofDigis(ECbmModuleId::kSts) << " STS digis.";
+ }
+
+ //Read MUCH digis
+ if (fDigiMan->IsPresent(ECbmModuleId::kMuch)) {
+ if (fbUseMuchBeamtimeDigi) {
+ fMuchBeamTimeDigis->clear();
+ for (Int_t i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kMuch); i++) {
+ const CbmMuchBeamTimeDigi* Digi = fDigiMan->Get<CbmMuchBeamTimeDigi>(i);
+ fMuchBeamTimeDigis->push_back(*Digi);
+ if (fSeedTimeDet.detId == ECbmModuleId::kMuch) { fSeedTimes->push_back(Digi->GetTime()); }
+ }
+ LOG(debug) << "Read: " << fDigiMan->GetNofDigis(ECbmModuleId::kMuch) << " MUCH digis.";
+ LOG(debug) << "In DigiManager: " << fMuchBeamTimeDigis->size() << " MUCH digis.";
+ }
+ else {
+ fMuchDigis->clear();
+ for (Int_t i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kMuch); i++) {
+ const CbmMuchDigi* Digi = fDigiMan->Get<CbmMuchDigi>(i);
+ fMuchDigis->push_back(*Digi);
+ if (fSeedTimeDet.detId == ECbmModuleId::kMuch) { fSeedTimes->push_back(Digi->GetTime()); }
+ }
+ LOG(debug) << "Read: " << fDigiMan->GetNofDigis(ECbmModuleId::kMuch) << " MUCH digis.";
+ LOG(debug) << "In DigiManager: " << fMuchDigis->size() << " MUCH digis.";
+ }
+ }
+
+ //Read TRD digis
+ if (fDigiMan->IsPresent(ECbmModuleId::kTrd)) {
+ fTrdDigis->clear();
+ for (Int_t i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kTrd); i++) {
+ const CbmTrdDigi* Digi = fDigiMan->Get<CbmTrdDigi>(i);
+ fTrdDigis->push_back(*Digi);
+ if (fSeedTimeDet.detId == ECbmModuleId::kTrd) { fSeedTimes->push_back(Digi->GetTime()); }
+ }
+ LOG(debug) << "Read: " << fDigiMan->GetNofDigis(ECbmModuleId::kTrd) << " TRD digis.";
+ LOG(debug) << "In DigiManager: " << fTrdDigis->size() << " TRD digis.";
+ }
+
+ //Read TOF digis
+ if (fDigiMan->IsPresent(ECbmModuleId::kTof)) {
+ fTofDigis->clear();
+ for (Int_t i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kTof); i++) {
+ const CbmTofDigi* Digi = fDigiMan->Get<CbmTofDigi>(i);
+ fTofDigis->push_back(*Digi);
+ if (fSeedTimeDet.detId == ECbmModuleId::kTof) { fSeedTimes->push_back(Digi->GetTime()); }
+ }
+ LOG(debug) << "Read: " << fDigiMan->GetNofDigis(ECbmModuleId::kTof) << " TOF digis.";
+ LOG(debug) << "In DigiManager: " << fTofDigis->size() << " TOF digis.";
+ }
+
+ //Read RICH digis
+ if (fDigiMan->IsPresent(ECbmModuleId::kRich)) {
+ fRichDigis->clear();
+ for (Int_t i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kRich); i++) {
+ const CbmRichDigi* Digi = fDigiMan->Get<CbmRichDigi>(i);
+ fRichDigis->push_back(*Digi);
+ if (fSeedTimeDet.detId == ECbmModuleId::kRich) { fSeedTimes->push_back(Digi->GetTime()); }
+ }
+ LOG(debug) << "Read: " << fDigiMan->GetNofDigis(ECbmModuleId::kRich) << " RICH digis.";
+ LOG(debug) << "In DigiManager: " << fRichDigis->size() << " RICH digis.";
+ }
+
+ //Read PSD digis
+ if (fDigiMan->IsPresent(ECbmModuleId::kPsd)) {
+ fPsdDigis->clear();
+ for (Int_t i = 0; i < fDigiMan->GetNofDigis(ECbmModuleId::kPsd); i++) {
+ const CbmPsdDigi* Digi = fDigiMan->Get<CbmPsdDigi>(i);
+ fPsdDigis->push_back(*Digi);
+ if (fSeedTimeDet.detId == ECbmModuleId::kPsd) { fSeedTimes->push_back(Digi->GetTime()); }
+ }
+ LOG(debug) << "Read: " << fDigiMan->GetNofDigis(ECbmModuleId::kPsd) << " PSD digis.";
+ LOG(debug) << "In DigiManager: " << fPsdDigis->size() << " PSD digis.";
+ }
+
+ if (fCopyTimer != nullptr) { fCopyTimer->Stop(); }
+
+ if (fSeedTimeDetList.size() > 0) { FillSeedTimesFromDetList(); }
+ //DumpSeedTimesFromDetList();
+
+ /// Call Algo ProcessTs method
+ fpAlgo->ProcessTs();
+
+ /// Save the resulting vector of events in TClonesArray
+ FillOutput();
+ LOG(debug2) << "CbmEvbuildRawTask::Exec => Done";
+
+ if (fTimer != nullptr) { fTimer->Stop(); }
+
+ // --- Timeslice log and statistics
+ timer.Stop();
+ stringstream logOut;
+ logOut << setw(20) << left << GetName() << " [";
+ logOut << fixed << setw(8) << setprecision(1) << right << timer.RealTime() * 1000. << " ms] ";
+ logOut << "TS " << fNofTs;
+ if (fEvents) logOut << ", events " << fEvents->GetEntriesFast();
+ LOG(info) << logOut.str();
+ fNofTs++;
+ fNofEvents += fEvents->GetEntriesFast();
+ fTime += timer.RealTime();
+}
+
+void CbmEvbuildRawTask::FillSeedTimesFromDetList()
+{
+ std::map<ECbmModuleId, UInt_t> DigiNumbers;
+ std::map<ECbmModuleId, UInt_t> DigiCounters;
+ fSeedTimes->clear();
+
+ for (RawEventBuilderDetector& system : fSeedTimeDetList) {
+ DigiNumbers[system.detId] = GetNofDigis(system.detId);
+ DigiCounters[system.detId] = 0;
+ }
+
+ do {
+ ECbmModuleId nextAddedSystem;
+ Double_t earliestTime = -1;
+
+ for (RawEventBuilderDetector& system : fSeedTimeDetList) {
+ if (DigiCounters[system.detId] < DigiNumbers[system.detId]) {
+ Double_t thisTime = GetDigiTime(system.detId, DigiCounters[system.detId]);
+ if (thisTime < earliestTime || earliestTime == -1) {
+ nextAddedSystem = system.detId;
+ earliestTime = thisTime;
+ }
+ }
+ }
+ if (earliestTime != -1) {
+ fSeedTimes->push_back(earliestTime);
+ DigiCounters[nextAddedSystem]++;
+ }
+ else {
+ break;
+ }
+ } while (true);
+}
+
+void CbmEvbuildRawTask::DumpSeedTimesFromDetList()
+{
+ std::ofstream timesUnsorted("digiTimesUnsorted.dat", std::ofstream::out);
+ timesUnsorted << std::setprecision(16);
+
+ for (RawEventBuilderDetector& system : fSeedTimeDetList) {
+ for (UInt_t i = 0; i < GetNofDigis(system.detId); i++) {
+ timesUnsorted << GetDigiTime(system.detId, i) << std::endl;
+ }
+ }
+ timesUnsorted.close();
+ LOG(info) << "Completed write of unsorted digi list.";
+
+ std::ofstream timesSorted("digiTimesSorted.dat", std::ofstream::out);
+ timesSorted << std::setprecision(16);
+
+ for (UInt_t i = 0; i < fSeedTimes->size(); i++) {
+ timesSorted << fSeedTimes->at(i) << std::endl;
+ }
+ timesSorted.close();
+ LOG(info) << "Completed DumpSeedTimesFromDetList(). Closing.";
+ exit(0); //terminate as this method should only be used for diagnostics
+}
+
+Double_t CbmEvbuildRawTask::GetDigiTime(ECbmModuleId _system, UInt_t _entry)
+{
+ switch (_system) {
+ case ECbmModuleId::kMuch:
+ if (fbUseMuchBeamtimeDigi) { return (fMuchBeamTimeDigis->at(_entry)).GetTime(); }
+ else {
+ return (fMuchDigis->at(_entry)).GetTime();
+ }
+ case ECbmModuleId::kSts: return (fStsDigis->at(_entry)).GetTime();
+ case ECbmModuleId::kTrd: return (fTrdDigis->at(_entry)).GetTime();
+ case ECbmModuleId::kTof: return (fTofDigis->at(_entry)).GetTime();
+ case ECbmModuleId::kRich: return (fRichDigis->at(_entry)).GetTime();
+ case ECbmModuleId::kPsd: return (fPsdDigis->at(_entry)).GetTime();
+ case ECbmModuleId::kHodo: return (fT0Digis->at(_entry)).GetTime();
+ default: break;
+ }
+ return -1;
+}
+
+UInt_t CbmEvbuildRawTask::GetNofDigis(ECbmModuleId _system)
+{
+ switch (_system) {
+ case ECbmModuleId::kMuch:
+ if (fbUseMuchBeamtimeDigi) { return fMuchBeamTimeDigis->size(); }
+ else {
+ return fMuchDigis->size();
+ }
+ case ECbmModuleId::kSts: return fStsDigis->size();
+ case ECbmModuleId::kTrd: return fTrdDigis->size();
+ case ECbmModuleId::kTof: return fTofDigis->size();
+ case ECbmModuleId::kRich: return fRichDigis->size();
+ case ECbmModuleId::kPsd: return fPsdDigis->size();
+ case ECbmModuleId::kHodo: return fT0Digis->size();
+ default: break;
+ }
+ return 0;
+}
+
+void CbmEvbuildRawTask::PrintTimings()
+{
+ if (fTimer == nullptr) { LOG(fatal) << "Trying to print timings but timer not set"; }
+ else {
+ Double_t rtime = fTimer->RealTime();
+ Double_t ctime = fTimer->CpuTime();
+ LOG(info) << "CbmEvbuildRawTask: Real time " << rtime << " s, CPU time " << ctime << " s";
+ }
+ if (fCopyTimer == nullptr) { LOG(fatal) << "Trying to print timings but timer not set"; }
+ else {
+ Double_t rtime = fCopyTimer->RealTime();
+ Double_t ctime = fCopyTimer->CpuTime();
+ LOG(info) << "CbmEvbuildRawTask (digi copy only): Real time " << rtime << " s, CPU time " << ctime << " s";
+ }
+}
+
+void CbmEvbuildRawTask::Finish()
+{
+ /// Call Algo finish method
+ fpAlgo->Finish();
+ if (fbFillHistos) { SaveHistos(); }
+ if (fbGetTimings) { PrintTimings(); }
+
+ std::cout << std::endl;
+ LOG(info) << "=====================================";
+ LOG(info) << GetName() << ": Run summary";
+ LOG(info) << "Time slices : " << fNofTs;
+ LOG(info) << "Events : " << fNofEvents;
+ LOG(info) << "Time / TS : " << fixed << setprecision(2) << 1000. * fTime / Double_t(fNofTs) << " ms";
+ LOG(info) << "=====================================";
+}
+
+void CbmEvbuildRawTask::FillOutput()
+{
+ /// Clear TClonesArray before usage.
+ fEvents->Delete();
+
+ /// Get vector reference from algo
+ std::vector<CbmRawEvent*> vEvents = fpAlgo->GetEventVector();
+
+ /// Move CbmEvent from temporary vector to TClonesArray
+ for (CbmRawEvent* event : vEvents) {
+ LOG(debug) << "Vector: " << event->ToString();
+ new ((*fEvents)[fEvents->GetEntriesFast()]) CbmRawEvent(std::move(*event));
+ LOG(debug) << "TClonesArray: " << static_cast<CbmRawEvent*>(fEvents->At(fEvents->GetEntriesFast() - 1))->ToString();
+ }
+ /// Clear event vector after usage
+ fpAlgo->ClearEventVector();
+}
+
+void CbmEvbuildRawTask::SaveHistos()
+{
+ if (fbWriteHistosToFairSink) {
+ if (!FairRootManager::Instance() || !FairRootManager::Instance()->GetSink()) {
+ LOG(error) << "No sink found";
+ return;
+ }
+ FairSink* sink = FairRootManager::Instance()->GetSink();
+ sink->WriteObject(dynamic_cast<TObject*>(fpAlgo->GetOutFolder()), nullptr);
+ }
+ else {
+
+ /// Obtain vector of pointers on each histo from the algo (+ optionally desired folder)
+ std::vector<std::pair<TNamed*, std::string>> vHistos = fpAlgo->GetHistoVector();
+
+ /// (Re-)Create ROOT file to store the histos
+ TDirectory* oldDir = NULL;
+ TFile* histoFile = NULL;
+ /// Store current directory position to allow restore later
+ oldDir = gDirectory;
+ /// open separate histo file in recreate mode
+ histoFile = new TFile(fsOutFileName, "RECREATE");
+ histoFile->cd();
+
+ /// Save all plots and create folders if needed
+ for (UInt_t uHisto = 0; uHisto < vHistos.size(); ++uHisto) {
+ /// Make sure we end up in chosen folder
+ const TString sFolder = vHistos[uHisto].second.data();
+ if (nullptr == gDirectory->Get(sFolder)) gDirectory->mkdir(sFolder);
+ gDirectory->cd(sFolder);
+
+ /// Write plot
+ vHistos[uHisto].first->Write();
+ histoFile->cd();
+ }
+
+ /// Restore original directory position
+ oldDir->cd();
+ histoFile->Close();
+ }
+}
+
+
+ClassImp(CbmEvbuildRawTask)
diff --git a/reco/eventbuilder/digis/CbmEvbuildRawTask.h b/reco/eventbuilder/digis/CbmEvbuildRawTask.h
new file mode 100644
index 0000000000000000000000000000000000000000..c2a3c29b94e4560a3de4e7f00cdb31be3e4e7474
--- /dev/null
+++ b/reco/eventbuilder/digis/CbmEvbuildRawTask.h
@@ -0,0 +1,177 @@
+/********************************************************************************
+ * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
+ * *
+ * This software is distributed under the terms of the *
+ * GNU Lesser General Public Licence (LGPL) version 3, *
+ * copied verbatim in the file "LICENSE" *
+ ********************************************************************************/
+#ifndef CBMEVBUILDRAWTASK_H
+#define CBMEVBUILDRAWTASK_H
+
+/// FAIRROOT headers
+#include "FairTask.h"
+
+/// FAIRSOFT headers (geant, boost, ...)
+
+/// C/C++ headers
+#include "CbmEvbuildRawAlgo.h"
+#include "CbmMuchBeamTimeDigi.h"
+#include "CbmMuchDigi.h"
+#include "CbmPsdDigi.h"
+#include "CbmRichDigi.h"
+#include "CbmStsDigi.h"
+#include "CbmTofDigi.h"
+#include "CbmTrdDigi.h"
+
+#include <tuple>
+
+#include <array>
+#include <map>
+#include <set>
+#include <vector>
+
+class CbmDigiManager;
+class RawEventBuilderDetector;
+class TClonesArray;
+class TStopwatch;
+
+enum class EOverlapModeRaw;
+
+class CbmEvbuildRawTask : public FairTask {
+public:
+ /** Default constructor **/
+ CbmEvbuildRawTask();
+
+ CbmEvbuildRawTask(const CbmEvbuildRawTask&) = delete;
+ CbmEvbuildRawTask operator=(const CbmEvbuildRawTask&) = delete;
+
+ /** Constructor with parameters (Optional) **/
+ // CbmEvbuildRawTask(Int_t verbose);
+
+ /** Destructor **/
+ ~CbmEvbuildRawTask();
+
+ /** Initiliazation of task at the beginning of a run **/
+ virtual InitStatus Init();
+
+ /** ReInitiliazation of task when the runID changes **/
+ virtual InitStatus ReInit();
+
+ /** Executed for each event. **/
+ virtual void Exec(Option_t*);
+
+ /** Load the parameter container from the runtime database **/
+ virtual void SetParContainers();
+
+ /** Finish task called at the end of the run **/
+ virtual void Finish();
+
+ /** Setters **/
+ void SetOutFilename(TString sNameIn) { fsOutFileName = sNameIn; }
+ void SetWriteHistosToFairSink(Bool_t var) { fbWriteHistosToFairSink = var; }
+
+ void SetFillHistos(Bool_t bFlag = kTRUE)
+ {
+ fbFillHistos = bFlag;
+ if (nullptr != fpAlgo) fpAlgo->SetFillHistos(fbFillHistos);
+ }
+ void SetReferenceDetector(RawEventBuilderDetector refDet)
+ {
+ if (nullptr != fpAlgo) fpAlgo->SetReferenceDetector(refDet);
+ }
+ void AddDetector(RawEventBuilderDetector selDet)
+ {
+ if (nullptr != fpAlgo) fpAlgo->AddDetector(selDet);
+ }
+ void RemoveDetector(RawEventBuilderDetector selDet)
+ {
+ if (nullptr != fpAlgo) fpAlgo->RemoveDetector(selDet);
+ }
+ void SetTriggerMinNumber(ECbmModuleId selDet, UInt_t uVal)
+ {
+ if (nullptr != fpAlgo) fpAlgo->SetTriggerMinNumber(selDet, uVal);
+ }
+ void SetTriggerMaxNumber(ECbmModuleId selDet, Int_t iVal)
+ {
+ if (nullptr != fpAlgo) fpAlgo->SetTriggerMaxNumber(selDet, iVal);
+ }
+ void SetTriggerWindow(ECbmModuleId det, Double_t dWinBeg, Double_t dWinEnd)
+ {
+ if (nullptr != fpAlgo) fpAlgo->SetTriggerWindow(det, dWinBeg, dWinEnd);
+ }
+ void SetTsParameters(Double_t dTsStartTime, Double_t dTsLength, Double_t dTsOverLength)
+ {
+ if (nullptr != fpAlgo) fpAlgo->SetTsParameters(dTsStartTime, dTsLength, dTsOverLength);
+ }
+ void SetEventOverlapMode(EOverlapModeRaw mode)
+ {
+ if (nullptr != fpAlgo) fpAlgo->SetEventOverlapMode(mode);
+ }
+ void SetIgnoreTsOverlap(Bool_t bFlagIn)
+ {
+ if (nullptr != fpAlgo) fpAlgo->SetIgnoreTsOverlap(bFlagIn);
+ }
+ void ChangeMuchBeamtimeDigiFlag(Bool_t bFlagIn = kFALSE)
+ {
+ if (nullptr != fpAlgo) fpAlgo->ChangeMuchBeamtimeDigiFlag(bFlagIn);
+ fbUseMuchBeamtimeDigi = bFlagIn;
+ }
+ void SetTimings(Bool_t bFlagIn = kTRUE)
+ {
+ if (nullptr != fpAlgo) fpAlgo->SetTimings(bFlagIn);
+ fbGetTimings = bFlagIn;
+ }
+
+ void PrintTimings();
+ void SetSeedTimeFiller(RawEventBuilderDetector seedDet);
+ void AddSeedTimeFillerToList(RawEventBuilderDetector seedDet);
+ void DumpSeedTimesFromDetList();
+ void SetSeedTimeWindow(Double_t beg, Double_t end) { fpAlgo->SetSeedTimeWindow(beg, end); }
+
+private:
+ void FillOutput();
+ void SaveHistos();
+
+ Bool_t fbUseMuchBeamtimeDigi = kTRUE; //! Switch between MUCH digi classes
+
+ CbmDigiManager* fDigiMan = nullptr;
+ const std::vector<CbmTofDigi>* fT0Digis = nullptr;
+ std::vector<CbmMuchDigi>* fMuchDigis = nullptr;
+ std::vector<CbmMuchBeamTimeDigi>* fMuchBeamTimeDigis = nullptr;
+ std::vector<CbmStsDigi>* fStsDigis = nullptr;
+ std::vector<CbmTrdDigi>* fTrdDigis = nullptr;
+ std::vector<CbmTofDigi>* fTofDigis = nullptr;
+ std::vector<CbmRichDigi>* fRichDigis = nullptr;
+ std::vector<CbmPsdDigi>* fPsdDigis = nullptr;
+ std::vector<Double_t>* fSeedTimes = nullptr;
+
+ std::vector<RawEventBuilderDetector> fSeedTimeDetList; //if multiple are desired
+ RawEventBuilderDetector fSeedTimeDet = kRawEventBuilderDetUndef; //single seed det
+
+ Double_t GetDigiTime(ECbmModuleId _system, UInt_t _entry);
+ UInt_t GetNofDigis(ECbmModuleId _system);
+
+ void FillSeedTimesFromDetList();
+
+ TStopwatch* fTimer = nullptr; //! is created when fbGetTimings is set before init
+ TStopwatch* fCopyTimer = nullptr; //! timing only for filling of std::vector<Digi> fields
+
+ CbmEvbuildRawAlgo* fpAlgo = nullptr;
+
+ TClonesArray* fEvents = nullptr; //! output container of CbmEvents
+
+ Bool_t fbFillHistos {kTRUE}; //! Switch ON/OFF filling of histograms
+ Bool_t fbWriteHistosToFairSink {kTRUE}; //! Write histos to FairRootManager instead of separate file
+ Bool_t fbGetTimings = kFALSE; //! Measure CPU time using stopwatch
+
+ /** Name of the histogram output file **/
+ TString fsOutFileName {"data/HistosEvtWin.root"};
+
+ Int_t fNofTs = 0;
+ Long64_t fNofEvents = 0;
+ Double_t fTime = 0.;
+
+ ClassDef(CbmEvbuildRawTask, 1);
+};
+
+#endif // CBMEVBUILDRAWTASK_H