diff --git a/algo/detectors/tof/Clusterizer.cxx b/algo/detectors/tof/Clusterizer.cxx index 6f504cc168aa498f2e79b39203f7a8b4a5882ec0..a2b48b0ef01f77fd2c223039f82e7b96ec552a6f 100644 --- a/algo/detectors/tof/Clusterizer.cxx +++ b/algo/detectors/tof/Clusterizer.cxx @@ -278,21 +278,28 @@ namespace cbm::algo::tof inputType& storDigi = input[chan]; if (0 == storDigi.size()) { return false; } + const double clusterTime = cluster.weightedTime / cluster.weightsSum; + const TofCell& cell = fParams.fChanPar[chan].cell; + + //Get starting position in buffer + auto i1 = (lastChanPos == nullptr) ? storDigi.begin() : (*lastChanPos)[chan]; + + //Find first digi in buffer that can possibly be part of the cluster. + //Relies on time-order to function properly. + i1 = std::lower_bound(i1, storDigi.end(), clusterTime - fParams.fdMaxTimeDist, + [this, cell](const auto& obj, double val) { + return obj.first->GetTime() + cell.sizeY * fParams.fPosYMaxScal / fParams.fSigVel < val; + }); + //Store last position in input channels to avoid unnecessary checks. //Relies on time-order to function properly. Can only be used in first-level AddNextChan() //calls, as time-order is not guaranteed in nested calls. - for (auto i1 = (lastChanPos == nullptr) ? storDigi.begin() : (*lastChanPos)[chan]; i1 < storDigi.end() - 1; i1++) { - - const CbmTofDigi* xDigiA = i1->first; - const double timeMax = - xDigiA->GetTime() + fParams.fChanPar[chan].cell.sizeY * fParams.fPosYMaxScal / fParams.fSigVel; + if (lastChanPos) { (*lastChanPos)[chan] = i1; } - if (timeMax < cluster.weightedTime / cluster.weightsSum - fParams.fdMaxTimeDist) { - if (lastChanPos) { (*lastChanPos)[chan]++; } - continue; - } - if (xDigiA->GetTime() > cluster.weightedTime / cluster.weightsSum + fParams.fdMaxTimeDist) { break; } + for (; i1 < storDigi.end() - 1; i1++) { + const CbmTofDigi* xDigiA = i1->first; + if (xDigiA->GetTime() > clusterTime + fParams.fdMaxTimeDist) { break; } auto i2 = i1; while (++i2 < storDigi.end()) { @@ -304,22 +311,26 @@ namespace cbm::algo::tof //D.Smith 8.9.23: Original continue statement, in general much slower (replaced by the conditions below). // The optimization relies on time-order. - //if (std::abs(time - cluster.weightedTime / cluster.weightsSum) >= fParams.fdMaxTimeDist) { continue; } + //if (std::abs(time - clusterTime) >= fParams.fdMaxTimeDist) { continue; } //Continue if digis are in the past of cluster time - if (cluster.weightedTime / cluster.weightsSum - time >= fParams.fdMaxTimeDist) { continue; } + if (time <= clusterTime - fParams.fdMaxTimeDist) { continue; } //Break if digis are in the future of cluster time - if (time - cluster.weightedTime / cluster.weightsSum >= fParams.fdMaxTimeDist) { break; } + if (time >= clusterTime + fParams.fdMaxTimeDist) { break; } const double timeDif = xDigiA->GetTime() - xDigiB->GetTime(); double posY = fParams.fSigVel * timeDif * 0.5; + + //Break if position is outside of the detector + if (std::abs(posY) > cell.sizeY * fParams.fPosYMaxScal) { break; } + 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) * fParams.fChanPar[chan].cell.sizeX; + 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.); @@ -337,10 +348,10 @@ namespace cbm::algo::tof break; } } - 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? + TofCell detcell = 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); + cluster.finalize(detcell, fParams); clustersOut.push_back(cluster); return true; } @@ -399,11 +410,8 @@ namespace cbm::algo::tof 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 timeDifN = (xDigiC->GetSide() == xDigiA->GetSide()) ? timeDifN = xDigiC->GetTime() - xDigiB->GetTime() + : timeDifN = xDigiC->GetTime() - xDigiA->GetTime(); double posYN = fParams.fSigVel * timeDifN * 0.5; if (xDigiC->GetSide() != 1.) { posYN *= -1.; }