From d25b4751de25f2cc05194e4c0b037b67e7cea123 Mon Sep 17 00:00:00 2001
From: Volker Friese <v.friese@gsi.de>
Date: Fri, 11 Nov 2022 11:41:59 +0000
Subject: [PATCH] t0 reconstruction
---
core/data/CbmEvent.h | 26 ++++-
macro/run/run_reco.C | 8 ++
reco/global/CMakeLists.txt | 3 +-
reco/global/CbmGlobalLinkDef.h | 1 +
reco/global/CbmRecoTzero.cxx | 162 ++++++++++++++++++++++++++
reco/global/CbmRecoTzero.h | 90 ++++++++++++++
sim/response/base/CbmDigitization.cxx | 2 +-
7 files changed, 284 insertions(+), 8 deletions(-)
create mode 100644 reco/global/CbmRecoTzero.cxx
create mode 100644 reco/global/CbmRecoTzero.h
diff --git a/core/data/CbmEvent.h b/core/data/CbmEvent.h
index c9032922b8..e71c4945cf 100644
--- a/core/data/CbmEvent.h
+++ b/core/data/CbmEvent.h
@@ -137,6 +137,13 @@ public:
**/
double GetStartTime() const { return fTimeStart; }
+
+ /** Get t0
+ ** @value Reconstructed t0 [ns]
+ **/
+ double GetTzero() const { return fTzero; }
+
+
/** Set event number
** @value Event number
**/
@@ -160,6 +167,12 @@ public:
void SetStartTime(double startTime) { fTimeStart = startTime; }
+ /** Set t0
+ ** @param tZero T0 measurement [ns]
+ **/
+ void SetTzero(double tZero) { fTzero = tZero; }
+
+
/** Set the STS track index array
** @brief Sets the index array for STS tracks.
** Old content will be overwritten.
@@ -204,12 +217,13 @@ public:
private:
/** Event meta data **/
- int32_t fNumber; ///< Event number
- double fTimeStart; ///< Event start time [ns]
- double fTimeEnd; ///< Event end time [ns]
- int32_t fNofData; ///< Number of data objects in the event
- CbmVertex fVertex; ///< Primary Vertex
- CbmMatch* fMatch; ///< Match object to MCEvent
+ int32_t fNumber = -1; ///< Event number
+ double fTimeStart = 0.; ///< Event start time [ns]
+ double fTimeEnd = 0.; ///< Event end time [ns]
+ double fTzero = -999999.; ///< T0 of event for TOF PID [ns]
+ int32_t fNofData = 0; ///< Number of data objects in the event
+ CbmVertex fVertex = {}; ///< Primary Vertex
+ CbmMatch* fMatch = nullptr; ///< Match object to MCEvent
/** Arrays of indices to data types **/
std::map<ECbmDataType, std::vector<uint32_t>> fIndexMap;
diff --git a/macro/run/run_reco.C b/macro/run/run_reco.C
index e49474e04b..bbb6675c92 100644
--- a/macro/run/run_reco.C
+++ b/macro/run/run_reco.C
@@ -29,6 +29,7 @@
#include "CbmPrimaryVertexFinder.h"
#include "CbmPsdHitProducer.h"
#include "CbmRecoSts.h"
+#include "CbmRecoTzero.h"
#include "CbmRichHitProducer.h"
#include "CbmRichReconstruction.h"
#include "CbmSetup.h"
@@ -459,6 +460,13 @@ void run_reco(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice =
}
// ----------------------------------------------------------------------
+
+ // ----- T0 reconstruction --------------------------------------------
+ CbmRecoTzero* recoT0 = new CbmRecoTzero();
+ run->AddTask(recoT0);
+ std::cout << "-I- : Added task " << recoT0->GetName() << std::endl;
+ // ----------------------------------------------------------------------
+
} //? event-based reco
else {
diff --git a/reco/global/CMakeLists.txt b/reco/global/CMakeLists.txt
index 9c39ec8611..f683ed80d9 100644
--- a/reco/global/CMakeLists.txt
+++ b/reco/global/CMakeLists.txt
@@ -17,7 +17,8 @@ set(SRCS
CbmFitGlobalTracksQa.cxx
CbmGlobalTrackFitterIdeal.cxx
CbmPVFinderIdeal.cxx
- CbmTrackMergerIdeal.cxx
+ CbmTrackMergerIdeal.cxx
+ CbmRecoTzero.cxx
#CbmTofMergerIdeal.cxx
)
diff --git a/reco/global/CbmGlobalLinkDef.h b/reco/global/CbmGlobalLinkDef.h
index 06ee7981de..756cd9d997 100644
--- a/reco/global/CbmGlobalLinkDef.h
+++ b/reco/global/CbmGlobalLinkDef.h
@@ -19,6 +19,7 @@
#pragma link C++ class CbmGlobalTrackFitterIdeal + ;
#pragma link C++ class CbmPVFinderIdeal + ;
#pragma link C++ class CbmTrackMergerIdeal + ;
+#pragma link C++ class CbmRecoTzero + ;
//#pragma link C++ class CbmTofMergerIdeal+;
diff --git a/reco/global/CbmRecoTzero.cxx b/reco/global/CbmRecoTzero.cxx
new file mode 100644
index 0000000000..2ce3b5bed3
--- /dev/null
+++ b/reco/global/CbmRecoTzero.cxx
@@ -0,0 +1,162 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Volker Friese [committer] */
+
+
+#include "CbmRecoTzero.h"
+
+#include "CbmEvent.h"
+#include "CbmTzdDigi.h"
+
+#include <TClonesArray.h>
+#include <TStopwatch.h>
+
+#include <cassert>
+#include <iomanip>
+#include <sstream>
+
+
+using std::fixed;
+using std::left;
+using std::right;
+using std::setprecision;
+using std::setw;
+
+
+// ----- Constructor ---------------------------------------------------
+CbmRecoTzero::CbmRecoTzero(const char* name) : FairTask(name) {}
+// -------------------------------------------------------------------------
+
+
+// ----- Destructor ----------------------------------------------------
+CbmRecoTzero::~CbmRecoTzero() {}
+// -------------------------------------------------------------------------
+
+
+// ----- Initialization ------------------------------------------------
+InitStatus CbmRecoTzero::Init()
+{
+
+ std::cout << std::endl;
+ LOG(info) << "==========================================================";
+ LOG(info) << GetName() << ": Initialisation";
+
+ // --- Get FairRootManager
+ FairRootManager* ioman = FairRootManager::Instance();
+ assert(ioman);
+
+ // --- Get TzdDigi array
+ fTzdDigis = ioman->InitObjectAs<const std::vector<CbmTzdDigi>*>("TzdDigi");
+ if (!fTzdDigis) {
+ LOG(error) << GetName() << ": No TzdDigi array!";
+ return kERROR;
+ }
+ LOG(info) << "--- Found branch TzdDigi";
+
+ // --- Get CbmEvent array
+ fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("Event"));
+ if (fEvents) { LOG(info) << "--- Found branch Event"; }
+ else {
+ fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
+ if (!fEvents) {
+ LOG(error) << GetName() << ": No Event nor CbmEvent branch! Task will be inactive.";
+ return kERROR;
+ }
+ LOG(info) << "--- Found branch CbmEvent";
+ }
+
+ LOG(info) << GetName() << ": Initialisation successful";
+ LOG(info) << "==========================================================";
+ std::cout << std::endl;
+ return kSUCCESS;
+}
+// -------------------------------------------------------------------------
+
+
+// ----- Public method Exec --------------------------------------------
+void CbmRecoTzero::Exec(Option_t*)
+{
+
+ // Timer
+ TStopwatch timer;
+ timer.Start();
+ CbmRecoTzeroMoniData tsMonitor {};
+
+ // Event loop
+ Int_t nEvents = fEvents->GetEntriesFast();
+ for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
+ CbmEvent* event = dynamic_cast<CbmEvent*>(fEvents->At(iEvent));
+ assert(event);
+ Int_t nDigis = fTzdDigis->size();
+ double tzero = -999999.;
+ switch (nDigis) {
+
+ // If there is no TZD digi, set t0 to -999999 (error code).
+ case 0: {
+ tzero = -999999.;
+ tsMonitor.fNumEvtsTzd0++;
+ break;
+ }
+
+ // If there is exactly one TZD digi, take the event time from there
+ case 1: {
+ tzero = fTzdDigis->at(0).GetTime();
+ tsMonitor.fNumEvtsTzd1++;
+ break;
+ }
+
+ // If there are more than one TZD digis, set t0 to -999999 (error code).
+ default: {
+ tzero = -999999.;
+ tsMonitor.fNumEvtsTzdn++;
+ break;
+ }
+ }
+
+ event->SetTzero(tzero);
+ tsMonitor.fNumEvents++;
+ }
+
+
+ // Timeslice monitor
+ timer.Stop();
+ tsMonitor.fExecTime = 1000. * timer.RealTime();
+ tsMonitor.fNumTs = 1;
+ std::stringstream logOut;
+ logOut << setw(20) << left << GetName() << " [";
+ logOut << fixed << setw(8) << setprecision(1) << right << timer.RealTime() * 1000. << " ms] ";
+ logOut << "TS " << fMonitor.fNumTs << ", events " << tsMonitor.fNumEvents;
+ logOut << " (1 TZD: " << tsMonitor.fNumEvtsTzd1;
+ logOut << " , 0 TZD: " << tsMonitor.fNumEvtsTzd0;
+ logOut << " , n TZD: " << tsMonitor.fNumEvtsTzdn << ")";
+ LOG(info) << logOut.str();
+
+ // Run monitor
+ fMonitor += tsMonitor;
+}
+// -------------------------------------------------------------------------
+
+
+// ----- Public method Finish ------------------------------------------
+void CbmRecoTzero::Finish()
+{
+ double tExec = fMonitor.fExecTime / double(fMonitor.fNumTs);
+ double evtsPerTs = double(fMonitor.fNumEvents) / double(fMonitor.fNumTs);
+ double fracTzd1 = 100. * double(fMonitor.fNumEvtsTzd1) / double(fMonitor.fNumEvents);
+ double fracTzd0 = 100. * double(fMonitor.fNumEvtsTzd0) / double(fMonitor.fNumEvents);
+ double fracTzdn = 100. * double(fMonitor.fNumEvtsTzdn) / double(fMonitor.fNumEvents);
+ std::cout << std::endl;
+ LOG(info) << "=====================================";
+ LOG(info) << GetName() << ": Run summary";
+ LOG(info) << "Time slices : " << fMonitor.fNumTs;
+ LOG(info) << "Exec time / TS : " << fixed << setprecision(2) << tExec << " ms";
+ LOG(info) << "Events / TS : " << fixed << setprecision(2) << evtsPerTs;
+ LOG(info) << "Fraction with 1 TZD : " << fixed << setprecision(2) << fracTzd1 << " %";
+ LOG(info) << "Fraction with 0 TZD : " << fixed << setprecision(2) << fracTzd0 << " %";
+ LOG(info) << "Fraction with n TZD : " << fixed << setprecision(2) << fracTzdn << " %";
+ LOG(info) << "=====================================";
+}
+// -------------------------------------------------------------------------
+
+
+ClassImp(CbmRecoTzero)
diff --git a/reco/global/CbmRecoTzero.h b/reco/global/CbmRecoTzero.h
new file mode 100644
index 0000000000..3a58bca296
--- /dev/null
+++ b/reco/global/CbmRecoTzero.h
@@ -0,0 +1,90 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Volker Friese [committer] */
+
+
+#ifndef CBMRECOTZERO_H
+#define CBMRECOTZERO_H 1
+
+#include "CbmTzdDigi.h"
+
+#include "FairTask.h"
+
+#include <vector>
+
+class TClonesArray;
+
+
+/** @struct CbmRecoTzeroMoniData
+ ** @brief Monitor data for T0 reconstruction
+ ** @author Volker Friese <v.friese@gsi.de>
+ ** @since 10.11.2022
+ **/
+struct CbmRecoTzeroMoniData {
+ size_t fNumTs = 0;
+ size_t fNumEvents = 0;
+ size_t fNumEvtsTzd0 = 0;
+ size_t fNumEvtsTzd1 = 0;
+ size_t fNumEvtsTzdn = 0;
+ double fExecTime = 0.;
+
+ CbmRecoTzeroMoniData& operator+=(const CbmRecoTzeroMoniData& other)
+ {
+ fNumTs += other.fNumTs;
+ fNumEvents += other.fNumEvents;
+ fNumEvtsTzd0 += other.fNumEvtsTzd0;
+ fNumEvtsTzd1 += other.fNumEvtsTzd1;
+ fNumEvtsTzdn += other.fNumEvtsTzdn;
+ fExecTime += other.fExecTime;
+ return *this;
+ }
+};
+
+/** @class CbmRecoTzero
+ ** @brief Task class for reconstruction of the event t0
+ ** @author Volker Friese <v.friese@gsi.de>
+ ** @since 10.11.2022
+ **
+ ** The current implementation reads the t0 information from the TdzDigi object. t0 is set to -1. if
+ ** no such object is in the event, and to -2. if there are several.
+ **/
+class CbmRecoTzero : public FairTask {
+
+public:
+ /** @brief Constructor
+ **
+ ** @param name Name of task
+ ** @param title Title of task
+ **/
+ CbmRecoTzero(const char* name = "RecoTzero");
+
+
+ /** @brief Destructor **/
+ virtual ~CbmRecoTzero();
+
+
+ /** @brief Initialisation **/
+ virtual InitStatus Init();
+
+
+ /** @brief Task execution **/
+ virtual void Exec(Option_t* opt);
+
+
+ /** @brief End-of-timeslice action **/
+ virtual void Finish();
+
+
+private:
+ // --- Data
+ const std::vector<CbmTzdDigi>* fTzdDigis = nullptr; ///< TZD digis
+ TClonesArray* fEvents = nullptr; ///< CbmEvent
+
+ // --- Monitor
+ CbmRecoTzeroMoniData fMonitor = {}; ///< Monitor data
+
+
+ ClassDef(CbmRecoTzero, 1);
+};
+
+#endif
diff --git a/sim/response/base/CbmDigitization.cxx b/sim/response/base/CbmDigitization.cxx
index 5109b7f44c..30ab4cb3d3 100644
--- a/sim/response/base/CbmDigitization.cxx
+++ b/sim/response/base/CbmDigitization.cxx
@@ -205,7 +205,7 @@ Int_t CbmDigitization::CreateDefaultDigitizers()
break;
case ECbmModuleId::kT0:
fDigitizers[system]->SetDigitizer(new CbmTzdDigitize());
- ss << "BMON ";
+ ss << "TZD ";
nDigis++;
break;
default: LOG(fatal) << fName << ": Unknown system " << system; break;
--
GitLab