11 #include <framework/logging/Logger.h>
12 #include <svd/reconstruction/SimpleClusterCandidate.h>
13 #include <vxd/geometry/GeoCache.h>
14 #include <svd/geometry/SensorInfo.h>
15 #include <framework/core/Environment.h>
17 #include <framework/datastore/StoreArray.h>
18 #include <svd/dataobjects/SVDRecoDigit.h>
30 SimpleClusterCandidate::SimpleClusterCandidate(
VxdID vxdID,
bool isUside,
int sizeHeadTail,
double cutSeed,
double cutAdjacent,
31 double cutCluster,
int timeAlgorithm)
34 , m_sizeHeadTail(sizeHeadTail)
36 , m_cutAdjacent(cutAdjacent)
37 , m_cutCluster(cutCluster)
38 , m_timeAlgorithm(timeAlgorithm)
39 , m_storeShaperDigitsName(
"SVDShaperDigits")
40 , m_storeRecoDigitsName(
"SVDRecoDigits")
45 , m_6SampleTimeError(0)
54 double cutCluster,
int timeAlgorithm, std::string storeShaperDigitsName, std::string storeRecoDigitsName)
57 , m_sizeHeadTail(sizeHeadTail)
59 , m_cutAdjacent(cutAdjacent)
60 , m_cutCluster(cutCluster)
61 , m_timeAlgorithm(timeAlgorithm)
62 , m_storeShaperDigitsName(storeShaperDigitsName)
63 , m_storeRecoDigitsName(storeRecoDigitsName)
68 , m_6SampleTimeError(0)
116 double pitch =
m_isUside ? info.getUPitch() : info.getVPitch();
120 int minStripCellID =
m_strips.at(0).cellID;
121 double minStripCharge =
m_strips.at(0).charge;
125 double weightSum = 0;
128 double stripPos =
m_isUside ? info.getUCellPosition(aStrip.cellID) : info.getVCellPosition(aStrip.cellID);
133 float tmp_sigmaSquared = aStrip.timeError * aStrip.timeError;
134 weightSum += tmp_sigmaSquared;
139 noise += aStrip.noise * aStrip.noise;
142 if ((
m_charge == 0) || (noise == 0)) {
157 if (clusterSize == 1) {
175 double centreCharge = (
m_charge - minStripCharge - maxStripCharge) / (clusterSize - 2);
176 if (centreCharge == 0) {
181 minStripCharge = (minStripCharge < centreCharge) ? minStripCharge : centreCharge;
182 maxStripCharge = (maxStripCharge < centreCharge) ? maxStripCharge : centreCharge;
183 double minPos =
m_isUside ? info.getUCellPosition(minStripCellID) : info.getVCellPosition(minStripCellID);
184 double maxPos =
m_isUside ? info.getUCellPosition(maxStripCellID) : info.getVCellPosition(maxStripCellID);
185 m_position = 0.5 * (minPos + maxPos + (maxStripCharge - minStripCharge) / centreCharge * pitch);
189 double landauHead = minStripCharge / centreCharge;
190 double landauTail = maxStripCharge / centreCharge;
192 0.5 * landauHead * landauHead +
193 0.5 * landauTail * landauTail);
215 B2WARNING(
"Something is wrong in the cluster creation, this cluster will not be created!");
236 B2FATAL(
"unrecognized timeAlgorithm, please check the input parameter and select a value among {0 (6-sample CoG), 1 (3-sample CoG), 2 (3-sample ELS)} to select the algorithm that you want to reconstruct the cluster time");
251 B2FATAL(
"unrecognized timeAlgorithm, please check the input parameter and select a value among {0 (6-sample CoG), 1 (3-sample CoG), 2 (3-sample ELS)} to select the algorithm that you want to reconstruct the cluster time");
261 auto begin = clustered3s.begin();
262 const auto end = clustered3s.end();
265 constexpr
auto stepSize = 16000. / 509;
266 auto retval = 0., norm = 0.;
267 for (
auto step = 0.; begin != end; ++begin, step += stepSize) {
268 norm +=
static_cast<double>(*begin);
269 retval +=
static_cast<double>(*begin) * step;
271 float rawtime = retval / norm;
281 const auto begin = clustered3s.begin();
284 constexpr
auto stepSize = 16000. / 509;
285 constexpr
auto tau = 55;
286 auto num = 2 * stepSize * std::exp(-4 * stepSize / tau) + std::exp(-2 * stepSize / tau) * stepSize / 2 * (*begin - std::exp(
287 -2 * stepSize / tau) * (*(begin + 2))) / (*begin + std::exp(-stepSize / tau) * (*(begin + 1)) / 2);
288 auto denom = 1 - std::exp(-4 * stepSize / tau) - (1 + std::exp(-2 * stepSize / tau) / 2) * (*begin - std::exp(
289 -2 * stepSize / tau) * (*(begin + 2))) / (*begin + std::exp(-stepSize / tau) * (*(begin + 1)) / 2);
290 float rawtime = - num / denom;
314 B2ERROR(
" you are asking fo the cluster samples for a cluster candidate with no strips, it make no sense to ask for the cluster time!");
330 if (!shaperdigit) B2ERROR(
"No shaperdigit for strip!?");
332 for (
int iSample = 0; iSample < static_cast<int>(APVsamples.size()); ++iSample)
333 returnSamples.at(iSample) += APVsamples.at(iSample);
336 return returnSamples;
346 if (clsSamples.size() < 3) B2ERROR(
"APV25 samples less than 3!?");
347 std::vector<float> Sum2bin(clsSamples.size() - 1, 0);
348 for (
int iBin = 0; iBin < static_cast<int>(clsSamples.size()) - 1; ++iBin)
349 Sum2bin.at(iBin) = clsSamples.at(iBin) + clsSamples.at(iBin + 1);
350 auto itSum = std::max_element(std::begin(Sum2bin), std::end(Sum2bin));
351 int ctrFrame = std::distance(std::begin(Sum2bin), itSum);
352 if (ctrFrame == 0) ctrFrame = 1;
353 std::vector<float> clustered3s = {clsSamples.at(ctrFrame - 1), clsSamples.at(ctrFrame), clsSamples.at(ctrFrame + 1)};
355 return std::make_pair(ctrFrame - 1, clustered3s);