diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt
index 6f425a623a31eb885bd18164b07c7ce17d798d0b..2841c8f7e2f44a43f9e988a7d8d537694afe5e15 100644
--- a/algo/CMakeLists.txt
+++ b/algo/CMakeLists.txt
@@ -35,6 +35,7 @@ set(SRCS
detectors/much/ReadoutConfig.cxx
detectors/much/Unpack.cxx
detectors/tof/HitFinder.cxx
+ detectors/tof/ClusterizerTof.cxx
detectors/tof/ReadoutConfig.cxx
detectors/tof/Unpack.cxx
detectors/bmon/ReadoutConfig.cxx
diff --git a/algo/detectors/tof/ClusterizerTof.cxx b/algo/detectors/tof/ClusterizerTof.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7f12af2f7dc79036e9aec95ad723ed2703de136f
--- /dev/null
+++ b/algo/detectors/tof/ClusterizerTof.cxx
@@ -0,0 +1,445 @@
+/* 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 "ClusterizerTof.h"
+
+// TOF Classes and includes
+#include "CbmTofDigi.h"
+
+// C++ Classes and includes
+#include <algorithm>
+#include <iomanip>
+#include <iostream>
+#include <map>
+
+namespace cbm::algo
+{
+
+ ClusterizerTof::resultType ClusterizerTof::operator()(inputType digisIn)
+ {
+ //std::vector<inputType> input = calibrateDigis(digisIn);
+ std::vector<inputType> input = chanSortDigis(digisIn);
+ return buildClusters(input);
+ }
+
+ std::vector<ClusterizerTof::inputType> ClusterizerTof::chanSortDigis(inputType& digisIn)
+ {
+ std::vector<inputType> result(fParams.fChanPar.size()); //[nbCh][nDigis]
+
+ for (size_t iDigi = 0; iDigi < digisIn.size(); iDigi++) {
+ CbmTofDigi* pDigi = digisIn[iDigi].first;
+ const int32_t index = digisIn[iDigi].second;
+ const double chan = pDigi->GetChannel();
+ result[chan].push_back(std::make_pair(new CbmTofDigi(*pDigi), index));
+ }
+ return result;
+ }
+
+ //This calibration cannot reproduce the non RPC-parallel version, as the digi indices
+ //will be invalidated whenever digis are discarded in the dead-time check.
+ std::vector<ClusterizerTof::inputType> ClusterizerTof::calibrateDigis(inputType& digisIn)
+ {
+ std::vector<inputType> result(fParams.fChanPar.size());
+ const int32_t numClWalkBinX = fParams.numClWalkBinX;
+ const double TOTMax = fParams.TOTMax;
+ const double TOTMin = fParams.TOTMin;
+
+ // channel deadtime map
+ std::map<int32_t, double> mChannelDeadTime;
+
+ for (size_t iDigi = 0; iDigi < digisIn.size(); iDigi++) {
+ CbmTofDigi* pDigi = digisIn[iDigi].first;
+ const double chan = pDigi->GetChannel();
+ const double side = pDigi->GetSide();
+ const double type = pDigi->GetType();
+
+ if (fParams.swapChannelSides && 5 != type && 8 != type) {
+ pDigi->SetAddress(pDigi->GetSm(), pDigi->GetRpc(), chan, (0 == side) ? 1 : 0, type);
+ }
+
+ // Check dead time
+ const int32_t addr = pDigi->GetAddress();
+ auto it = mChannelDeadTime.find(addr);
+ if (it != mChannelDeadTime.end() && pDigi->GetTime() <= it->second) {
+ it->second = pDigi->GetTime() + fParams.channelDeadtime;
+ continue;
+ }
+ mChannelDeadTime[addr] = pDigi->GetTime() + fParams.channelDeadtime;
+
+ // Create calibrated digi
+ CbmTofDigi* pCalDigi = new CbmTofDigi(*pDigi);
+
+ // calibrate Digi Time
+ pCalDigi->SetTime(pCalDigi->GetTime() - fParams.fChanPar[chan].fvCPTOff[side]);
+
+ // subtract Offset
+ double charge = pCalDigi->GetTot() - fParams.fChanPar[chan].fvCPTotOff[side];
+ if (charge < 0.001) charge = 0.001;
+
+ // calibrate Digi ToT
+ pCalDigi->SetTot(charge * fParams.fChanPar[chan].fvCPTotGain[side]);
+
+ // walk correction
+ const double chargeBinSize = (TOTMax - TOTMin) / numClWalkBinX;
+ int32_t iWx = (int32_t)((pCalDigi->GetTot() - TOTMin) / chargeBinSize);
+
+ if (0 > iWx) { iWx = 0; }
+ if (iWx >= numClWalkBinX) { iWx = numClWalkBinX - 1; }
+
+ std::vector<double>& walk = fParams.fChanPar[chan].fvCPWalk[side];
+ const double dDTot = (pCalDigi->GetTot() - TOTMin) / chargeBinSize - (double) iWx - 0.5;
+ double dWT = walk[iWx];
+
+ // linear interpolation to next bin //memory leak??? (D.Smith 10.8.23: Clarify this comment!)
+ if (dDTot > 0) {
+ if (iWx < numClWalkBinX - 1) { dWT += dDTot * (walk[iWx + 1] - walk[iWx]); }
+ }
+ else {
+ if (0 < iWx) { dWT -= dDTot * (walk[iWx - 1] - walk[iWx]); }
+ }
+ pCalDigi->SetTime(pCalDigi->GetTime() - dWT); // calibrate Digi Time
+
+ const int32_t index = digisIn[iDigi].second;
+ result[chan].push_back(std::make_pair(pCalDigi, index));
+ }
+
+ // D.Smith 10.8.23: Sorting below might not be needed. Check with P.A.
+ /// Sort the buffers of hits due to the time offsets applied
+ for (size_t chan = 0; chan < fParams.fChanPar.size(); chan++) {
+ std::sort(result[chan].begin(), result[chan].end(),
+ [](const std::pair<CbmTofDigi*, int32_t>& a, const std::pair<CbmTofDigi*, int32_t>& b) -> bool {
+ return a.first->GetTime() < b.first->GetTime();
+ });
+ }
+
+ return result;
+ }
+
+
+ ClusterizerTof::resultType ClusterizerTof::buildClusters(std::vector<inputType>& input)
+ {
+ // Hit variables
+ TofCluster cluster;
+ std::vector<TofCluster> clustersOut;
+
+ // Reference cell of a cluster
+ TofCell* cell = nullptr;
+
+ // Last Channel Temp variables
+ int32_t lastChan = -1;
+ double lastPosY = 0.0;
+ double lastTime = 0.0;
+
+ const size_t numChan = fParams.fChanPar.size();
+
+ for (int32_t chan = 0; chan < numChan; chan++) {
+ if (fParams.fDeadStrips & (1 << chan)) { continue; } // skip over dead channels
+
+ inputType& storDigi = input[chan];
+
+ while (1 < storDigi.size()) {
+ while (storDigi[0].first->GetSide() == storDigi[1].first->GetSide()) {
+ // Not one Digi of each end!
+ uint8_t offset = 0;
+
+ /* D.Smith 14.8.23: This condition is never needed due to time-sorted input (check with PAL and VF)
+ // use digi that is closer to the last one
+ if (storDigi.size() > 2 && storDigi[2].first->GetSide() != storDigi[0].first->GetSide()
+ && storDigi[2].first->GetTime() - storDigi[0].first->GetTime()
+ <= storDigi[2].first->GetTime() - storDigi[1].first->GetTime()) {
+ offset++;
+ }
+*/
+ storDigi.erase(storDigi.begin() + offset);
+ if (2 > storDigi.size()) { break; }
+ } // same condition side end
+ if (2 > storDigi.size()) { break; }
+
+ // 2 Digis = both sides present
+ cell = &fParams.fChanPar[chan].cell;
+ CbmTofDigi* xDigiA = storDigi[0].first;
+ CbmTofDigi* xDigiB = storDigi[1].first;
+
+ // use local coordinates, (0,0,0) is in the center of counter ?
+ ROOT::Math::XYZVector pos(((-(double) numChan / 2. + (double) chan) + 0.5) * cell->sizeX, 0., 0.);
+
+ double timeDif = xDigiA->GetTime() - xDigiB->GetTime();
+
+ pos.SetY(fParams.fSigVel * timeDif * 0.5); // A is the top side, B is the bottom side
+ if (xDigiA->GetSide() != 1.) { pos.SetY(-pos.Y()); } // B is the bottom side, A is the top side
+
+ while (storDigi.size() > 2 && std::abs(pos.Y()) > cell->sizeY * fParams.fPosYMaxScal) {
+
+ CbmTofDigi* xDigiC = storDigi[2].first;
+
+ double timeDifN = 0;
+ if (xDigiC->GetSide() == xDigiA->GetSide()) { timeDifN = xDigiC->GetTime() - xDigiB->GetTime(); }
+ else {
+ timeDifN = xDigiC->GetTime() - xDigiA->GetTime();
+ }
+
+ double posYN = fParams.fSigVel * timeDifN * 0.5;
+ if (xDigiC->GetSide() != 1.) { posYN *= -1.; }
+
+ if (std::abs(posYN) >= std::abs(pos.Y())) { break; }
+ pos.SetY(posYN);
+
+ if (xDigiC->GetSide() == xDigiA->GetSide()) {
+ xDigiA = xDigiC;
+ storDigi.erase(storDigi.begin());
+ }
+ else {
+ xDigiB = xDigiC;
+
+ //D.Smith 14.8.23: Shouldn't we erase the "+1" element here instead?
+ //Doesn't seem to have an effect on the result regardless.
+ //"+1" would be better for an iterator based implementation.
+
+ //D.Smith 25.8.23: Actually, for some strange reason I don't fully understand, this
+ //erase operation doesn't matter at all, even though this branch is sometimes called.
+ //Perhaps this is data dependent? Erasing the zero element also works. This would
+ //be best for iterators.
+ storDigi.erase(storDigi.begin() + 2);
+ }
+
+ } //while loop end
+
+ if (std::abs(pos.Y()) > cell->sizeY * fParams.fPosYMaxScal) { // remove both digis
+ storDigi.erase(storDigi.begin());
+ continue;
+ }
+ // The "Strip" time is the mean time between each end
+ double time = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
+
+ // Weight is the total charge => sum of both ends ToT
+ double totSum = xDigiA->GetTot() + xDigiB->GetTot();
+
+ // 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.fdMaxTimeDist && lastChan == chan - 1
+ && std::abs(pos.Y() - lastPosY) < fParams.fdMaxSpaceDist)) {
+
+ cluster.normalize(fParams.fTimeRes);
+ cluster.finalize(*cell, fParams);
+ clustersOut.push_back(cluster);
+ cluster.reset();
+ }
+ }
+ else {
+ // first fired strip in this RPC
+ cluster.reset();
+ }
+ cluster.add(pos, time, totSum, totSum, storDigi[0].second, storDigi[1].second);
+
+ storDigi.erase(storDigi.begin());
+ storDigi.erase(storDigi.begin());
+
+ lastChan = chan;
+ lastPosY = pos.Y();
+ lastTime = time;
+ if (AddNextChan(input, lastChan, cluster, clustersOut)) { cluster.reset(); }
+ } // while( 1 < storDigi.size() )
+ storDigi.clear(); //D.Smith 11.8.23: In rare cases, a single digi remains and is deleted here.
+ } // for( int32_t chan = 0; chan < iNbCh; chan++ )
+
+ // Now check if another hit/cluster is started
+ // and save it if it's the case.
+ if (0 < cluster.numChan()) {
+ cluster.normalize(fParams.fTimeRes);
+ cluster.finalize(*cell, fParams);
+ clustersOut.push_back(cluster);
+ }
+ //std::cout << "hits " << fiNbHits << std::endl;
+ return clustersOut;
+ }
+
+
+ bool ClusterizerTof::AddNextChan(std::vector<inputType>& input, int32_t lastChan, TofCluster& cluster,
+ std::vector<TofCluster>& clustersOut)
+ {
+ //D.Smith 25.8.23: Why are "C" digis (position "2") not considered here?
+
+ //const int32_t iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
+ size_t numChan = fParams.fChanPar.size();
+ int32_t chan = lastChan + 1;
+ inputType& storDigi = input[chan];
+
+ while (fParams.fDeadStrips & (1 << chan)) {
+ chan++;
+ if (chan >= numChan) { return false; }
+ }
+
+ if (chan == numChan) { return false; }
+ if (0 == storDigi.size()) { return false; }
+
+ bool addedHit = false;
+ for (size_t i1 = 0; i1 < storDigi.size() - 1; i1++) {
+ if (addedHit) { break; }
+
+ size_t i2 = i1;
+ while (!addedHit && ++i2 < storDigi.size()) {
+
+ const CbmTofDigi* xDigiA = storDigi[i1].first;
+ const CbmTofDigi* xDigiB = storDigi[i2].first;
+
+ if (xDigiA->GetSide() == xDigiB->GetSide()) { continue; }
+ const double time = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
+
+ if (std::abs(time - cluster.weightedTime / cluster.weightsSum) >= fParams.fdMaxTimeDist) { continue; }
+
+ TofCell* cell = &fParams.fChanPar[chan].cell;
+
+ const double timeDif = xDigiA->GetTime() - xDigiB->GetTime();
+ double posY = fParams.fSigVel * timeDif * 0.5;
+
+ if (1 != xDigiA->GetSide()) { posY *= -1.; }
+
+ if (std::abs(posY - cluster.weightedPos.Y() / cluster.weightsSum) >= fParams.fdMaxSpaceDist) { continue; }
+
+ // append digi pair to current cluster
+ const double posX = ((-(double) numChan / 2. + chan) + 0.5) * cell->sizeX;
+ const double totSum = xDigiA->GetTot() + xDigiB->GetTot();
+
+ ROOT::Math::XYZVector pos(posX, posY, 0.);
+ cluster.add(pos, time, totSum, totSum, storDigi[i1].second, storDigi[i2].second);
+
+ // remove selected digis from pool
+ storDigi.erase(storDigi.begin() + i2);
+ storDigi.erase(storDigi.begin() + i1);
+
+ if (AddNextChan(input, chan, cluster, clustersOut)) {
+ return true; // signal hit was already added
+ }
+ addedHit = true;
+ }
+ }
+
+ TofCell cell = fParams.fChanPar[0].cell; //D.Smith 17.8.23: This is equivalent to using iDetId, see below
+ //D.Smith 10.8.23: Why pass iDetId here and not iChId?
+ cluster.normalize(fParams.fTimeRes);
+ cluster.finalize(cell, fParams);
+ clustersOut.push_back(cluster);
+ return true;
+ }
+
+
+ /*
+ ClusterizerTof::resultType ClusterizerTof::buildClusters(inputType& input)
+ {
+ // Intermediate storage variables
+ std::vector<std::vector<CbmTofDigi*>>& digisExp = input.first; //[nbCh][nDigis]
+ std::vector<std::vector<int32_t>>& digisInd = input.second; //[nbCh][nDigis]
+
+ // 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;
+
+ for (int32_t chan = 0; chan < fParams.fChanPar.size(); chan++) {
+
+ std::vector<CbmTofDigi*>& storDigiExp = digisExp[chan];
+ std::vector<int32_t>& storDigiInd = digisInd[chan];
+ ClusterizerTofChanPar& chanPar = fParams.fChanPar[chan];
+
+ auto digiExpIt = storDigiExp.begin();
+ auto digiIndIt = storDigiInd.begin();
+
+ while (1 < std::distance(digiExpIt, storDigiExp.end())) {
+ while ((*digiExpIt)->GetSide() == (*std::next(digiExpIt, 1))->GetSide()) {
+ // Not one Digi of each end!
+ numSameSide++;
+ digiExpIt++;
+ digiIndIt++;
+ if (2 > std::distance(digiExpIt, storDigiExp.end())) break;
+ }
+ if (2 > std::distance(digiExpIt, storDigiExp.end())) break; // 2 Digis = both sides present
+
+ TofCell* channelInfo = &chanPar.cell;
+ CbmTofDigi* xDigiA = *digiExpIt;
+ CbmTofDigi* xDigiB = *std::next(digiExpIt, 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)
+ digiExpIt++;
+ digiIndIt++;
+ 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, *digiIndIt, *std::next(digiIndIt, 1));
+ digiExpIt += 2;
+ digiIndIt += 2;
+
+ lastChan = chan;
+ lastPosY = pos.Y();
+ lastTime = time;
+ } // while (1 < storDigiExp.size())
+ } // 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/ClusterizerTof.h b/algo/detectors/tof/ClusterizerTof.h
new file mode 100644
index 0000000000000000000000000000000000000000..33f254da5c019cf974f2dc207c271723c1bcd78d
--- /dev/null
+++ b/algo/detectors/tof/ClusterizerTof.h
@@ -0,0 +1,194 @@
+/* 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 */
+
+/*
+ This algo was based on CbmTofSimpClusterizer, which can be used only for simulation of the main setup and
+ is as the name implies a simplified solution.
+ A later step will be the replacement with a version based on CbmTofEventClusterizer, which is the version
+ currently maintained, based on what we learned from real data at mCBM.
+ This step will be required to apply the algo to real/online data and to prepare
+ our simulations for first CBM beam
+*/
+
+#ifndef CLUSTERIZERTOF_H
+#define CLUSTERIZERTOF_H
+
+// TOF Classes and includes
+class CbmTofDigi;
+
+// ROOT Classes and includes
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+// C++ Classes and includes
+#include <memory>
+#include <vector>
+
+#include <cmath>
+
+namespace cbm::algo
+{
+ struct TofCell {
+ double sizeX, sizeY;
+ ROOT::Math::XYZVector pos;
+ ROOT::Math::Rotation3D rotation;
+ };
+
+ struct ClusterizerTofChanPar {
+ std::vector<double> fvCPTOff; //[nbSide]
+ std::vector<double> fvCPTotGain; //[nbSide]
+ std::vector<double> fvCPTotOff; //[nbSide]
+ std::vector<std::vector<double>> fvCPWalk; //[nbSide][nbWalkBins]
+ int32_t address; //unique address
+ TofCell cell;
+ };
+
+ struct ClusterizerTofRpcPar {
+ uint32_t fDeadStrips;
+ double fPosYMaxScal;
+ double fdMaxTimeDist;
+ double fdMaxSpaceDist;
+ double fSigVel;
+ double fTimeRes;
+ int32_t numClWalkBinX;
+ double TOTMax;
+ double TOTMin;
+ bool swapChannelSides;
+ double channelDeadtime;
+ double fCPTOffYBinWidth;
+ double fCPTOffYRange;
+ std::vector<double> fCPTOffY; //[nBin]
+ std::vector<ClusterizerTofChanPar> 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.5, 0.5, 0.5); //D.Smith 15.8.23: Set this back to zero eventually.
+ 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 ClusterizerTofRpcPar& 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(0.5, 0.5, 0.5);
+
+ //store results
+ //globalErr = rotMatrix(hitErrLocal);
+ globalErr = hitErrLocal;
+ channel = par.fChanPar.size() / 2 + floor(weightedPos.X() / trafoCell.sizeX);
+
+ // D.Smith 15.8.23: This channel correction doesn't seem to be needed
+ // if (channel < 0) channel = 0;
+ // if (channel > par.fChanPar.size() - 1) channel = par.fChanPar.size() - 1;
+ detId = par.fChanPar[channel].address;
+
+ //Calibrate hit time if applicable
+ if (par.fCPTOffYRange == 0) { return; }
+ const double dBin = (weightedPos.Y() + par.fCPTOffYRange) / par.fCPTOffYBinWidth;
+ const int i1 = floor(dBin);
+ const int i2 = ceil(dBin);
+ weightedTime -= par.fCPTOffY[i1] + (dBin - i1) * (par.fCPTOffY[i2] - par.fCPTOffY[i1]);
+ }
+ };
+
+ class ClusterizerTof {
+ public:
+ typedef std::vector<TofCluster> resultType;
+ typedef std::vector<std::pair<CbmTofDigi*, int32_t>> inputType;
+
+ /**
+ ** @brief Constructor.
+ **/
+ ClusterizerTof() {};
+
+ /**
+ ** @brief Destructor.
+ **/
+ ~ClusterizerTof() {};
+
+ /**
+ ** @brief Build clusters out of ToF Digis and store the resulting info in a TofHit.
+ **/
+ resultType operator()(inputType digisIn);
+
+ /** @brief Set the parameter container
+ ** @param params Vectorer to parameter container
+ **/
+ void SetParams(std::unique_ptr<ClusterizerTofRpcPar> params) { fParams = *(std::move(params)); }
+
+ private:
+ ClusterizerTofRpcPar fParams = {}; ///< Parameter container
+
+ std::vector<inputType> calibrateDigis(inputType& digisIn);
+ std::vector<inputType> chanSortDigis(inputType& digisIn);
+
+ resultType buildClusters(std::vector<inputType>& input);
+
+ bool AddNextChan(std::vector<inputType>& input, int32_t iLastChan, TofCluster& cluster,
+ std::vector<TofCluster>& clustersOut);
+
+ int32_t numSameSide; // Digis quality
+ };
+
+} /* namespace cbm::algo */
+
+#endif // CLUSTERIZERTOF_H
diff --git a/reco/tasks/CMakeLists.txt b/reco/tasks/CMakeLists.txt
index 4743544232b5b10e28bb4d63ff4421a13bdbc5a6..300bf6ed0f3ef5dbf2a93cc3d079302550ce083f 100644
--- a/reco/tasks/CMakeLists.txt
+++ b/reco/tasks/CMakeLists.txt
@@ -15,6 +15,7 @@ set(SRCS
CbmTaskMakeRecoEvents.cxx
CbmTaskTriggerDigi.cxx
CbmTaskTofHitFinder.cxx
+ CbmTaskTofClusterizer.cxx
CbmTaskUnpack.cxx
CbmTaskUnpackXpu.cxx
)
diff --git a/reco/tasks/CbmRecoTasksLinkDef.h b/reco/tasks/CbmRecoTasksLinkDef.h
index 16e011359c06c17317bcf9b373304d58bec362f6..8548dc9c4608bf4df73d7b4026568537187bffc5 100644
--- a/reco/tasks/CbmRecoTasksLinkDef.h
+++ b/reco/tasks/CbmRecoTasksLinkDef.h
@@ -18,6 +18,7 @@
#pragma link C++ class CbmTaskEventsCloneInToOut + ;
#pragma link C++ class CbmTaskMakeRecoEvents + ;
#pragma link C++ class CbmTaskTofHitFinder + ;
+#pragma link C++ class CbmTaskTofClusterizer + ;
#pragma link C++ class CbmTaskTriggerDigi + ;
#pragma link C++ class CbmTaskUnpack + ;
#pragma link C++ class CbmTaskUnpackXpu + ;
diff --git a/reco/tasks/CbmTaskTofClusterizer.cxx b/reco/tasks/CbmTaskTofClusterizer.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7e90676aeb616eea7f229eb15efc6219098d6238
--- /dev/null
+++ b/reco/tasks/CbmTaskTofClusterizer.cxx
@@ -0,0 +1,1082 @@
+/* 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 "CbmTaskTofClusterizer.h"
+
+// TOF Classes and includes
+#include "CbmBmonDigi.h" // in cbmdata/bmon
+#include "CbmDigiManager.h"
+#include "CbmEvent.h"
+#include "CbmMatch.h"
+#include "CbmTofAddress.h" // in cbmdata/tof
+#include "CbmTofCell.h" // in tof/TofData
+#include "CbmTofCreateDigiPar.h" // in tof/TofTools
+#include "CbmTofDetectorId_v12b.h" // in cbmdata/tof
+#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
+#include "CbmTsEventHeader.h"
+
+// FAIR classes and includes
+#include "FairRootFileSink.h"
+#include "FairRootManager.h"
+#include "FairRunAna.h"
+#include "FairRuntimeDb.h"
+#include <Logger.h>
+
+// ROOT Classes and includes
+#include "TClonesArray.h"
+#include "TGeoManager.h"
+#include "TGeoPhysicalNode.h"
+#include "TH2.h"
+#include "TProfile.h"
+#include "TROOT.h"
+#include "TStopwatch.h"
+#include "TVector3.h"
+
+// C++ Classes and includes
+#include <iomanip>
+#include <vector>
+
+// Globals
+
+
+/************************************************************************************/
+CbmTaskTofClusterizer::CbmTaskTofClusterizer() : CbmTaskTofClusterizer("TestbeamClusterizer", 0, 0) {}
+
+CbmTaskTofClusterizer::CbmTaskTofClusterizer(const char* name, int32_t verbose, bool writeDataInOut)
+ : FairTask(TString(name), verbose)
+ , fTofId(NULL)
+ , fDigiPar(NULL)
+ , fDigiBdfPar(NULL)
+ , fTsHeader(NULL)
+ , fDigiMan(nullptr)
+ , fEventsColl(nullptr)
+ , fbWriteHitsInOut(writeDataInOut)
+ , fbWriteDigisInOut(writeDataInOut)
+ , fTofHitsColl(NULL)
+ , fTofDigiMatchColl(NULL)
+ , fTofHitsCollOut(NULL)
+ , fTofDigiMatchCollOut(NULL)
+ , fiNbHits(0)
+ , fStorDigi()
+ , fvCPTOff()
+ , fvCPTotGain()
+ , fvCPTotOff()
+ , fvCPWalk()
+ , fvCPTOffY()
+ , fvCPTOffYBinWidth()
+ , fvCPTOffYRange()
+ , fvDeadStrips()
+ , fCalMode(0)
+ , fDutId(-1)
+ , fPosYMaxScal(1.)
+ , fTotMax(100.)
+ , fTotMin(0.)
+ , fTotOff(0.)
+ , fTotMean(0.)
+ , fMaxTimeDist(1.)
+ , fdChannelDeadtime(0.)
+ , fCalParFileName("")
+ , fdTOTMax(50.)
+ , fdTOTMin(0.)
+ , fdTTotMean(2.)
+ , fdMaxTimeDist(0.)
+ , fdMaxSpaceDist(0.)
+ , fdModifySigvel(1.0)
+ , fdEvent(0)
+ , fbSwapChannelSides(false)
+ , fiOutputTreeEntry(0)
+ , fiFileIndex(0)
+{
+}
+
+CbmTaskTofClusterizer::~CbmTaskTofClusterizer() {}
+
+/************************************************************************************/
+// FairTasks inherited functions
+InitStatus CbmTaskTofClusterizer::Init()
+{
+ LOG(info) << "CbmTaskTofClusterizer initializing... expect Digis in ns units! ";
+ 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 CbmTaskTofClusterizer::SetParContainers()
+{
+ LOG(info) << "=> Get the digi parameters for tof";
+ FairRunAna* ana = FairRunAna::Instance();
+ FairRuntimeDb* rtdb = ana->GetRuntimeDb();
+
+ fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
+
+ LOG(info) << "found " << fDigiPar->GetNrOfModules() << " cells ";
+ fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
+}
+
+void CbmTaskTofClusterizer::Exec(Option_t* option)
+{
+ if (fTofCalDigiVecOut) fTofCalDigiVecOut->clear();
+ if (fEventsColl) {
+ if (NULL != fTsHeader)
+ LOG(info) << "New Ts " << iNbTs << ", size " << fEventsColl->GetSize() << " at " << fTsHeader->GetTsStartTime()
+ << " with " << fEventsColl->GetEntriesFast() << " events, " << fDigiMan->GetNofDigis(ECbmModuleId::kTof)
+ << " TOF digis + " << fDigiMan->GetNofDigis(ECbmModuleId::kT0) << " T0 digis ";
+
+ TStopwatch timerTs;
+ timerTs.Start();
+
+ iNbTs++;
+ fiHitStart = 0;
+ int32_t iNbHits = 0;
+ int32_t iNbCalDigis = 0;
+ fTofDigiMatchCollOut->Delete(); // costly, FIXME
+ fTofHitsCollOut->Delete(); // costly, FIXME
+
+ for (int32_t iEvent = 0; iEvent < fEventsColl->GetEntriesFast(); iEvent++) {
+ CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEventsColl->At(iEvent));
+ fTofDigiVec.clear();
+ LOG(debug) << "TS event " << iEvent << " with " << tEvent->GetNofData(ECbmDataType::kT0Digi) << " T0 and "
+ << tEvent->GetNofData(ECbmDataType::kTofDigi) << " Tof digis ";
+
+ for (int32_t iDigi = 0; iDigi < tEvent->GetNofData(ECbmDataType::kT0Digi); iDigi++) {
+ int32_t iDigiIndex = static_cast<int32_t>(tEvent->GetIndex(ECbmDataType::kT0Digi, iDigi));
+ CbmTofDigi tDigi(fDigiMan->Get<CbmBmonDigi>(iDigiIndex));
+ if (tDigi.GetType() != 5) {
+ tDigi.SetAddress(0, 0, 0, 0, 5); // convert to Tof Address
+ }
+ if (tDigi.GetSide() == 1) { // HACK for May22 setup
+ tDigi.SetAddress(tDigi.GetSm(), tDigi.GetRpc(), tDigi.GetChannel() + 6, 0, tDigi.GetType());
+ }
+ fTofDigiVec.push_back(tDigi);
+ }
+ for (int32_t iDigi = 0; iDigi < tEvent->GetNofData(ECbmDataType::kTofDigi); iDigi++) {
+ int32_t iDigiIndex = static_cast<int32_t>(tEvent->GetIndex(ECbmDataType::kTofDigi, iDigi));
+ const CbmTofDigi* tDigi = fDigiMan->Get<CbmTofDigi>(iDigiIndex);
+ fTofDigiVec.push_back(CbmTofDigi(*tDigi));
+ }
+
+ ExecEvent(option, tEvent);
+
+ // --- In event-by-event mode: copy caldigis, hits and matches to output array and register them to event
+ uint uDigi0 = fTofCalDigiVecOut->size(); //starting index of current event
+ LOG(debug) << "Append " << fTofCalDigiVec->size() << " CalDigis to event " << tEvent->GetNumber();
+ for (uint32_t index = 0; index < fTofCalDigiVec->size(); index++) {
+ CbmTofDigi* tDigi = &(fTofCalDigiVec->at(index));
+ tEvent->AddData(ECbmDataType::kTofCalDigi, iNbCalDigis);
+ fTofCalDigiVecOut->push_back(CbmTofDigi(*tDigi));
+ iNbCalDigis++;
+ }
+
+ int iHit0 = iNbHits; //starting index of current event
+ LOG(debug) << "Assign " << fTofHitsColl->GetEntriesFast() << " hits to event " << tEvent->GetNumber();
+ for (int32_t index = 0; index < fTofHitsColl->GetEntriesFast(); index++) {
+ CbmTofHit* pHit = (CbmTofHit*) fTofHitsColl->At(index);
+ new ((*fTofHitsCollOut)[iNbHits]) CbmTofHit(*pHit);
+ tEvent->AddData(ECbmDataType::kTofHit, iNbHits);
+
+ CbmMatch* pDigiMatch = (CbmMatch*) fTofDigiMatchColl->At(index);
+
+ // update content of match object to TS array
+ CbmMatch* mDigiMatch = new CbmMatch();
+ for (int32_t iLink = 0; iLink < pDigiMatch->GetNofLinks(); iLink++) {
+ CbmLink Link = pDigiMatch->GetLink(iLink);
+ Link.SetIndex(Link.GetIndex() + uDigi0);
+ mDigiMatch->AddLink(Link);
+ }
+
+ new ((*fTofDigiMatchCollOut)[iNbHits]) CbmMatch(*mDigiMatch);
+ delete mDigiMatch;
+ iNbHits++;
+ }
+
+ fiHitStart += iNbHits - iHit0;
+ fTofDigiVec.clear();
+ fTofCalDigiVec->clear();
+ fTofHitsColl->Delete(); //Clear("C");
+ fTofDigiMatchColl->Delete(); //Clear("C");
+ }
+
+ /// PAL: add TS statistics for monitoring and perf evaluation
+ timerTs.Stop();
+ LOG(debug) << GetName() << "::Exec: real time=" << timerTs.RealTime() << " CPU time=" << timerTs.CpuTime();
+ fProcessTime += timerTs.RealTime();
+ fuNbDigis += fDigiMan->GetNofDigis(ECbmModuleId::kTof) // TOF
+ + fDigiMan->GetNofDigis(ECbmModuleId::kT0); // T0
+ fuNbHits += fiHitStart;
+
+ std::stringstream logOut;
+ logOut << std::setw(20) << std::left << GetName() << " [";
+ logOut << std::fixed << std::setw(8) << std::setprecision(1) << std::right << timerTs.RealTime() * 1000. << " ms] ";
+ logOut << "TS " << iNbTs;
+ logOut << ", events " << fEventsColl->GetEntriesFast();
+ logOut << ", hits " << fiHitStart << ", time/1k-hit " << std::setprecision(4)
+ << timerTs.RealTime() * 1e6 / fiHitStart << " [ms]";
+ LOG(info) << logOut.str();
+ }
+ else {
+ // fTofDigisColl=fTofRawDigisColl;
+ // (VF) This does not work here. The digi manager does not foresee to add
+ // new data to the input array. So, I here copy the input digis into
+ // the array fTofDigisColl. Not very efficient, but temporary only, until
+ // also the internal data representations are changed to std::vectors.
+
+ fTofDigiVec.clear();
+ if (NULL != fT0DigiVec) { // 2022 data
+ for (int32_t iDigi = 0; iDigi < (int) (fT0DigiVec->size()); iDigi++) {
+ CbmTofDigi tDigi = fT0DigiVec->at(iDigi);
+ if (tDigi.GetType() != 5)
+ LOG(fatal) << "Wrong T0 type " << tDigi.GetType() << ", Addr 0x" << std::hex << tDigi.GetAddress();
+ if (tDigi.GetSide() == 1) { // HACK for May22 setup
+ tDigi.SetAddress(tDigi.GetSm(), tDigi.GetRpc(), tDigi.GetChannel() + 6, 0, tDigi.GetType());
+ }
+ fTofDigiVec.push_back(CbmTofDigi(tDigi));
+ }
+ }
+ for (int32_t iDigi = 0; iDigi < fDigiMan->GetNofDigis(ECbmModuleId::kTof); iDigi++) {
+ const CbmTofDigi* tDigi = fDigiMan->Get<CbmTofDigi>(iDigi);
+ fTofDigiVec.push_back(CbmTofDigi(*tDigi));
+
+ if (iDigi == 10000) break; // Only use 10000 digis for now. D.Smith
+ }
+ ExecEvent(option);
+ }
+}
+
+void CbmTaskTofClusterizer::ExecEvent(Option_t* /*option*/, CbmEvent* tEvent)
+{
+ // Clear output arrays
+ fTofCalDigiVec->clear();
+ fTofHitsColl->Delete(); // Computationally costly!, but hopefully safe
+ fTofDigiMatchColl->Delete();
+ FairRootFileSink* sink = (FairRootFileSink*) FairRootManager::Instance()->GetSink();
+ if (sink) { fiOutputTreeEntry = sink->GetOutTree()->GetEntries(); }
+
+ // Check for corruption (D.Smith: disabled )
+ // if (fTofDigiVec.size() > 20. * fDigiBdfPar->GetNbDet()) { // FIXME constant in code, skip this event
+ // LOG(warn) << "Skip processing of Tof DigiEvent with " << fTofDigiVec.size() << " digis ";
+ // return;
+ // }
+
+ fiNbHits = 0;
+
+ BuildClusters();
+}
+
+/************************************************************************************/
+void CbmTaskTofClusterizer::Finish()
+{
+ LOG(info) << "Finish with " << fdEvent << " processed events";
+
+ /// PAL: add run statistics for monitoring and perf evaluation
+ LOG(info) << "=====================================";
+ LOG(info) << GetName() << ": Finish run";
+ LOG(info) << GetName() << ": Run summary ";
+ LOG(info) << GetName() << ": Processing time : " << std::fixed << std::setprecision(3) << fProcessTime;
+ LOG(info) << GetName() << ": Nr of events : " << fdEvent;
+ LOG(info) << GetName() << ": Nr of input digis : " << fuNbDigis;
+ LOG(info) << GetName() << ": Nr of produced hits : " << fuNbHits;
+ LOG(info) << GetName() << ": Nr of hits / event : " << std::fixed << std::setprecision(2)
+ << (fdEvent > 0 ? fuNbHits / fdEvent : 0);
+ LOG(info) << GetName() << ": Nr of hits / digis : " << std::fixed << std::setprecision(2)
+ << (fuNbDigis > 0 ? fuNbHits / (double) fuNbDigis : 0);
+ LOG(info) << "=====================================";
+}
+
+void CbmTaskTofClusterizer::Finish(double calMode)
+{
+ SetCalMode(calMode);
+ Finish();
+}
+
+/************************************************************************************/
+// Functions common for all clusters approximations
+bool CbmTaskTofClusterizer::RegisterInputs()
+{
+ FairRootManager* fManager = FairRootManager::Instance();
+
+ if (NULL == fManager) {
+ LOG(error) << "CbmTaskTofClusterizer::RegisterInputs => Could not find "
+ "FairRootManager!!!";
+ return false;
+ }
+
+ fEventsColl = dynamic_cast<TClonesArray*>(fManager->GetObject("Event"));
+ if (NULL == fEventsColl) fEventsColl = dynamic_cast<TClonesArray*>(fManager->GetObject("CbmEvent"));
+
+ if (NULL == fEventsColl) LOG(info) << "CbmEvent not found in input file, assume eventwise input";
+
+ fDigiMan = CbmDigiManager::Instance();
+ fDigiMan->Init();
+ if (!fDigiMan->IsPresent(ECbmModuleId::kTof)) {
+ LOG(error) << GetName() << ": No Tof digi input!";
+ return false;
+ }
+ if (fDigiMan->IsPresent(ECbmModuleId::kT0)) {
+ LOG(info) << GetName() << ": separate T0 digi input!";
+ //fT0DigiVec = fManager->InitObjectAs<std::vector<CbmTofDigi> const*>("TzdDigi");
+ }
+ else {
+ LOG(info) << "No separate T0 digi input found.";
+ } // if( ! fT0DigiVec )
+
+ fTsHeader = fManager->InitObjectAs<CbmTsEventHeader const*>("EventHeader."); //for data
+ if (NULL == fTsHeader) { LOG(info) << "CbmTaskTofClusterizer::RegisterInputs => Could not get TsHeader Object"; }
+
+ return true;
+}
+bool CbmTaskTofClusterizer::RegisterOutputs()
+{
+ FairRootManager* rootMgr = FairRootManager::Instance();
+ rootMgr->InitSink();
+
+ fTofCalDigiVec = new std::vector<CbmTofDigi>();
+ fTofHitsColl = new TClonesArray("CbmTofHit", 100);
+ fTofDigiMatchColl = new TClonesArray("CbmMatch", 100);
+
+ if (NULL == fEventsColl) {
+ // Flag check to control whether digis are written in output root file
+ rootMgr->RegisterAny("TofCalDigi", fTofCalDigiVec, fbWriteDigisInOut);
+
+ // Flag check to control whether digis are written in output root file
+ rootMgr->Register("TofHit", "Tof", fTofHitsColl, fbWriteHitsInOut);
+ fTofHitsCollOut = fTofHitsColl;
+
+ rootMgr->Register("TofHitCalDigiMatch", "Tof", fTofDigiMatchColl, fbWriteHitsInOut);
+ LOG(info) << Form("Register DigiMatch in TS mode at %p with %d", fTofDigiMatchColl, (int) fbWriteHitsInOut);
+ fTofDigiMatchCollOut = fTofDigiMatchColl;
+ }
+ else { // CbmEvent - mode
+ fTofCalDigiVecOut = new std::vector<CbmTofDigi>();
+ fTofHitsCollOut = new TClonesArray("CbmTofHit", 10000);
+ fTofDigiMatchCollOut = new TClonesArray("CbmMatch", 10000);
+
+ rootMgr->RegisterAny("TofCalDigi", fTofCalDigiVecOut, fbWriteDigisInOut);
+ rootMgr->Register("TofHit", "Tof", fTofHitsCollOut, fbWriteHitsInOut);
+ rootMgr->Register("TofHitCalDigiMatch", "Tof", fTofDigiMatchCollOut, fbWriteHitsInOut);
+ LOG(info) << Form("Register DigiMatch in event mode at %p with %d ", fTofDigiMatchCollOut, (int) fbWriteHitsInOut);
+ }
+ return true;
+}
+bool CbmTaskTofClusterizer::InitParameters()
+{
+
+ // Initialize the TOF GeoHandler
+ bool isSimulation = false;
+ LOG(info) << "CbmTaskTofClusterizer::InitParameters - Geometry, Mapping, ... ??";
+
+ // Get Base Container
+ FairRun* ana = FairRun::Instance();
+ FairRuntimeDb* rtdb = ana->GetRuntimeDb();
+
+ CbmTofGeoHandler geoHandler;
+ int32_t iGeoVersion = geoHandler.Init(isSimulation);
+ if (k14a > iGeoVersion) {
+ LOG(error) << "CbmTaskTofClusterizer::InitParameters => Only compatible "
+ "with geometries after v14a !!!";
+ return false;
+ }
+ if (iGeoVersion == k14a) fTofId = new CbmTofDetectorId_v14a();
+ else
+ fTofId = new CbmTofDetectorId_v21a();
+
+ // 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) << "CbmTaskTofClusterizer::InitParameters => Could not obtain "
+ "the CbmTofDigiPar ";
+ return false;
+ }
+
+ fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
+ if (0 == fDigiBdfPar) {
+ LOG(error) << "CbmTaskTofClusterizer::InitParameters => Could not obtain "
+ "the CbmTofDigiBdfPar ";
+ return false;
+ }
+
+ rtdb->initContainers(ana->GetRunId());
+
+ LOG(info) << "CbmTaskTofClusterizer::InitParameter: currently " << fDigiPar->GetNrOfModules() << " digi cells ";
+
+ fdMaxTimeDist = fDigiBdfPar->GetMaxTimeDist(); // in ns
+ fdMaxSpaceDist = fDigiBdfPar->GetMaxDistAlongCh(); // in cm
+
+ if (fMaxTimeDist != fdMaxTimeDist) {
+ fdMaxTimeDist = fMaxTimeDist; // modify default
+ fdMaxSpaceDist = fdMaxTimeDist * fDigiBdfPar->GetSignalSpeed()
+ * 0.5; // cut consistently on positions (with default signal velocity)
+ }
+
+ LOG(info) << " BuildCluster with MaxTimeDist " << fdMaxTimeDist << ", MaxSpaceDist " << fdMaxSpaceDist;
+
+ fvDeadStrips.resize(fDigiBdfPar->GetNbDet());
+ return true;
+}
+/************************************************************************************/
+bool CbmTaskTofClusterizer::InitCalibParameter()
+{
+ // dimension and initialize calib parameter
+ int32_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
+
+ if (fTotMean != 0.) fdTTotMean = fTotMean; // adjust target mean for TOT
+
+ fvCPTOff.resize(iNbSmTypes);
+ fvCPTotGain.resize(iNbSmTypes);
+ fvCPTotOff.resize(iNbSmTypes);
+ fvCPWalk.resize(iNbSmTypes);
+ fvCPTOffY.resize(iNbSmTypes);
+ fvCPTOffYBinWidth.resize(iNbSmTypes);
+ fvCPTOffYRange.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);
+ fvCPTotOff[iSmType].resize(iNbSm * iNbRpc);
+ fvCPWalk[iSmType].resize(iNbSm * iNbRpc);
+ fvCPTOffY[iSmType].resize(iNbSm * iNbRpc);
+ fvCPTOffYBinWidth[iSmType].resize(iNbSm * iNbRpc);
+ fvCPTOffYRange[iSmType].resize(iNbSm * iNbRpc);
+ for (int32_t iSm = 0; iSm < iNbSm; iSm++) {
+ for (int32_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
+
+ fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc] = 0.; // initialize
+ fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc] = 0.; // initialize
+
+ /* D.Smith 23.8.23: For testing hit time calibration. Please remove when done.
+ fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc] = 1.; // initialize
+ fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc] = 200.; // initialize
+ fvCPTOffY[iSmType][iSm * iNbRpc + iRpc] = std::vector<double>(10000, 1000.);
+ for( size_t i = 0; i < fvCPTOffY[iSmType][iSm * iNbRpc + iRpc].size(); i++ )
+ {
+ fvCPTOffY[iSmType][iSm * iNbRpc + iRpc][i] = 1000.+500.*i;
+ }
+*/
+ int32_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
+ fvCPTOff[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
+ fvCPTotGain[iSmType][iSm * iNbRpc + iRpc].resize(iNbChan);
+ fvCPTotOff[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);
+ fvCPTotOff[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.; //initialize
+ fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 1.; //initialize
+ fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh][iSide] = 0.; //initialize
+ fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide].resize(nbClWalkBinX);
+ for (int32_t iWx = 0; iWx < nbClWalkBinX; iWx++) {
+ fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh][iSide][iWx] = 0.;
+ }
+ }
+ }
+ }
+ }
+ }
+ LOG(info) << "CbmTaskTofClusterizer::InitCalibParameter: defaults set";
+
+ if (fCalMode <= 0) { return true; } // Skip calibration from histograms in mode zero
+
+ /// Save old global file and folder pointer to avoid messing with FairRoot
+ // <= To prevent histos from being sucked in by the param file of the TRootManager!
+ TFile* oldFile = gFile;
+ TDirectory* oldDir = gDirectory;
+
+ LOG(info) << "CbmTaskTofClusterizer::InitCalibParameter: read histos from "
+ << "file " << fCalParFileName;
+
+ // read parameter from histos
+ if (fCalParFileName.IsNull()) return true;
+
+ TFile* calParFile = new TFile(fCalParFileName, "READ");
+ if (NULL == calParFile) {
+ if (fCalMode % 10 == 9) { //modify reference file name
+ int iCalMode = fCalParFileName.Index("tofClust") - 3;
+ fCalParFileName(iCalMode) = '3';
+ LOG(info) << "Modified CalFileName = " << fCalParFileName;
+ calParFile = new TFile(fCalParFileName, "update");
+ if (NULL == calParFile)
+ LOG(fatal) << "CbmTaskTofClusterizer::InitCalibParameter: "
+ << "file " << fCalParFileName << " does not exist!";
+
+ return true;
+ }
+ LOG(fatal) << "Calibration parameter file not existing!";
+ }
+
+ for (int32_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
+ int32_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
+ int32_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
+ TProfile* hSvel = (TProfile*) gDirectory->FindObjectAny(Form("cl_SmT%01d_Svel", iSmType));
+ for (int32_t iSm = 0; iSm < iNbSm; iSm++)
+ for (int32_t iRpc = 0; iRpc < iNbRpc; iRpc++) {
+
+ std::vector<std::vector<double>>& vCPTotGain = fvCPTotGain[iSmType][iSm * iNbRpc + iRpc];
+ std::vector<std::vector<double>>& vCPTOff = fvCPTOff[iSmType][iSm * iNbRpc + iRpc];
+ std::vector<std::vector<double>>& vCPTotOff = fvCPTotOff[iSmType][iSm * iNbRpc + iRpc];
+ std::vector<std::vector<std::vector<double>>>& vCPWalk = fvCPWalk[iSmType][iSm * iNbRpc + iRpc];
+
+ std::vector<double>& vCPTOffY = fvCPTOffY[iSmType][iSm * iNbRpc + iRpc];
+ double& vCPTOffYBinWidth = fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc];
+ double& vCPTOffYRange = fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc];
+
+ // update default parameter
+ if (NULL != hSvel) {
+ double Vscal = hSvel->GetBinContent(iSm * iNbRpc + iRpc + 1);
+ if (Vscal == 0.) Vscal = 1.;
+ Vscal *= fdModifySigvel; //1.03; // testing the effect of wrong signal velocity, FIXME
+ fDigiBdfPar->SetSigVel(iSmType, iSm, iRpc, fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc) * Vscal);
+ LOG(info) << "Modify " << iSmType << iSm << iRpc << " Svel by " << Vscal << " to "
+ << fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
+ }
+ TH2D* htempPos_pfx =
+ (TH2D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Pos_pfx", iSmType, iSm, iRpc));
+ TH2D* htempTOff_pfx =
+ (TH2D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_TOff_pfx", iSmType, iSm, iRpc));
+ TH1D* htempTot_Mean =
+ (TH1D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Mean", iSmType, iSm, iRpc));
+ TH1D* htempTot_Off =
+ (TH1D*) gDirectory->FindObjectAny(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Off", iSmType, iSm, iRpc));
+ if (NULL != htempPos_pfx && NULL != htempTOff_pfx && NULL != htempTot_Mean && NULL != htempTot_Off) {
+ int32_t iNbCh = fDigiBdfPar->GetNbChan(iSmType, iRpc);
+ //int32_t iNbinTot = htempTot_Mean->GetNbinsX();//not used any more
+ for (int32_t iCh = 0; iCh < iNbCh; iCh++) {
+ for (int32_t iSide = 0; iSide < 2; iSide++) {
+ double TotMean = htempTot_Mean->GetBinContent(iCh * 2 + 1 + iSide); //nh +1 empirical(?)
+ if (0.001 < TotMean) { vCPTotGain[iCh][iSide] *= fdTTotMean / TotMean; }
+ vCPTotOff[iCh][iSide] = htempTot_Off->GetBinContent(iCh * 2 + 1 + iSide);
+ }
+ double YMean = ((TProfile*) htempPos_pfx)->GetBinContent(iCh + 1); //nh +1 empirical(?)
+ double TMean = ((TProfile*) htempTOff_pfx)->GetBinContent(iCh + 1);
+ if (std::isnan(YMean) || std::isnan(TMean)) {
+ LOG(warn) << "Invalid calibration for TSRC " << iSmType << iSm << iRpc << iCh << ", use default!";
+ continue;
+ }
+ double dTYOff = YMean / fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
+ if (5 == iSmType || 8 == iSmType) dTYOff = 0.; // no valid Y positon for pad counters
+ vCPTOff[iCh][0] += -dTYOff + TMean;
+ vCPTOff[iCh][1] += +dTYOff + TMean;
+
+ if (5 == iSmType || 8 == iSmType) { // for PAD counters
+ vCPTOff[iCh][1] = vCPTOff[iCh][0];
+ vCPTotGain[iCh][1] = vCPTotGain[iCh][0];
+ vCPTotOff[iCh][1] = vCPTotOff[iCh][0];
+ }
+ //if(iSmType==0 && iSm==4 && iRpc==2 && iCh==26)
+ if (iSmType == 0 && iSm == 2 && iRpc == 0)
+ //if (iSmType == 9 && iSm == 0 && iRpc == 0 && iCh == 10) // select specific channel
+ LOG(info) << "InitCalibParameter:"
+ << " TSRC " << iSmType << iSm << iRpc << iCh
+ << Form(": YMean %6.3f, TYOff %6.3f, TMean %6.3f", YMean, dTYOff, TMean) << " -> "
+ << Form(" TOff %f, %f, TotG %f, %f ", vCPTOff[iCh][0], vCPTOff[iCh][1], vCPTotGain[iCh][0],
+ vCPTotGain[iCh][1]);
+
+ TH1D* htempWalk0 = (TH1D*) gDirectory->FindObjectAny(
+ Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh));
+ TH1D* htempWalk1 = (TH1D*) gDirectory->FindObjectAny(
+ Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh));
+ if (NULL == htempWalk0 && NULL == htempWalk1) { // regenerate Walk histos
+ int iSide = 0;
+ htempWalk0 = new TH1D(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh),
+ Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]",
+ iSmType, iSm, iRpc, iCh, iSide),
+ nbClWalkBinX, fdTOTMin, fdTOTMax);
+ iSide = 1;
+ htempWalk1 = new TH1D(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh),
+ Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]",
+ iSmType, iSm, iRpc, iCh, iSide),
+ nbClWalkBinX, fdTOTMin, fdTOTMax);
+ }
+ if (NULL != htempWalk0 && NULL != htempWalk1) { // reinitialize Walk array
+ LOG(debug) << "Initialize Walk correction for "
+ << Form(" SmT%01d_sm%03d_rpc%03d_Ch%03d", iSmType, iSm, iRpc, iCh);
+ if (htempWalk0->GetNbinsX() != nbClWalkBinX)
+ LOG(error) << "CbmTaskTofClusterizer::InitCalibParameter: "
+ "Inconsistent Walk histograms";
+ for (int32_t iBin = 0; iBin < nbClWalkBinX; iBin++) {
+ vCPWalk[iCh][0][iBin] = htempWalk0->GetBinContent(iBin + 1);
+ vCPWalk[iCh][1][iBin] = htempWalk1->GetBinContent(iBin + 1);
+ if (iSmType == 0 && iSm == 0 && iRpc == 2 && iCh == 15) // debugging
+ LOG(info) << Form("Read new SmT%01d_sm%03d_rpc%03d_Ch%03d bin %d cen %f walk %f %f", iSmType, iSm,
+ iRpc, iCh, iBin, htempWalk0->GetBinCenter(iBin + 1), vCPWalk[iCh][0][iBin],
+ vCPWalk[iCh][1][iBin]);
+ if (5 == iSmType || 8 == iSmType) { // Pad structure, enforce consistency
+ if (vCPWalk[iCh][1][iBin] != vCPWalk[iCh][0][iBin]) {
+ LOG(fatal) << "Inconsisten walk values for TSRC " << iSmType << iSm << iRpc << iCh << ", Bin "
+ << iBin << ": " << vCPWalk[iCh][0][iBin] << ", " << vCPWalk[iCh][1][iBin];
+ }
+ vCPWalk[iCh][1][iBin] = vCPWalk[iCh][0][iBin];
+ htempWalk1->SetBinContent(iBin + 1, vCPWalk[iCh][1][iBin]);
+ }
+ }
+ }
+ else {
+ LOG(info) << "No Walk histograms for TSRC " << iSmType << iSm << iRpc << iCh;
+ }
+ }
+ // look for TcorY corrections
+ LOG(info) << "Check for TCorY in " << gDirectory->GetName();
+ TH1* hTCorY =
+ (TH1*) gDirectory->FindObjectAny(Form("calXY_SmT%d_sm%03d_rpc%03d_TOff_z_all_TcorY", iSmType, iSm, iRpc));
+ if (NULL != hTCorY) {
+ vCPTOffYBinWidth = hTCorY->GetBinWidth(0);
+ vCPTOffYRange = hTCorY->GetXaxis()->GetXmax();
+ LOG(info) << "Initialize TCorY: TSR " << iSmType << iSm << iRpc << ", Bwid " << vCPTOffYBinWidth
+ << ", Range " << vCPTOffYRange;
+ vCPTOffY.resize(hTCorY->GetNbinsX());
+ for (int iB = 0; iB < hTCorY->GetNbinsX(); iB++) {
+ vCPTOffY[iB] = hTCorY->GetBinContent(iB + 1);
+ }
+ }
+ }
+ else {
+ LOG(warning) << " Calibration histos " << Form("cl_SmT%01d_sm%03d_rpc%03d_XXX", iSmType, iSm, iRpc)
+ << " not found. ";
+ }
+ }
+ }
+
+ /// Restore old global file and folder pointer to avoid messing with FairRoot
+ // <= To prevent histos from being sucked in by the param file of the TRootManager!
+ gFile = oldFile;
+ gDirectory = oldDir;
+ LOG(info) << "CbmTaskTofClusterizer::InitCalibParameter: initialization done";
+ return true;
+}
+
+bool CbmTaskTofClusterizer::InitAlgos()
+{
+ // Needed as external TOT values might be different than ones used for input histo reading (TO DO: FIX)
+ if (fTotMax != 0.) fdTOTMax = fTotMax;
+ if (fTotMin != 0.) fdTOTMin = fTotMin;
+ LOG(info) << "ToT init to Min " << fdTOTMin << " Max " << fdTOTMax;
+
+ /// Go to Top volume of the geometry in the GeoManager to make sure our nodes are found
+ gGeoManager->CdTop();
+
+ int32_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes();
+
+ //Prepare storage vectors
+ fStorDigi.resize(iNbSmTypes);
+ for (int32_t iSmType = 0; iSmType < iNbSmTypes; iSmType++) {
+ int32_t iNbSm = fDigiBdfPar->GetNbSm(iSmType);
+ int32_t iNbRpc = fDigiBdfPar->GetNbRpc(iSmType);
+ fStorDigi[iSmType].resize(iNbSm * iNbRpc);
+ }
+
+ // Create map with unique detector IDs and fill (needed only for dead strip array)
+ std::map<uint32_t, uint32_t> detIdIndexMap;
+ for (int32_t ind = 0; ind < fDigiBdfPar->GetNbDet(); ind++) {
+ int32_t iUniqueId = fDigiBdfPar->GetDetUId(ind);
+ detIdIndexMap[iUniqueId] = ind;
+ }
+
+ // 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<cbm::algo::ClusterizerTofRpcPar> par(new cbm::algo::ClusterizerTofRpcPar());
+
+ const int32_t iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
+ const int32_t iDetIndx = detIdIndexMap[iDetId];
+ par->fDeadStrips = fvDeadStrips[iDetIndx];
+ par->fPosYMaxScal = fPosYMaxScal;
+ par->fdMaxTimeDist = fdMaxTimeDist;
+ par->fdMaxSpaceDist = fdMaxSpaceDist;
+ par->fCPTOffYBinWidth = fvCPTOffYBinWidth[iSmType][iSm * iNbRpc + iRpc];
+ par->fCPTOffY = fvCPTOffY[iSmType][iSm * iNbRpc + iRpc];
+ par->fCPTOffYRange = fvCPTOffYRange[iSmType][iSm * iNbRpc + iRpc];
+ par->fSigVel = fDigiBdfPar->GetSigVel(iSmType, iSm, iRpc);
+ par->fTimeRes = 0.08;
+ par->numClWalkBinX = nbClWalkBinX;
+ par->TOTMax = fdTOTMax;
+ par->TOTMin = fdTOTMin;
+ par->swapChannelSides = fbSwapChannelSides;
+ par->channelDeadtime = fdChannelDeadtime;
+
+ int32_t iNbChan = fDigiBdfPar->GetNbChan(iSmType, iRpc);
+ par->fChanPar.resize(iNbChan);
+ for (int32_t iCh = 0; iCh < iNbChan; iCh++) {
+
+ const int32_t address = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, 0, iSmType);
+ CbmTofCell* channelInfo = fDigiPar->GetCell(address);
+ if (channelInfo == nullptr) { continue; } //D.Smith 17.8.23: Needed as channel info sometimes not defined
+
+ gGeoManager->FindNode(channelInfo->GetX(), channelInfo->GetY(), channelInfo->GetZ());
+ const double* tra_ptr = gGeoManager->MakePhysicalNode()->GetMatrix()->GetTranslation();
+
+ //init Tof cell
+ cbm::algo::TofCell& cell = par->fChanPar[iCh].cell;
+ cell.pos = ROOT::Math::XYZVector(tra_ptr[0], tra_ptr[1], tra_ptr[2]);
+
+ /* why doesn't this variant work?
+ 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 = fvCPTOff[iSmType][iSm * iNbRpc + iRpc][iCh];
+ par->fChanPar[iCh].fvCPTotGain = fvCPTotGain[iSmType][iSm * iNbRpc + iRpc][iCh];
+ par->fChanPar[iCh].fvCPTotOff = fvCPTotOff[iSmType][iSm * iNbRpc + iRpc][iCh];
+ par->fChanPar[iCh].fvCPWalk = fvCPWalk[iSmType][iSm * iNbRpc + iRpc][iCh];
+ par->fChanPar[iCh].address = address;
+ }
+ fAlgo[iSmType][iSm * iNbRpc + iRpc].SetParams(std::move(par));
+ }
+ }
+ }
+ return true;
+}
+
+
+/************************************************************************************/
+bool CbmTaskTofClusterizer::BuildClusters()
+{
+ gGeoManager->CdTop();
+
+ if (fTofDigiVec.empty()) {
+ LOG(info) << " No RawDigis defined ! Check! ";
+ return false;
+ }
+ LOG(info) << "Build clusters from " << fTofDigiVec.size() << " digis in event " << fdEvent + 1;
+
+ int32_t iNbTofDigi = fTofDigiVec.size();
+
+ if (bAddBeamCounterSideDigi) {
+ // Duplicate type "5" - digis
+ for (int32_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
+ CbmTofDigi* pDigi = &(fTofDigiVec.at(iDigInd));
+ if (pDigi->GetType() == 5) { // || pDigi->GetType() == 8) {
+ if (pDigi->GetSide() == 1) {
+ bAddBeamCounterSideDigi = false; // disable for current data set
+ LOG(info) << "Start counter digi duplication disabled";
+ break;
+ }
+ fTofDigiVec.push_back(CbmTofDigi(*pDigi));
+ CbmTofDigi* pDigiN = &(fTofDigiVec.back());
+ pDigiN->SetAddress(pDigi->GetSm(), pDigi->GetRpc(), pDigi->GetChannel(), (0 == pDigi->GetSide()) ? 1 : 0,
+ pDigi->GetType());
+ LOG(debug) << "Duplicated digi " << fTofDigiVec.size() << " with address 0x" << std::hex
+ << pDigiN->GetAddress();
+ }
+ }
+ iNbTofDigi = fTofDigiVec.size();
+ }
+
+ //D.Smith 10.8.23: This entire if(true){...} block doesn't seem to actually do anything
+ //for run 2391 data if enabled. Check with PA whether it can be dropped (might be obsolete).
+ //The original purpose was to correct missing hits, which might have been fixed.
+ if (true) {
+ for (int32_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++) {
+ CbmTofDigi* pDigi = &(fTofDigiVec.at(iDigInd));
+ const double SmType = pDigi->GetType();
+ const double Sm = pDigi->GetSm();
+ const double Rpc = pDigi->GetRpc();
+ const double Chan = pDigi->GetChannel();
+ const double Side = pDigi->GetSide();
+ int32_t iDetIndx = fDigiBdfPar->GetDetInd(pDigi->GetAddress() & DetMask);
+
+ if (fDigiBdfPar->GetNbDet() - 1 < iDetIndx || iDetIndx < 0) {
+ LOG(debug) << Form(" Wrong DetIndx %d >< %d ", iDetIndx, fDigiBdfPar->GetNbDet());
+ break;
+ }
+ CbmTofDigi* pDigi2Min = NULL;
+ double dTDifMin = dDoubleMax;
+
+ if (fDutId > -1) {
+ for (int32_t iDigI2 = 0; iDigI2 < iNbTofDigi; iDigI2++) {
+ CbmTofDigi* pDigi2 = &(fTofDigiVec.at(iDigI2));
+ const double SmType2 = pDigi2->GetType();
+ const double Sm2 = pDigi2->GetSm();
+ const double Rpc2 = pDigi2->GetRpc();
+ const double Chan2 = pDigi2->GetChannel();
+ const double Side2 = pDigi2->GetSide();
+
+ if (iDetIndx == fDigiBdfPar->GetDetInd(pDigi2->GetAddress())) {
+ if (Side != Side2) {
+ if (Chan == Chan2) {
+ double dTDif = TMath::Abs(pDigi->GetTime() - pDigi2->GetTime());
+ if (dTDif < dTDifMin) {
+ dTDifMin = dTDif;
+ pDigi2Min = pDigi2;
+ }
+ }
+ else if (TMath::Abs(Chan - Chan2)
+ == 1) { // opposite side missing, neighbouring channel has hit on opposite side // FIXME
+ // check that same side digi of neighbouring channel is absent
+ int32_t iDigI3 = 0;
+ for (; iDigI3 < iNbTofDigi; iDigI3++) {
+ CbmTofDigi* pDigi3 = &(fTofDigiVec.at(iDigI3));
+ if (pDigi3->GetSide() == Side && Chan2 == pDigi3->GetChannel()) break;
+ }
+ if (iDigI3 == iNbTofDigi) // same side neighbour did not fire
+ {
+ int32_t iCorMode = 0; // Missing hit correction mode
+ switch (iCorMode) {
+ case 0: // no action
+ break;
+ case 1: // shift found hit
+ if (Side == 0) pDigi2->SetAddress(Sm, Rpc, Chan, 1 - Side, SmType);
+ else
+ pDigi->SetAddress(Sm2, Rpc2, Chan2, 1 - Side2, SmType2);
+
+ break;
+ case 2: // insert missing hits
+ fTofDigiVec.push_back(CbmTofDigi(*pDigi));
+ CbmTofDigi* pDigiN = &(fTofDigiVec.back());
+ pDigiN->SetAddress(Sm, Rpc, Chan2, Side, SmType);
+ pDigiN->SetTot(pDigi2->GetTot());
+
+ fTofDigiVec.push_back(CbmTofDigi(*pDigi2));
+ CbmTofDigi* pDigiN2 = &(fTofDigiVec.back());
+ pDigiN2->SetAddress(Sm2, Rpc2, Chan, Side2, SmType2);
+ pDigiN2->SetTot(pDigi->GetTot());
+
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ if (pDigi2Min != NULL) {
+ CbmTofDetectorInfo xDetInfo(ECbmModuleId::kTof, SmType, Sm, Rpc, 0, Chan);
+ int32_t iChId = fTofId->SetDetectorInfo(xDetInfo);
+ CbmTofCell* cell = fDigiPar->GetCell(iChId);
+ if (NULL == cell) {
+ LOG(warning) << Form("BuildClusters: invalid ChannelInfo for 0x%08x", iChId);
+ continue;
+ }
+ if (fDigiBdfPar->GetSigVel(SmType, Sm, Rpc) * dTDifMin * 0.5 < fPosYMaxScal * cell->GetSizey()) {
+ //check consistency
+ if (8 == SmType || 5 == SmType) {
+ if (pDigi->GetTime() != pDigi2Min->GetTime()) {
+ LOG(warning) << " BuildClusters: Inconsistent duplicated digis in event " << fdEvent + 1
+ << ", Ind: " << iDigInd; // << "CTyp: " << pDigi->GetCounterType;
+ LOG(warning) << " " << pDigi->ToString();
+ LOG(warning) << " " << pDigi2Min->ToString();
+
+ pDigi2Min->SetTot(pDigi->GetTot());
+ pDigi2Min->SetTime(pDigi->GetTime());
+ }
+ }
+ }
+ }
+ }
+ } // true end
+
+ // Calibrate RawDigis
+ *fTofCalDigiVec = CalibRawDigis(fTofDigiVec);
+ // *fTofCalDigiVec = fTofDigiVec; //disable calibration for now (done in algo). Linked indices will differ!
+
+ // Then loop over the digis array and store the Digis in separate vectors for
+ // each RPC modules
+ for (int32_t iDigInd = 0; iDigInd < fTofCalDigiVec->size(); iDigInd++) {
+ CbmTofDigi* pDigi = &(fTofCalDigiVec->at(iDigInd));
+
+ // These are doubles in the digi class
+ const double SmType = pDigi->GetType();
+ const double Sm = pDigi->GetSm();
+ const double Rpc = pDigi->GetRpc();
+ const int NbRpc = fDigiBdfPar->GetNbRpc(SmType);
+ fStorDigi[SmType][Sm * NbRpc + Rpc].push_back(std::make_pair(pDigi, iDigInd));
+ }
+ // inspect digi array (D.Smith: Block removed for cleanup. Could be in monitoring class in principle.)
+
+ for (uint32_t iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); 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<std::pair<CbmTofDigi*, int32_t>>& digiExp = fStorDigi[iSmType][rpcIdx];
+
+ //call cluster finder
+ std::vector<cbm::algo::TofCluster> clusters = fAlgo[iSmType][rpcIdx](digiExp);
+
+ //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, fiNbHits, cluster.weightedTime, cluster.weightedTimeErr,
+ cluster.vDigiIndRef.size(), int32_t(cluster.weightsSum * 10.));
+
+ fiNbHits++;
+ //if (event) event->AddData(ECbmDataType::kTofHit, hitIndex);
+
+ CbmMatch* digiMatch = new ((*fTofDigiMatchColl)[hitIndex]) CbmMatch();
+ for (uint32_t i = 0; i < cluster.vDigiIndRef.size(); i++) {
+ double tot = fTofCalDigiVec->at(cluster.vDigiIndRef.at(i)).GetTot();
+ digiMatch->AddLink(CbmLink(tot, cluster.vDigiIndRef.at(i), fiOutputTreeEntry, fiFileIndex));
+ }
+ }
+ // digiExp.clear();
+ // digiInd.clear();
+ }
+ }
+ }
+ std::cout << "hits " << fiNbHits << std::endl;
+
+ fdEvent++;
+ return true;
+}
+
+/* No longer used. Keeping for reference / discussion of comments
+void CbmTaskTofClusterizer::StoreHit(TofCluster& cluster, int32_t iSmType, int32_t iChId, int32_t iNbCh, int32_t iSm,
+ int32_t iRpc)
+{
+ double* tra_ptr = fCellIdGeoMap[iChId]->GetMatrix()->GetTranslation();
+ double* rot_ptr = fCellIdGeoMap[iChId]->GetMatrix()->GetRotationMatrix();
+
+ TofCell cell;
+ cell.pos = ROOT::Math::XYZVector(tra_ptr[0], tra_ptr[1], tra_ptr[2]);
+ cell.rotation = ROOT::Math::Rotation3D(&rot_ptr[0], &rot_ptr[9]);
+ cell.sizeX = fDigiPar->GetCell(iChId)->GetSizex();
+
+ //D.Smith 10.8.23: Replacing this with the lines above seems to work, although the exact
+ //cell might be a different one when called from AddNextChan().
+ //int32_t iChm = floor(cluster.weightedPos.X() / fChannelInfo->GetSizex()) + iNbCh / 2;
+
+ //D.Smith 15.8.23: Can we drop the condition?
+ if (5 != iSmType) { cluster.finalize(cell, iNbCh, ClusterizerTofRpcPar()); }
+
+ TVector3 hitPos(cluster.globalPos.X(), cluster.globalPos.Y(), cluster.globalPos.Z());
+ TVector3 hitPosErr(cluster.globalErr.X(), cluster.globalErr.Y(), cluster.globalErr.Z());
+
+ const int32_t iDetId = CbmTofAddress::GetUniqueAddress(iSm, iRpc, cluster.channel, 0, iSmType);
+
+ if (cluster.vDigiIndRef.size() < 2) {
+ LOG(warning) << "Digi refs for Hit " << fiNbHits << ": vDigiIndRef.size()";
+ }
+ CbmTofHit* pHit = new CbmTofHit(iDetId, hitPos,
+ hitPosErr, //local detector coordinates
+ fiNbHits, // this number is used as reference!!
+ cluster.weightedTime,
+ cluster.vDigiIndRef.size(), // number of linked digis = 2*CluSize
+ int32_t(cluster.weightsSum * 10.)); //channel -> Tot
+ pHit->SetTimeError(dTimeRes);
+
+ // output hit
+ new ((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit);
+ pHit->Delete();
+
+ CbmMatch* digiMatch = new ((*fTofDigiMatchColl)[fiNbHits]) CbmMatch();
+
+ for (size_t i = 0; i < cluster.vDigiIndRef.size(); i++) {
+ double dTot = fTofCalDigiVec->at(cluster.vDigiIndRef.at(i)).GetTot();
+ digiMatch->AddLink(CbmLink(dTot, cluster.vDigiIndRef.at(i), fiOutputTreeEntry, fiFileIndex));
+ }
+ fiNbHits++;
+}
+*/
+
+std::vector<CbmTofDigi> CbmTaskTofClusterizer::CalibRawDigis(const std::vector<CbmTofDigi>& digiVec)
+{
+ // channel deadtime map
+ std::map<int32_t, double> mChannelDeadTime;
+ std::vector<CbmTofDigi> calDigiVecOut;
+
+ for (int32_t iDigi = 0; iDigi < digiVec.size(); iDigi++) {
+
+ CbmTofDigi pDigi = digiVec.at(iDigi);
+ const double SmType = pDigi.GetType();
+ const double Sm = pDigi.GetSm();
+ const double Rpc = pDigi.GetRpc();
+ const double Chan = pDigi.GetChannel();
+ const double Side = pDigi.GetSide();
+ const int NbRpc = fDigiBdfPar->GetNbRpc(SmType);
+ const uint32_t rpcIdx = Sm * NbRpc + Rpc;
+
+ if (fbSwapChannelSides && 5 != SmType && 8 != SmType) {
+ pDigi.SetAddress(Sm, Rpc, Chan, (0 == Side) ? 1 : 0, SmType);
+ }
+
+ // Check dead time
+ const int32_t iAddr = pDigi.GetAddress();
+ auto it = mChannelDeadTime.find(iAddr);
+ if (it != mChannelDeadTime.end() && pDigi.GetTime() <= it->second) {
+ it->second = pDigi.GetTime() + fdChannelDeadtime;
+ continue;
+ }
+ mChannelDeadTime[iAddr] = pDigi.GetTime() + fdChannelDeadtime;
+
+ // Create calibrated digi
+ CbmTofDigi pCalDigi(pDigi);
+
+ // calibrate Digi Time
+ pCalDigi.SetTime(pCalDigi.GetTime() - fvCPTOff[SmType][rpcIdx][Chan][Side]);
+
+ // subtract Offset
+ double dTot = pCalDigi.GetTot() - fvCPTotOff[SmType][rpcIdx][Chan][Side];
+ if (dTot < 0.001) dTot = 0.001;
+
+ // calibrate Digi ToT
+ pCalDigi.SetTot(dTot * fvCPTotGain[SmType][rpcIdx][Chan][Side]);
+
+ // walk correction
+ std::vector<double>& walk = fvCPWalk[SmType][rpcIdx][Chan][Side];
+ const double dTotBinSize = (fdTOTMax - fdTOTMin) / nbClWalkBinX;
+ int32_t iWx = (int32_t)((pCalDigi.GetTot() - fdTOTMin) / dTotBinSize);
+
+ if (0 > iWx) { iWx = 0; }
+ if (iWx >= nbClWalkBinX) { iWx = nbClWalkBinX - 1; }
+
+ const double dDTot = (pCalDigi.GetTot() - fdTOTMin) / dTotBinSize - (double) iWx - 0.5;
+ double dWT = walk[iWx];
+
+ // linear interpolation to next bin //memory leak??? (D.Smith 10.8.23: Clarify this comment!)
+ if (dDTot > 0) {
+ if (iWx < nbClWalkBinX - 1) { dWT += dDTot * (walk[iWx + 1] - walk[iWx]); }
+ }
+ else {
+ if (0 < iWx) { dWT -= dDTot * (walk[iWx - 1] - walk[iWx]); }
+ }
+ pCalDigi.SetTime(pCalDigi.GetTime() - dWT); // calibrate Digi Time
+ calDigiVecOut.push_back(pCalDigi);
+ }
+
+ // D.Smith 10.8.23: Sorting might not be needed. Check with P.A.
+ LOG(debug) << "CbmTaskTofClusterizer::BuildClusters: Sort " << calDigiVecOut.size() << " calibrated digis ";
+ /// Sort the buffers of hits due to the time offsets applied
+ std::sort(calDigiVecOut.begin(), calDigiVecOut.end(),
+ [](const CbmTofDigi& a, const CbmTofDigi& b) -> bool { return a.GetTime() < b.GetTime(); });
+
+ return calDigiVecOut;
+}
+
+void CbmTaskTofClusterizer::SetDeadStrips(int32_t iDet, uint32_t ival)
+{
+ if (fvDeadStrips.size() < static_cast<size_t>(iDet + 1)) fvDeadStrips.resize(iDet + 1);
+ fvDeadStrips[iDet] = ival;
+}
diff --git a/reco/tasks/CbmTaskTofClusterizer.h b/reco/tasks/CbmTaskTofClusterizer.h
new file mode 100644
index 0000000000000000000000000000000000000000..d5e9df9cd4d064e17fee69ca48357140f782b953
--- /dev/null
+++ b/reco/tasks/CbmTaskTofClusterizer.h
@@ -0,0 +1,227 @@
+/* Copyright (C) 2023 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 */
+
+#ifndef CBMTASKTOFCLUSTERIZER_H
+#define CBMTASKTOFCLUSTERIZER_H 1
+
+// TOF Classes and includes
+// Input/Output
+class CbmEvent;
+// Geometry
+class CbmTofDetectorId;
+class CbmTofDigiPar;
+class CbmTofDigiBdfPar;
+class CbmTofCell;
+class CbmDigiManager;
+class CbmTsEventHeader;
+
+#include "CbmMatch.h"
+#include "CbmTofDigi.h"
+
+#include "tof/ClusterizerTof.h"
+
+// ROOT Classes and includes
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+// FAIR classes and includes
+#include "FairTask.h"
+
+// ROOT Classes and includes
+class TClonesArray;
+
+// C++ Classes and includes
+#include <vector>
+
+class CbmTaskTofClusterizer : public FairTask {
+
+public:
+ /**
+ ** @brief Constructor.
+ **/
+ CbmTaskTofClusterizer();
+
+ /**
+ ** @brief Constructor.
+ **/
+ CbmTaskTofClusterizer(const char* name, int32_t verbose = 1, bool writeDataInOut = true);
+ /**
+ ** @brief Destructor.
+ **/
+ virtual ~CbmTaskTofClusterizer();
+
+ /**
+ ** @brief Inherited from FairTask.
+ **/
+ virtual InitStatus Init();
+
+ /**
+ ** @brief Inherited from FairTask.
+ **/
+ virtual void SetParContainers();
+
+ /**
+ ** @brief Inherited from FairTask.
+ **/
+ virtual void Exec(Option_t* option);
+ virtual void ExecEvent(Option_t* option, CbmEvent* tEvent = NULL);
+
+ /**
+ ** @brief Inherited from FairTask.
+ **/
+ virtual void Finish();
+ virtual void Finish(double calMode);
+
+ inline void SetCalMode(int32_t iMode) { fCalMode = iMode; }
+ inline void SetDutId(int32_t Id) { fDutId = Id; }
+
+ inline void PosYMaxScal(double val) { fPosYMaxScal = val; }
+ inline void SetTotMax(double val) { fTotMax = val; }
+ inline void SetTotMin(double val) { fTotMin = val; }
+ inline void SetTotMean(double val) { fTotMean = val; }
+ inline void SetMaxTimeDist(double val) { fMaxTimeDist = val; }
+ inline void SetChannelDeadtime(double val) { fdChannelDeadtime = val; }
+ inline void SetCalParFileName(TString CalParFileName) { fCalParFileName = CalParFileName; }
+
+ inline int GetNbHits() { return fiNbHits; }
+ inline double GetTotMean() { return fTotMean; }
+
+ void SwapChannelSides(bool bSwap) { fbSwapChannelSides = bSwap; }
+ void SetFileIndex(int32_t iIndex) { fiFileIndex = iIndex; }
+ void SetWriteDigisInOut(bool bDigis) { fbWriteDigisInOut = bDigis; }
+ void SetWriteHitsInOut(bool bHits) { fbWriteHitsInOut = bHits; }
+ void SetDeadStrips(int32_t iDet, uint32_t ival);
+
+protected:
+private:
+ int32_t iNbTs = 0;
+ int fiHitStart = 0;
+ bool bAddBeamCounterSideDigi = true;
+ const double dDoubleMax = 1.E300;
+ const int32_t nbClWalkBinX = 50; // was 100 (11.10.2018)
+ const int32_t nbClWalkBinY = 41; // choose odd number to have central bin symmetric around 0
+ const int32_t DetMask = 0x001FFFFF; // geo v21a
+
+ std::vector<CbmTofDigi>* fT0DigiVec = nullptr; //! T0 Digis
+
+
+ /**
+ ** @brief Copy constructor.
+ **/
+ CbmTaskTofClusterizer(const CbmTaskTofClusterizer&);
+ /**
+ ** @brief Copy operator.
+ **/
+ CbmTaskTofClusterizer& operator=(const CbmTaskTofClusterizer&);
+
+ // 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 Apply calibration to digis
+ **/
+ std::vector<CbmTofDigi> CalibRawDigis(const std::vector<CbmTofDigi>& digiVec);
+
+ /**
+ ** @brief Build clusters out of ToF Digis and store the resulting info in a TofHit.
+ **/
+ bool BuildClusters();
+
+ /**
+ ** @brief Create one algo object for each RPC
+ **/
+ bool InitAlgos();
+
+ // Hit finder algo
+ std::map<uint32_t, std::map<uint32_t, cbm::algo::ClusterizerTof>> fAlgo = {}; //[nbType][nbSm*nbRpc]
+
+ // ToF geometry variables
+ CbmTofDetectorId* fTofId;
+ CbmTofDigiPar* fDigiPar;
+ CbmTofDigiBdfPar* fDigiBdfPar;
+
+ const CbmTsEventHeader* fTsHeader;
+
+ // Input variables
+ std::vector<CbmTofDigi> fTofDigiVec {}; //! TOF Digis
+ CbmDigiManager* fDigiMan; // TOF Input Digis
+ TClonesArray* fEventsColl; // CBMEvents (time based)
+
+ // Output variables
+ bool fbWriteHitsInOut;
+ bool fbWriteDigisInOut;
+ std::vector<CbmTofDigi>* fTofCalDigiVec = nullptr; //! // Calibrated TOF Digis
+ TClonesArray* fTofHitsColl; // TOF hits
+ TClonesArray* fTofDigiMatchColl; // TOF Digi Links
+ std::vector<CbmTofDigi>* fTofCalDigiVecOut = nullptr; //! // Calibrated TOF Digis
+ TClonesArray* fTofHitsCollOut; // TOF hits
+ TClonesArray* fTofDigiMatchCollOut; // TOF Digi Links
+ int32_t fiNbHits; // Index of the CbmTofHit TClonesArray
+
+ // Intermediate storage variables (digi, index)
+ std::vector<std::vector<std::vector<std::pair<CbmTofDigi*, int32_t>>>> fStorDigi; //[nbType][nbSm*nbRpc][nDigis]
+
+ //Calibration parameters
+ 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<double>>>> fvCPTotOff; //[nSMT][nRpc][nCh][nbSide]
+ std::vector<std::vector<std::vector<std::vector<std::vector<double>>>>>
+ fvCPWalk; //[nSMT][nRpc][nCh][nbSide][nbWalkBins]
+ std::vector<std::vector<std::vector<double>>> fvCPTOffY; //[nSMT][nRpc][nBin]
+ std::vector<std::vector<double>> fvCPTOffYBinWidth; //[nSMT][nRpc]
+ std::vector<std::vector<double>> fvCPTOffYRange; //[nSMT][nRpc]
+
+ std::vector<uint32_t> fvDeadStrips; //[nbDet]
+
+ // Calib
+ int32_t fCalMode;
+ int32_t fDutId;
+
+ double fPosYMaxScal;
+ double fTotMax;
+ double fTotMin;
+ double fTotOff;
+ double fTotMean;
+ double fMaxTimeDist;
+ double fdChannelDeadtime;
+
+ TString fCalParFileName; // name of the file name with Calibration Parameters
+
+ double fdTOTMax;
+ double fdTOTMin;
+ double fdTTotMean;
+
+ double fdMaxTimeDist; // Isn't this just a local variable? Why make it global and preset?!?
+ double fdMaxSpaceDist; // Isn't this just a local variable? Why make it global and preset?!?
+ double fdModifySigvel;
+
+ double fdEvent;
+
+ double fProcessTime = 0.0;
+ uint64_t fuNbDigis = 0;
+ uint64_t fuNbHits = 0;
+
+ bool fbSwapChannelSides;
+ int32_t fiOutputTreeEntry;
+ int32_t fiFileIndex;
+
+ ClassDef(CbmTaskTofClusterizer, 1);
+};
+
+#endif // CBMTASKTOFCLUSTERIZER_H