diff --git a/algo/detectors/tof/Clusterizer.cxx b/algo/detectors/tof/Clusterizer.cxx
index 68b462701d1eea61129cbbaf53f147fc689e19aa..9b23eb627ed0c9857b30873af38c458d0790b926 100644
--- a/algo/detectors/tof/Clusterizer.cxx
+++ b/algo/detectors/tof/Clusterizer.cxx
@@ -20,7 +20,8 @@ namespace cbm::algo::tof
{
//std::vector<inputType> input = calibrateDigis(digisIn);
std::vector<inputType> input = chanSortDigis(digisIn);
- return buildClusters(input);
+ //return buildClusters(input);
+ return buildClustersIter(input); //Iterator based implementation
}
std::vector<Clusterizer::inputType> Clusterizer::chanSortDigis(inputType& digisIn)
@@ -330,119 +331,116 @@ namespace cbm::algo::tof
}
- /*
- Clusterizer::resultType Clusterizer::buildClusters(inputType& input)
+ Clusterizer::resultType Clusterizer::buildClustersIter(std::vector<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;
+ TofCell* cell = nullptr;
// Last Channel Temp variables
int32_t lastChan = -1;
double lastPosY = 0.0;
double lastTime = 0.0;
- //reset counter
- numSameSide = 0;
+ const size_t numChan = fParams.fChanPar.size();
- for (int32_t chan = 0; chan < fParams.fChanPar.size(); chan++) {
+ for (int32_t chan = 0; chan < numChan; chan++) {
+ if (fParams.fDeadStrips & (1 << chan)) { continue; } // skip over dead channels
- std::vector<CbmTofDigi*>& storDigiExp = digisExp[chan];
- std::vector<int32_t>& storDigiInd = digisInd[chan];
- ClusterizerChanPar& chanPar = fParams.fChanPar[chan];
+ inputType& storDigi = input[chan];
+ auto digiIt = storDigi.begin();
- auto digiExpIt = storDigiExp.begin();
- auto digiIndIt = storDigiInd.begin();
+ while (1 < std::distance(digiIt, storDigi.end())) {
- 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;
+ while (digiIt->first->GetSide() == std::next(digiIt, 1)->first->GetSide()) { // Not one Digi of each end!
+ digiIt++;
+ if (2 > std::distance(digiIt, storDigi.end())) { break; }
}
- if (2 > std::distance(digiExpIt, storDigiExp.end())) break; // 2 Digis = both sides present
+ if (2 > std::distance(digiIt, storDigi.end())) { break; }
- TofCell* channelInfo = &chanPar.cell;
- CbmTofDigi* xDigiA = *digiExpIt;
- CbmTofDigi* xDigiB = *std::next(digiExpIt, 1);
+ // 2 Digis = both sides present
+ cell = &fParams.fChanPar[chan].cell;
+ CbmTofDigi* xDigiA = digiIt->first;
+ CbmTofDigi* xDigiB = std::next(digiIt, 1)->first;
- // The "Strip" time is the mean time between each end
- const double time = 0.5 * (xDigiA->GetTime() + xDigiB->GetTime());
+ // 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.);
- // Weight is the total charge => sum of both ends ToT
- const double totSum = xDigiA->GetTot() + xDigiB->GetTot();
+ double timeDif = xDigiA->GetTime() - xDigiB->GetTime();
- if (nullptr == trafoCell) {
- trafoCell = channelInfo;
- iTrafoCell = chan;
- }
+ 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
- // 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.);
+ while (std::distance(digiIt, storDigi.end()) > 2 && std::abs(pos.Y()) > cell->sizeY * fParams.fPosYMaxScal) {
- 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);
+ CbmTofDigi* xDigiC = std::next(digiIt, 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; }
+ else {
+ xDigiB = xDigiC;
+ }
+ digiIt++;
}
- 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++;
+ if (std::abs(pos.Y()) > cell->sizeY * fParams.fPosYMaxScal) { // remove both digis
+ digiIt++;
continue;
- } // Pair leads to hit oustide of strip limits
+ }
+ // 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.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);
+ 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();
}
}
- cluster.add(pos, time, totSum, totSum, *digiIndIt, *std::next(digiIndIt, 1));
- digiExpIt += 2;
- digiIndIt += 2;
+ cluster.add(pos, time, totSum, totSum, storDigi[0].second, storDigi[1].second);
+ digiIt += 2;
lastChan = chan;
lastPosY = pos.Y();
lastTime = time;
- } // while (1 < storDigiExp.size())
- } // for( int32_t chan = 0; chan < iNbCh; chan++ )
+ 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++ )
- // Save last cluster if it exists
+ // Now check if another hit/cluster is started
+ // and save it if it's the case.
if (0 < cluster.numChan()) {
- cluster.normalize(fParams.SysTimeRes);
- cluster.finalize(*trafoCell, iTrafoCell, fParams);
+ cluster.normalize(fParams.fTimeRes);
+ cluster.finalize(*cell, fParams);
clustersOut.push_back(cluster);
}
+ //std::cout << "hits " << fiNbHits << std::endl;
return clustersOut;
}
-*/
+
} // namespace cbm::algo::tof
diff --git a/algo/detectors/tof/Clusterizer.h b/algo/detectors/tof/Clusterizer.h
index ac2776002e2d37e804d626d3a7965dd9694e278c..9fe5bdc807ed1c6292cf1d11612e412b87be7863 100644
--- a/algo/detectors/tof/Clusterizer.h
+++ b/algo/detectors/tof/Clusterizer.h
@@ -177,6 +177,8 @@ namespace cbm::algo::tof
std::vector<inputType> chanSortDigis(inputType& digisIn);
resultType buildClusters(std::vector<inputType>& input);
+ resultType buildClustersIter(std::vector<inputType>& input);
+
bool AddNextChan(std::vector<inputType>& input, int32_t iLastChan, TofCluster& cluster,
std::vector<TofCluster>& clustersOut);