diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt
index 5f807e9cbf5b63df71177908170fb289c2a18aea..6e80b666e2d6cb2d0f2b5ad9c39529ba29ca2aa8 100644
--- a/algo/CMakeLists.txt
+++ b/algo/CMakeLists.txt
@@ -10,6 +10,7 @@ set(SRCS
detectors/sts/UnpackSts.cxx
detectors/much/MuchReadoutConfig.cxx
detectors/much/UnpackMuch.cxx
+ detectors/tof/HitFinderTof.cxx
)
add_library(Algo SHARED ${SRCS})
@@ -24,7 +25,7 @@ target_include_directories(Algo
${CMAKE_CURRENT_SOURCE_DIR}/detectors/tof
)
-target_link_libraries(Algo PUBLIC OnlineData INTERFACE FairLogger::FairLogger external::fles_ipc)
+target_link_libraries(Algo PUBLIC OnlineData ROOT::GenVector INTERFACE FairLogger::FairLogger external::fles_ipc)
target_compile_definitions(Algo PUBLIC NO_ROOT)
install(TARGETS Algo DESTINATION lib)
diff --git a/algo/detectors/tof/HitFinderTof.cxx b/algo/detectors/tof/HitFinderTof.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..dfa2b493710280223d8d2c0dc1dda9d052ad8674
--- /dev/null
+++ b/algo/detectors/tof/HitFinderTof.cxx
@@ -0,0 +1,186 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Pierre-Alain Loizeau */
+
+#include "HitFinderTof.h"
+
+// TOF Classes and includes
+#include "CbmTofDigi.h"
+
+// C++ Classes and includes
+#include <iomanip>
+#include <iostream>
+
+namespace cbm::algo
+{
+ HitFinderTof::resultType HitFinderTof::operator()(std::vector<CbmTofDigi> digisIn,
+ const std::vector<int32_t>& digiIndexIn)
+ {
+ // Intermediate storage variables
+ std::vector<std::vector<CbmTofDigi*>> digisExp; //[nbCh][nDigis]
+ std::vector<std::vector<int32_t>> digisInd; //[nbCh][nDigis]
+
+ digisExp.resize(fParams.fChanPar.size());
+ digisInd.resize(fParams.fChanPar.size());
+
+ { // digi calibration
+ const int32_t numClWalkBinX = fParams.numClWalkBinX;
+ const double TOTMax = fParams.TOTMax;
+ const double TOTMin = fParams.TOTMin;
+
+ for (size_t iDigi = 0; iDigi < digisIn.size(); iDigi++) {
+ CbmTofDigi* pDigi = &digisIn[iDigi];
+ assert(pDigi);
+
+ // These are doubles in the digi class
+ const double chan = pDigi->GetChannel();
+ const double side = pDigi->GetSide();
+ const double charge = pDigi->GetTot();
+ const double time = pDigi->GetTime();
+
+ digisExp[chan].push_back(pDigi);
+ digisInd[chan].push_back(digiIndexIn[iDigi]);
+
+ // calibrate Digi Time
+ pDigi->SetTime(time - fParams.fChanPar[chan].fvCPTOff[side]);
+
+ // calibrate Digi ToT
+ pDigi->SetTot(charge * fParams.fChanPar[chan].fvCPTotGain[side]);
+
+ { // walk correction
+ const double totBinSize = (TOTMax - TOTMin) / 2. / numClWalkBinX;
+ int32_t iWx = (int32_t)((charge - TOTMin / 2.) / totBinSize);
+ if (0 > iWx) iWx = 0;
+ if (iWx > (numClWalkBinX - 1)) iWx = numClWalkBinX - 1;
+
+ std::vector<double>& cpWalk = fParams.fChanPar[chan].fvCPWalk[side];
+ double dWT = cpWalk[iWx];
+ const double dDTot = (charge - TOTMin / 2.) / totBinSize - (double) iWx - 0.5;
+
+ if (dDTot > 0) {
+ if (iWx < numClWalkBinX - 1) { // linear interpolation to next bin
+ dWT += dDTot * (cpWalk[iWx + 1] - cpWalk[iWx]);
+ }
+ }
+ else {
+ if (0 < iWx) { // linear interpolation to next bin
+ dWT -= dDTot * (cpWalk[iWx - 1] - cpWalk[iWx]);
+ }
+ }
+ pDigi->SetTime(time - dWT);
+ }
+ }
+ } // digi calibration
+
+ // Hit variables
+ TofCluster cluster;
+ std::vector<TofCluster> clustersOut;
+
+ // Reference cell of a cluster
+ TofCell* trafoCell = nullptr;
+ int32_t iTrafoCell = -1;
+
+ // Last Channel Temp variables
+ int32_t lastChan = -1;
+ double lastPosY = 0.0;
+ double lastTime = 0.0;
+
+ //reset counter
+ numSameSide = 0;
+
+ cluster.reset(); // Don't spread clusters over RPCs
+
+ for (int32_t chan = 0; chan < fParams.fChanPar.size(); chan++) {
+
+ std::vector<CbmTofDigi*>& storDigiExp = digisExp[chan];
+ std::vector<int32_t>& storDigiInd = digisInd[chan];
+ HitFinderTofChanPar& chanPar = fParams.fChanPar[chan];
+
+ while (1 < storDigiExp.size()) {
+ while ((storDigiExp[0])->GetSide() == (storDigiExp[1])->GetSide()) {
+ // Not one Digi of each end!
+ numSameSide++;
+ storDigiExp.erase(storDigiExp.begin());
+ storDigiInd.erase(storDigiInd.begin());
+ if (2 > storDigiExp.size()) break;
+ }
+ if (2 > storDigiExp.size()) break; // 2 Digis = both sides present
+
+ TofCell* channelInfo = &chanPar.cell;
+ CbmTofDigi* xDigiA = storDigiExp[0];
+ CbmTofDigi* xDigiB = storDigiExp[1];
+
+ // The "Strip" time is the mean time between each end
+ const double time = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
+
+ // Weight is the total charge => sum of both ends ToT
+ const double totSum = xDigiA->GetTot() + xDigiB->GetTot();
+
+ if (nullptr == trafoCell) {
+ trafoCell = channelInfo;
+ iTrafoCell = chan;
+ }
+
+ // switch to local coordinates, (0,0,0) is in the center of counter ?
+ // It is assumed that the cell size is the same for all channels within one rpc!
+ ROOT::Math::XYZVector pos(((double) (-iTrafoCell + chan)) * trafoCell->sizeX, 0., 0.);
+
+ if (1 == xDigiA->GetSide()) {
+ // 0 is the top side, 1 is the bottom side
+ pos.SetY(fParams.CPSigPropSpeed * (xDigiA->GetTime() - xDigiB->GetTime()) / 2.0);
+ }
+ else {
+ // 0 is the bottom side, 1 is the top side
+ pos.SetY(fParams.CPSigPropSpeed * (xDigiB->GetTime() - xDigiA->GetTime()) / 2.0);
+ }
+
+ if (channelInfo->sizeY / 2.0 < pos.Y() || -1 * channelInfo->sizeY / 2.0 > pos.Y()) {
+ // if outside of strip limits, the pair is bad => try to remove one end and check the next pair
+ // (if another possibility exist)
+ storDigiExp.erase(storDigiExp.begin());
+ storDigiInd.erase(storDigiInd.begin());
+ continue;
+ } // Pair leads to hit oustide of strip limits
+
+ // Now check if a hit/cluster is already started
+ if (0 < cluster.numChan()) {
+ // a cluster is already started => check distance in space/time
+ // For simplicity, just check along strip direction for now
+ // and break cluster when a not fired strip is found
+ if (!(std::abs(time - lastTime) < fParams.maxTimeDist && lastChan == chan - 1
+ && std::abs(pos.Y() - lastPosY) < fParams.maxSpaceDist)) {
+ // simple error scaling with TOT
+ // weightedTimeErr = TMath::Sqrt( weightedTimeErr ) * SysTimeRes / weightsSum;
+ // OR harcoded value: weightedTimeErr = SysTimeRes;
+ cluster.normalize(fParams.SysTimeRes);
+
+ // Save Hit and start a new one
+ cluster.finalize(*trafoCell, iTrafoCell, fParams);
+ clustersOut.push_back(cluster);
+ cluster.reset();
+ }
+ }
+ cluster.add(pos, time, totSum, totSum, storDigiInd[0], storDigiInd[1]);
+ storDigiExp.erase(storDigiExp.begin());
+ storDigiExp.erase(storDigiExp.begin());
+ storDigiInd.erase(storDigiInd.begin());
+ storDigiInd.erase(storDigiInd.begin());
+
+ lastChan = chan;
+ lastPosY = pos.Y();
+ lastTime = time;
+ } // while (1 < storDigiExp.size())
+ storDigiExp.clear();
+ storDigiInd.clear();
+ } // for( int32_t chan = 0; chan < iNbCh; chan++ )
+
+ // Save last cluster if it exists
+ if (0 < cluster.numChan()) {
+ cluster.normalize(fParams.SysTimeRes);
+ cluster.finalize(*trafoCell, iTrafoCell, fParams);
+ clustersOut.push_back(cluster);
+ }
+
+ return clustersOut;
+ }
+} /* namespace cbm::algo */
diff --git a/algo/detectors/tof/HitFinderTof.h b/algo/detectors/tof/HitFinderTof.h
new file mode 100644
index 0000000000000000000000000000000000000000..cb21482196601d10006cd63a76fdf36ac38024bc
--- /dev/null
+++ b/algo/detectors/tof/HitFinderTof.h
@@ -0,0 +1,164 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Pierre-Alain Loizeau */
+
+#ifndef HITFINDERTOF_H
+#define HITFINDERTOF_H
+
+// TOF Classes and includes
+class CbmTofDigi;
+
+// ROOT Classes and includes
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+// C++ Classes and includes
+#include <vector>
+
+#include <cmath>
+
+namespace cbm::algo
+{
+ struct TofCell {
+ double sizeX, sizeY;
+ ROOT::Math::XYZVector pos;
+ ROOT::Math::Rotation3D rotation;
+ };
+
+ struct HitFinderTofChanPar {
+ std::vector<double> fvCPTOff; //[nbSide]
+ std::vector<double> fvCPTotGain; //[nbSide]
+ std::vector<std::vector<double>> fvCPWalk; //[nbSide][nbWalkBins]
+ int32_t address; //unique address
+ TofCell cell;
+ };
+
+ struct HitFinderTofRpcPar {
+ double FeeTimeRes;
+ double SysTimeRes;
+ double CPSigPropSpeed;
+ double outTimeFactor;
+ double TOTMax;
+ double TOTMin;
+ double maxTimeDist;
+ double maxSpaceDist;
+ double gapSize;
+ int32_t numGaps;
+ int32_t numClWalkBinX;
+ std::vector<HitFinderTofChanPar> fChanPar = {};
+ };
+
+ struct TofCluster {
+ //temporary values
+ std::vector<int32_t> vDigiIndRef;
+ ROOT::Math::XYZVector weightedPos = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
+ double weightedTime = 0.0;
+ double weightedTimeErr = 0.0;
+ double weightsSum = 0.0;
+
+ //after finalization
+ ROOT::Math::XYZVector globalPos = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
+ ROOT::Math::XYZVector globalErr = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
+ int32_t channel = 0;
+ int32_t detId = 0;
+
+ int32_t numChan() { return vDigiIndRef.size() / 2; }
+
+ void reset()
+ {
+ vDigiIndRef.clear();
+ weightedPos = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
+ globalPos = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
+ globalErr = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
+ weightedTime = 0.0;
+ weightedTimeErr = 0.0;
+ weightsSum = 0.0;
+ channel = 0;
+ detId = 0;
+ }
+
+ void add(ROOT::Math::XYZVector pos, double time, double timeErr, double weight, int32_t digiIndA, int32_t digiIndB)
+ {
+ vDigiIndRef.push_back(digiIndA);
+ vDigiIndRef.push_back(digiIndB);
+ weightedPos += pos * weight;
+ weightedTime += time * weight;
+ weightedTimeErr += timeErr * weight;
+ weightsSum += weight;
+ }
+
+ void normalize(double timeErr)
+ {
+ // a/=b is translated to a := a*(1/b) in the ROOT::Math::XYZVector class, which has a different
+ // rounding behavior than a := (a/b). In rare cases this leads to 1.000... becoming 0.999... inside
+ // the floor() operation in the finalize() function of this class, and results in a different
+ // channel being associated with the cluster. To reproduce the output of the old hit finder, we
+ // divide element-wise instead. Perhaps floor() should be replaced by round().
+ //////// weightedPos /= weightsSum;
+
+ weightedPos.SetXYZ(weightedPos.X() / weightsSum, weightedPos.Y() / weightsSum, weightedPos.Z() / weightsSum);
+ weightedTime /= weightsSum;
+ weightedTimeErr = timeErr;
+ }
+
+ void finalize(const TofCell& trafoCell, const int32_t iTrafoCell, const HitFinderTofRpcPar& par)
+ {
+ // prepare local->global trafo
+ ROOT::Math::Rotation3D rotMatrix = trafoCell.rotation;
+
+ // get offset from weighted cluster position by rotation to master frame
+ ROOT::Math::XYZVector hitposOffset = rotMatrix(weightedPos);
+
+ //get hit position by adding offset to cell coordinates
+ globalPos = trafoCell.pos + hitposOffset;
+
+ // Simple errors, not properly done at all for now
+ // Right way of doing it should take into account the weight distribution
+ // and real system time resolution
+ ROOT::Math::XYZVector hitErrLocal;
+ hitErrLocal.SetX(trafoCell.sizeX / std::sqrt(12.0)); // Single strips approximation
+ hitErrLocal.SetY(par.FeeTimeRes * par.CPSigPropSpeed); // Use the electronics resolution
+ hitErrLocal.SetZ(par.numGaps * par.gapSize / 10.0 // Change gap size in cm
+ / std::sqrt(12.0)); // Use full RPC thickness as "Channel" Z size
+
+ //store results
+ globalErr = rotMatrix(hitErrLocal);
+ channel = iTrafoCell + floor(weightedPos.X() / trafoCell.sizeX);
+ detId = par.fChanPar[channel].address;
+ weightedTime *= par.outTimeFactor;
+ }
+ };
+
+ class HitFinderTof {
+ public:
+ typedef std::vector<TofCluster> resultType;
+
+ /**
+ ** @brief Constructor.
+ **/
+ HitFinderTof() {};
+
+ /**
+ ** @brief Destructor.
+ **/
+ ~HitFinderTof() {};
+
+ /**
+ ** @brief Build clusters out of ToF Digis and store the resulting info in a TofHit.
+ **/
+ resultType operator()(std::vector<CbmTofDigi> digisIn, const std::vector<int32_t>& digiIndexIn);
+
+ /** @brief Set the parameter container
+ ** @param params Vectorer to parameter container
+ **/
+ void SetParams(std::unique_ptr<HitFinderTofRpcPar> params) { fParams = *(std::move(params)); }
+
+ private:
+ HitFinderTofRpcPar fParams = {}; ///< Parameter container
+
+ int32_t numSameSide; // Digis quality
+ };
+
+} /* namespace cbm::algo */
+
+#endif // HITFINDERTOF_H
diff --git a/reco/tasks/CMakeLists.txt b/reco/tasks/CMakeLists.txt
index 49ea82ea7382560ac1a17e08a032511a1bc2a6ab..6c7302a21cc694f3fe6713dff4e69a2af41d8f03 100644
--- a/reco/tasks/CMakeLists.txt
+++ b/reco/tasks/CMakeLists.txt
@@ -3,7 +3,6 @@
set(INCLUDE_DIRECTORIES
${CMAKE_CURRENT_SOURCE_DIR}
- ${CMAKE_SOURCE_DIR}/algo/evbuild
)
set(SRCS
@@ -14,6 +13,7 @@ set(SRCS
CbmTaskEventsCloneInToOut.cxx
CbmTaskMakeRecoEvents.cxx
CbmTaskTriggerDigi.cxx
+ CbmTaskTofHitFinder.cxx
CbmTaskUnpack.cxx
)
@@ -27,6 +27,7 @@ set(LIBRARY_NAME CbmRecoTasks)
set(LINKDEF ${LIBRARY_NAME}LinkDef.h)
set(PUBLIC_DEPENDENCIES
CbmData
+ CbmTofBase
FairRoot::Base
Algo
ROOT::Core
diff --git a/reco/tasks/CbmRecoTasksLinkDef.h b/reco/tasks/CbmRecoTasksLinkDef.h
index 32def2aa29305cd57e458a56ab5b3526ac2d3850..980fb920edc0e77f102f7470cd8bbe2ff346d3d9 100644
--- a/reco/tasks/CbmRecoTasksLinkDef.h
+++ b/reco/tasks/CbmRecoTasksLinkDef.h
@@ -18,6 +18,7 @@
#pragma link C++ class CbmTaskDigiEventQa + ;
#pragma link C++ class CbmTaskEventsCloneInToOut + ;
#pragma link C++ class CbmTaskMakeRecoEvents + ;
+#pragma link C++ class CbmTaskTofHitFinder + ;
#pragma link C++ class CbmTaskTriggerDigi + ;
#pragma link C++ class CbmTaskUnpack + ;
diff --git a/reco/tasks/CbmTaskTofHitFinder.cxx b/reco/tasks/CbmTaskTofHitFinder.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..079c9ae224fe2a3e77738324f3f857bf796e5906
--- /dev/null
+++ b/reco/tasks/CbmTaskTofHitFinder.cxx
@@ -0,0 +1,472 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Pierre-Alain Loizeau, Norbert Herrmann */
+
+#include "CbmTaskTofHitFinder.h"
+
+// TOF Classes and includes
+#include "CbmTofAddress.h" // in cbmdata/tof
+#include "CbmTofCell.h" // in tof/TofData
+#include "CbmTofCreateDigiPar.h"
+#include "CbmTofDetectorId_v14a.h" // in cbmdata/tof
+#include "CbmTofDetectorId_v21a.h" // in cbmdata/tof
+#include "CbmTofDigi.h" // in cbmdata/tof
+#include "CbmTofDigiBdfPar.h" // in tof/TofParam
+#include "CbmTofDigiPar.h" // in tof/TofParam
+#include "CbmTofGeoHandler.h" // in tof/TofTools
+#include "CbmTofHit.h" // in cbmdata/tof
+
+// CBMroot classes and includes
+#include "CbmDigiManager.h"
+#include "CbmEvent.h"
+#include "CbmMatch.h"
+
+// FAIR classes and includes
+#include "FairEventHeader.h" // from CbmStsDigitize, for GetEventInfo
+#include "FairMCEventHeader.h" // from CbmStsDigitize, for GetEventInfo
+#include "FairRootManager.h"
+#include "FairRunAna.h"
+#include "FairRunSim.h" // from CbmStsDigitize, for GetEventInfo
+#include "FairRuntimeDb.h"
+#include <Logger.h>
+
+// ROOT Classes and includes
+#include "TClonesArray.h"
+
+// C++ Classes and includes
+#include <iomanip>
+#include <iostream>
+
+using std::fixed;
+using std::left;
+using std::pair;
+using std::right;
+using std::setprecision;
+using std::setw;
+using std::stringstream;
+
+using cbm::algo::HitFinderTofRpcPar;
+using cbm::algo::TofCell;
+using cbm::algo::TofCluster;
+
+const int32_t numClWalkBinX = 20;
+const double TTotMean = 2.E4;
+
+/************************************************************************************/
+CbmTaskTofHitFinder::CbmTaskTofHitFinder() : FairTask("CbmTaskTofHitFinder"), fGeoHandler(new CbmTofGeoHandler()) {}
+
+CbmTaskTofHitFinder::CbmTaskTofHitFinder(const char* name, int32_t verbose)
+ : FairTask(TString(name), verbose)
+ , fGeoHandler(new CbmTofGeoHandler())
+{
+}
+
+CbmTaskTofHitFinder::~CbmTaskTofHitFinder()
+{
+ if (fGeoHandler) delete fGeoHandler;
+}
+
+InitStatus CbmTaskTofHitFinder::Init()
+{
+ fDigiMan = CbmDigiManager::Instance(), fDigiMan->Init();
+ if (false == RegisterInputs()) return kFATAL;
+ if (false == RegisterOutputs()) return kFATAL;
+ if (false == InitParameters()) return kFATAL;
+ if (false == InitCalibParameter()) return kFATAL;
+ if (false == InitAlgos()) return kFATAL;
+ return kSUCCESS;
+}
+
+void CbmTaskTofHitFinder::SetParContainers()
+{
+ LOG(info) << " CbmTaskTofHitFinder => Get the digi parameters for tof";
+
+ // Get Base Container
+ FairRunAna* ana = FairRunAna::Instance();
+ FairRuntimeDb* rtdb = ana->GetRuntimeDb();
+
+ fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
+ LOG(info) << " CbmTaskTofHitFinder::SetParContainers found " << fDigiPar->GetNrOfModules() << " cells ";
+
+ fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
+}
+
+void CbmTaskTofHitFinder::Exec(Option_t* /*option*/)
+{
+ // Start timer counter
+ fTimer.Start();
+
+ fTofHitsColl->Clear("C");
+ fTofDigiMatchColl->Delete();
+
+ fStart.Set();
+
+ // --- Local variables
+ int32_t nDigisAll = CbmDigiManager::GetNofDigis(ECbmModuleId::kTof);
+ int32_t nEvents = 0;
+ int32_t nDigis = 0;
+ int32_t nHits = 0;
+ CbmEvent* event = nullptr;
+ pair<int32_t, int32_t> result;
+
+ // --- Time-slice mode: process entire digi array
+ if (!fEvents) {
+ result = BuildClusters(nullptr);
+ nDigis = result.first;
+ nHits = result.second;
+ }
+
+ // --- Event-based mode: read and process event after event
+ else {
+ nEvents = fEvents->GetEntriesFast();
+ for (int32_t iEvent = 0; iEvent < nEvents; iEvent++) {
+ event = dynamic_cast<CbmEvent*>(fEvents->At(iEvent));
+ assert(event);
+ result = BuildClusters(event);
+ nDigis += result.first;
+ nHits += result.second;
+ } //# events
+ } //? event mode
+
+ fStop.Set();
+ fTimer.Stop();
+
+ // --- Timeslice log
+ stringstream logOut;
+ logOut << setw(20) << left << GetName() << " [";
+ logOut << fixed << setw(8) << setprecision(1) << right << fTimer.RealTime() * 1000. << " ms] ";
+ logOut << "TS " << fiNofTs;
+ if (fEvents) logOut << ", events " << nEvents;
+ logOut << ", digis " << nDigis << " / " << nDigisAll;
+ logOut << ", hits " << nHits;
+ LOG(info) << logOut.str();
+
+ // --- Update Counters
+ fiNofTs++;
+ fiNofEvents += nEvents;
+ fNofDigisAll += nDigisAll;
+ fNofDigisUsed += nDigis;
+ fdNofHitsTot += nHits;
+ fdTimeTot += fTimer.RealTime();
+}
+
+void CbmTaskTofHitFinder::Finish()
+{
+ // Screen log
+ std::cout << std::endl;
+ LOG(info) << "=====================================";
+ LOG(info) << GetName() << ": Run summary";
+ LOG(info) << "Time slices : " << fiNofTs;
+ LOG(info) << "Digis / TS : " << fixed << setprecision(2) << fNofDigisAll / double(fiNofTs);
+ LOG(info) << "Hits / TS : " << fixed << setprecision(2) << fdNofHitsTot / double(fiNofTs);
+ LOG(info) << "Time / TS : " << fixed << setprecision(2) << 1000. * fdTimeTot / double(fiNofTs) << " ms";
+ if (fEvents) {
+ double unusedFrac = 100. * (1. - fNofDigisUsed / fNofDigisAll);
+ LOG(info) << "Digis outside events : " << fNofDigisAll - fNofDigisUsed << " = " << unusedFrac << " %";
+ LOG(info) << "Events : " << fiNofEvents;
+ LOG(info) << "Events / TS : " << fixed << setprecision(2) << double(fiNofEvents) / double(fiNofTs);
+ LOG(info) << "Digis / event : " << fixed << setprecision(2) << fNofDigisUsed / double(fiNofEvents);
+ LOG(info) << "Hits / event : " << fixed << setprecision(2) << fdNofHitsTot / double(fiNofEvents);
+ }
+ LOG(info) << "=====================================\n";
+}
+
+/************************************************************************************/
+// Functions common for all clusters approximations
+bool CbmTaskTofHitFinder::RegisterInputs()
+{
+ FairRootManager* fManager = FairRootManager::Instance();
+
+ // --- Check input branch (TofDigiExp). If not present, set task inactive.
+ if (!fDigiMan->IsPresent(ECbmModuleId::kTof)) {
+ LOG(error) << GetName() << ": No TofDigi input array present; "
+ << "task will be inactive.";
+ return kERROR;
+ }
+
+ // --- Look for event branch
+ fEvents = dynamic_cast<TClonesArray*>(fManager->GetObject("Event"));
+ if (fEvents) LOG(info) << GetName() << ": Found Event branch; run event-based";
+ else {
+ fEvents = dynamic_cast<TClonesArray*>(fManager->GetObject("CbmEvent"));
+ if (fEvents) LOG(info) << GetName() << ": Found CbmEvent branch; run event-based";
+ else
+ LOG(info) << GetName() << ": No event branch found; run time-based";
+ }
+ return true;
+}
+
+bool CbmTaskTofHitFinder::RegisterOutputs()
+{
+ FairRootManager* rootMgr = FairRootManager::Instance();
+ fTofHitsColl = new TClonesArray("CbmTofHit");
+
+ // Flag check to control whether digis are written in output root file
+ rootMgr->Register("TofHit", "Tof", fTofHitsColl, IsOutputBranchPersistent("TofHit"));
+
+ fTofDigiMatchColl = new TClonesArray("CbmMatch", 100);
+ rootMgr->Register("TofHitDigiMatch", "Tof", fTofDigiMatchColl, IsOutputBranchPersistent("TofHitDigiMatch"));
+
+ return true;
+}
+
+bool CbmTaskTofHitFinder::InitParameters()
+{
+ if (fDigiBdfPar->UseExpandedDigi() == false) {
+ LOG(fatal) << " Compressed Digis not implemented ... ";
+ return false;
+ }
+
+ // Initialize the TOF GeoHandler
+ bool isSimulation = false;
+ int32_t iGeoVersion = fGeoHandler->Init(isSimulation);
+ LOG(info) << "CbmTaskTofHitFinder::InitParameters with GeoVersion " << iGeoVersion;
+
+ if (k14a > iGeoVersion) {
+ LOG(error) << "CbmTaskTofHitFinder::InitParameters => Only compatible "
+ "with geometries after v12b !!!";
+ return false;
+ }
+
+ if (nullptr != fTofId)
+ LOG(info) << "CbmTaskTofHitFinder::InitParameters with GeoVersion " << fGeoHandler->GetGeoVersion();
+ else {
+ switch (iGeoVersion) {
+ case k14a: fTofId = new CbmTofDetectorId_v14a(); break;
+ case k21a: fTofId = new CbmTofDetectorId_v21a(); break;
+ default: LOG(error) << "CbmTaskTofHitFinder::InitParameters => Invalid geometry!!!" << iGeoVersion; return false;
+ }
+ }
+
+ LOG(info) << "=> Get the digi parameters for tof";
+ FairRunAna* ana = FairRunAna::Instance();
+ FairRuntimeDb* rtdb = ana->GetRuntimeDb();
+
+ // create digitization parameters from geometry file
+ CbmTofCreateDigiPar* tofDigiPar = new CbmTofCreateDigiPar("TOF Digi Producer", "TOF task");
+ LOG(info) << "Create DigiPar ";
+ tofDigiPar->Init();
+
+ fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
+ if (0 == fDigiPar) {
+ LOG(error) << "CbmTofSimpleClusterizer::InitParameters => Could not obtain "
+ "the CbmTofDigiPar ";
+ return false;
+ }
+ return true;
+}
+
+/************************************************************************************/
+bool CbmTaskTofHitFinder::InitCalibParameter()
+{
+ // dimension and initialize calib parameter
+ int32_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
+
+ fvCPSigPropSpeed.resize(iNbSmTypes);
+ for (int32_t iT = 0; iT < iNbSmTypes; iT++) {
+ int32_t iNbRpc = fDigiBdfPar->GetNbRpc(iT);
+ fvCPSigPropSpeed[iT].resize(iNbRpc);
+ for (int32_t iRpc = 0; iRpc < iNbRpc; iRpc++)
+ if (0.0 < fDigiBdfPar->GetSigVel(iT, 0, iRpc)) fvCPSigPropSpeed[iT][iRpc] = fDigiBdfPar->GetSigVel(iT, 0, iRpc);
+ else
+ fvCPSigPropSpeed[iT][iRpc] = fDigiBdfPar->GetSignalSpeed();
+ }
+
+ // Other calibration constants can be set here. Currently set to default values
+ fvCPTOff.resize(iNbSmTypes);
+ fvCPTotGain.resize(iNbSmTypes);
+ fvCPWalk.resize(iNbSmTypes);
+ for (int32_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
+ int32_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
+ int32_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
+ fvCPTOff[iSmType].resize(iNbSm * iNbRpc);
+ fvCPTotGain[iSmType].resize(iNbSm * iNbRpc);
+ fvCPWalk[iSmType].resize(iNbSm * iNbRpc);
+ for (int32_t iSm = 0; iSm < iNbSm; iSm++) {
+ for (int32_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
+ int32_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
+ fvCPTOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
+ fvCPTotGain[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
+ fvCPWalk[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
+ int32_t nbSide = 2 - fDigiBdfPar->GetChanType(iSmType, iRpc);
+ for (int32_t iCh = 0; iCh < iNbChan; iCh++) {
+ fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
+ fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
+ fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh].resize(nbSide);
+ for (int32_t iSide = 0; iSide < nbSide; iSide++) {
+ fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 0.;
+ fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 1.;
+ fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide].resize(numClWalkBinX);
+ for (int32_t iWx = 0; iWx < numClWalkBinX; iWx++) {
+ fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide][iWx] = 0.;
+ }
+ }
+ }
+ }
+ }
+ }
+ LOG(info) << "CbmTaskTofHitFinder::InitCalibParameter: defaults set";
+ LOG(info) << "CbmTaskTofHitFinder::InitCalibParameter: initialization done";
+ return true;
+}
+
+bool CbmTaskTofHitFinder::InitAlgos()
+{
+ /// Go to Top volume of the geometry in the GeoManager to make sure our nodes are found
+ gGeoManager->CdTop();
+
+ //Prepare storage vectors
+ int32_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
+ fStorDigiExp.resize(iNbSmTypes);
+ fStorDigiInd.resize(iNbSmTypes);
+
+ for (int32_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
+ int32_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
+ int32_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
+ fStorDigiExp[iSmType].resize(iNbSm * iNbRpc);
+ fStorDigiInd[iSmType].resize(iNbSm * iNbRpc);
+ }
+
+ // Create one algorithm per RPC and configure it with parameters
+ for (int32_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
+ int32_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
+ int32_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
+ for (int32_t iSm = 0; iSm < iNbSm; iSm++) {
+ for (int32_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
+ std::unique_ptr<HitFinderTofRpcPar> par(new HitFinderTofRpcPar());
+ par->FeeTimeRes = fDigiBdfPar->GetFeeTimeRes();
+ par->SysTimeRes = 0.080;
+ par->CPSigPropSpeed = fvCPSigPropSpeed[iSmType][iRpc];
+ par->outTimeFactor = 1.0;
+ par->numClWalkBinX = numClWalkBinX;
+ par->TOTMax = 5.E4;
+ par->TOTMin = 2.E4;
+ par->maxTimeDist = fDigiBdfPar->GetMaxTimeDist();
+ par->maxSpaceDist = fDigiBdfPar->GetMaxDistAlongCh();
+ par->numGaps = fDigiBdfPar->GetNbGaps(iSmType, iRpc);
+ par->gapSize = fDigiBdfPar->GetGapSize(iSmType, iRpc);
+ int32_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
+ par->fChanPar.resize(iNbChan);
+ for (int32_t iCh = 0; iCh < iNbChan; iCh++) {
+
+ //get channel info from iSmType, iSm, iRpc, iCh
+ CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, iSmType, iSm, iRpc, 0, iCh);
+ const int32_t iChId = fTofId->SetDetectorInfo(xDetInfo);
+ CbmTofCell* channelInfo = fDigiPar->GetCell(iChId);
+
+ //init Tof cell
+ TofCell& cell = par->fChanPar[iCh].cell;
+ cell.pos.SetX(channelInfo->GetX());
+ cell.pos.SetY(channelInfo->GetY());
+ cell.pos.SetZ(channelInfo->GetZ());
+ cell.sizeX = channelInfo->GetSizex();
+ cell.sizeY = channelInfo->GetSizey();
+
+ // prepare local->global trafo
+ gGeoManager->FindNode(channelInfo->GetX(), channelInfo->GetY(), channelInfo->GetZ());
+ double* rot_ptr = gGeoManager->GetCurrentMatrix()->GetRotationMatrix();
+ cell.rotation = ROOT::Math::Rotation3D(&rot_ptr[0], &rot_ptr[9]);
+
+ par->fChanPar[iCh].fvCPTOff = std::move(fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh]);
+ par->fChanPar[iCh].fvCPTotGain = std::move(fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh]);
+ par->fChanPar[iCh].fvCPWalk = std::move(fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh]);
+ par->fChanPar[iCh].address = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, 0, iSmType);
+ }
+ fAlgo[iSmType][iSm * iNbRpc + iRpc].SetParams(std::move(par));
+ }
+ }
+ }
+ return true;
+}
+
+pair<int32_t, int32_t> CbmTaskTofHitFinder::BuildClusters(CbmEvent* event)
+{
+ // --- MC Event info (input file, entry number, start time)
+ int32_t iInputNr = 0;
+ int32_t iEventNr = 0;
+ double dEventTime = 0.;
+ GetEventInfo(iInputNr, iEventNr, dEventTime);
+
+ // Local variables
+ int32_t nDigis = 0;
+ int32_t nHits = 0;
+
+ nDigis = (event ? event->GetNofData(ECbmDataType::kTofDigi) : fDigiMan->GetNofDigis(ECbmModuleId::kTof));
+ LOG(debug) << "Number of TOF digis: " << nDigis;
+
+ // Loop over the digis array and store the Digis in separate vectors for each RPC modules
+ for (int32_t iDigi = 0; iDigi < nDigis; iDigi++) {
+ const uint32_t digiIndex = (event ? event->GetIndex(ECbmDataType::kTofDigi, iDigi) : iDigi);
+ assert(fDigiMan->Get<CbmTofDigi>(digiIndex));
+ CbmTofDigi* pDigi = new CbmTofDigi(*(fDigiMan->Get<CbmTofDigi>(digiIndex)));
+ assert(pDigi);
+
+ // These are doubles in the digi class
+ const double smType = pDigi->GetType();
+ const double smId = pDigi->GetSm();
+ const double rpcId = pDigi->GetRpc();
+ const int32_t numRpc = fDigiBdfPar->GetNbRpc(smType);
+ fStorDigiExp[smType][smId * numRpc + rpcId].push_back(*pDigi);
+ fStorDigiInd[smType][smId * numRpc + rpcId].push_back(digiIndex);
+ }
+
+ const uint32_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
+
+ for (uint32_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
+ const uint32_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
+ const uint32_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
+ for (uint32_t iSm = 0; iSm < iNbSm; iSm++) {
+ for (uint32_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
+ //read digis
+ const uint32_t rpcIdx = iSm * iNbRpc + iRpc;
+ std::vector<CbmTofDigi>& digiExp = fStorDigiExp[iSmType][rpcIdx];
+ std::vector<int32_t>& digiInd = fStorDigiInd[iSmType][rpcIdx];
+
+ //call cluster finder
+ std::vector<TofCluster> clusters = fAlgo[iSmType][rpcIdx](digiExp, digiInd);
+
+ //Store hits and match
+ for (auto const& cluster : clusters) {
+ const int32_t hitIndex = fTofHitsColl->GetEntriesFast();
+ TVector3 hitpos = TVector3(cluster.globalPos.X(), cluster.globalPos.Y(), cluster.globalPos.Z());
+ TVector3 hiterr = TVector3(cluster.globalErr.X(), cluster.globalErr.Y(), cluster.globalErr.Z());
+ new ((*fTofHitsColl)[hitIndex])
+ CbmTofHit(cluster.detId, hitpos, hiterr, nHits, cluster.weightedTime, cluster.weightedTimeErr, 0, 0);
+ nHits++;
+ if (event) event->AddData(ECbmDataType::kTofHit, hitIndex);
+ CbmMatch* digiMatch = new CbmMatch();
+ for (uint32_t i = 0; i < cluster.vDigiIndRef.size(); i++) {
+ digiMatch->AddLink(CbmLink(0., cluster.vDigiIndRef.at(i), iEventNr, iInputNr));
+ }
+ new ((*fTofDigiMatchColl)[hitIndex]) CbmMatch(*digiMatch);
+ delete digiMatch;
+ }
+ }
+ }
+ }
+
+ return std::make_pair(nDigis, nHits);
+}
+
+void CbmTaskTofHitFinder::GetEventInfo(int32_t& inputNr, int32_t& eventNr, double& eventTime)
+{
+ // --- In a FairRunAna, take the information from FairEventHeader
+ if (FairRunAna::Instance()) {
+ FairEventHeader* event = FairRunAna::Instance()->GetEventHeader();
+ inputNr = event->GetInputFileId();
+ eventNr = event->GetMCEntryNumber();
+ eventTime = event->GetEventTime();
+ }
+
+ // --- In a FairRunSim, the input number and event time are always zero;
+ // --- only the event number is retrieved.
+ else {
+ if (!FairRunSim::Instance()) LOG(fatal) << GetName() << ": neither SIM nor ANA run.";
+ FairMCEventHeader* event = FairRunSim::Instance()->GetMCEventHeader();
+ inputNr = 0;
+ eventNr = event->GetEventID();
+ eventTime = 0.;
+ }
+}
+
+ClassImp(CbmTaskTofHitFinder)
diff --git a/reco/tasks/CbmTaskTofHitFinder.h b/reco/tasks/CbmTaskTofHitFinder.h
new file mode 100644
index 0000000000000000000000000000000000000000..e192587d494f351552416af7527d772ec7b55f94
--- /dev/null
+++ b/reco/tasks/CbmTaskTofHitFinder.h
@@ -0,0 +1,163 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Pierre-Alain Loizeau, Norbert Herrmann */
+
+/** @class CbmTaskTofHitFinder
+ ** @brief Simple Cluster building and hit producing for CBM ToF using Digis as input
+ ** @author Pierre-Alain Loizeau <loizeau@physi.uni-heidelberg.de>
+ ** @version 1.0
+ **/
+#ifndef CBMTOFHITFINDER_H
+#define CBMTOFHITFINDER_H 1
+
+// TOF Classes and includes
+// Input/Output
+class CbmTofDigi;
+class CbmMatch;
+// Geometry
+class CbmTofGeoHandler;
+class CbmTofDetectorId;
+class CbmTofDigiPar;
+class CbmTofDigiBdfPar;
+class CbmTofCell;
+class CbmDigiManager;
+class CbmEvent;
+
+#include "HitFinderTof.h"
+
+// FAIR classes and includes
+#include "FairTask.h"
+
+// ROOT Classes and includes
+class TClonesArray;
+#include "TStopwatch.h"
+#include "TTimeStamp.h"
+
+// C++ Classes and includes
+#include <vector>
+
+class CbmTaskTofHitFinder : public FairTask {
+public:
+ /**
+ ** @brief Constructor.
+ **/
+ CbmTaskTofHitFinder();
+
+ /**
+ ** @brief Constructor.
+ **/
+ CbmTaskTofHitFinder(const char* name, int32_t verbose = 1);
+ /**
+ ** @brief Destructor.
+ **/
+ virtual ~CbmTaskTofHitFinder();
+
+ /**
+ ** @brief Inherited from FairTask.
+ **/
+ virtual InitStatus Init();
+
+ /**
+ ** @brief Inherited from FairTask.
+ **/
+ virtual void SetParContainers();
+
+ /**
+ ** @brief Inherited from FairTask.
+ **/
+ virtual void Exec(Option_t* option);
+
+ /**
+ ** @brief Inherited from FairTask.
+ **/
+ virtual void Finish();
+
+protected:
+private:
+ /**
+ ** @brief Copy constructor.
+ **/
+ CbmTaskTofHitFinder(const CbmTaskTofHitFinder&);
+ /**
+ ** @brief Copy operator.
+ **/
+ CbmTaskTofHitFinder& operator=(const CbmTaskTofHitFinder&);
+
+ // Functions common for all clusters approximations
+ /**
+ ** @brief Recover pointer on input TClonesArray: TofPoints, TofDigis...
+ **/
+ bool RegisterInputs();
+ /**
+ ** @brief Create and register output TClonesArray of Tof Hits.
+ **/
+ bool RegisterOutputs();
+ /**
+ ** @brief Initialize other parameters not included in parameter classes.
+ **/
+ bool InitParameters();
+ /**
+ ** @brief Initialize other parameters not included in parameter classes.
+ **/
+ bool InitCalibParameter();
+
+ /**
+ ** @brief Create one algo object for each RPC
+ **/
+ bool InitAlgos();
+
+ /**
+ ** @brief Build clusters out of ToF Digis and store the resulting info in a TofHit.
+ **/
+ std::pair<int32_t, int32_t> BuildClusters(CbmEvent* event);
+
+ /**
+ ** @brief Retrieve event info from run manager to properly fill the CbmLink objects.
+ **/
+ void GetEventInfo(int32_t& inputNr, int32_t& eventNr, double& eventTime);
+
+ // ToF geometry variables
+ CbmTofGeoHandler* fGeoHandler;
+ CbmTofDetectorId* fTofId = nullptr;
+ CbmTofDigiPar* fDigiPar = nullptr;
+ CbmTofDigiBdfPar* fDigiBdfPar = nullptr;
+
+ // Input variables
+ CbmDigiManager* fDigiMan = nullptr;
+ TClonesArray* fEvents = nullptr;
+
+ // Output variables
+ TClonesArray* fTofHitsColl = nullptr; // TOF hits
+ TClonesArray* fTofDigiMatchColl = nullptr; // TOF Digis
+
+ // Hit finder algo
+ std::map<uint32_t, std::map<uint32_t, cbm::algo::HitFinderTof>> fAlgo = {}; //[nbType][nbSm*nbRpc]
+
+ // Intermediate storage variables
+ std::vector<std::vector<std::vector<CbmTofDigi>>> fStorDigiExp; //[nbType][nbSm*nbRpc][nbCh][nDigis]
+ std::vector<std::vector<std::vector<int32_t>>> fStorDigiInd; //[nbType][nbSm*nbRpc][nbCh][nDigis]
+
+ //Calibration parameters
+ std::vector<std::vector<double>> fvCPSigPropSpeed; //[nSMT][nRpc]
+ std::vector<std::vector<std::vector<std::vector<double>>>> fvCPTOff; //[nSMT][nRpc][nCh][nbSide]
+ std::vector<std::vector<std::vector<std::vector<double>>>> fvCPTotGain; //[nSMT][nRpc][nCh][nbSide]
+ std::vector<std::vector<std::vector<std::vector<std::vector<double>>>>>
+ fvCPWalk; //[nSMT][nRpc][nCh][nbSide][nbWalkBins]
+
+ // Control
+ TTimeStamp fStart;
+ TTimeStamp fStop;
+
+ // --- Run counters
+ TStopwatch fTimer; ///< ROOT timer
+ int32_t fiNofTs = 0; ///< Number of processed timeslices
+ int32_t fiNofEvents = 0; ///< Total number of events processed
+ double fNofDigisAll = 0.; ///< Total number of TOF digis in input
+ double fNofDigisUsed = 0.; ///< Total number of Tof Digis processed
+ double fdNofHitsTot = 0.; ///< Total number of hits produced
+ double fdTimeTot = 0.; ///< Total execution time
+
+ ClassDef(CbmTaskTofHitFinder, 1);
+};
+
+#endif // CBMTOFHITFINDER_H