t0 reconstruction

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;

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 {

reco/global/CMakeLists.txt

@@ -17,7 +17,8 @@ set(SRCS
   CbmFitGlobalTracksQa.cxx 
   CbmGlobalTrackFitterIdeal.cxx 
   CbmPVFinderIdeal.cxx       
-  CbmTrackMergerIdeal.cxx    
+  CbmTrackMergerIdeal.cxx  
+  CbmRecoTzero.cxx  
   #CbmTofMergerIdeal.cxx 
   )
 

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

reco/global/CbmRecoTzero.cxx

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

reco/global/CbmRecoTzero.h

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

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;