59 if (fhc_vec.empty())
return;
61 const float chiCut = 1
e-3;
63 float chiSqr = std::numeric_limits<float>::max(), dchiSqr = std::numeric_limits<float>::max();
66 int startTime = fhc_vec.front().startTick;
67 int endTime = fhc_vec.back().stopTick;
69 int roiSize = endTime - startTime;
74 std::vector<float>
y(roiSize);
75 std::vector<float> p(3 * fhc_vec.size());
76 std::vector<float> plimmin(3 * fhc_vec.size());
77 std::vector<float> plimmax(3 * fhc_vec.size());
78 std::vector<float> perr(3 * fhc_vec.size());
83 for (
size_t ih = 0; ih < fhc_vec.size(); ih++) {
84 float const peakMean =
85 fhc_vec[ih].hitCenter - 0.5 - (float)startTime;
86 float const peakWidth = fhc_vec[ih].hitSigma;
87 float const amplitude = fhc_vec[ih].hitHeight;
88 float const meanLowLim = fmax(peakMean -
fPeakRange * peakWidth, 0.);
89 float const meanHiLim = fmin(peakMean +
fPeakRange * peakWidth, (
float)roiSize);
90 p[0 + nParams] = amplitude;
91 p[1 + nParams] = peakMean;
92 p[2 + nParams] = peakWidth;
94 plimmin[0 + nParams] = amplitude * 0.1;
95 plimmin[1 + nParams] = meanLowLim;
96 plimmin[2 + nParams] = std::max(
fMinWidth, 0.1 * peakWidth);
98 plimmax[0 + nParams] = amplitude *
fAmpRange;
99 plimmax[1 + nParams] = meanHiLim;
106 for (
size_t idx = 0; idx < size_t(roiSize); idx++) {
107 y[idx] = signal[startTime + idx];
113 fitResult =
fMarqFitAlg->gshf::MarqFitAlg::mrqdtfit(
114 lambda, &p[0], &plimmin[0], &plimmax[0], &
y[0], nParams, roiSize, chiSqr, dchiSqr);
116 if (fitResult || (trial > 100))
break;
117 }
while (fabs(dchiSqr) >= chiCut);
121 fMarqFitAlg->gshf::MarqFitAlg::cal_perr(&p[0], &
y[0], nParams, roiSize, &perr[0]);
123 NDF = roiSize - nParams;
124 chi2PerNDF = chiSqr / NDF;
126 for (
size_t i = 0; i < fhc_vec.size(); i++) {
137 PeakFitParams_t mhpp;
138 mhpp.peakAmplitude = p[parIdx + 0];
139 mhpp.peakAmplitudeError = perr[parIdx + 0];
141 p[parIdx + 1] + 0.5 + float(startTime);
142 mhpp.peakCenterError = perr[parIdx + 1];
143 mhpp.peakSigma =
std::abs(p[parIdx + 2]);
144 mhpp.peakSigmaError = perr[parIdx + 2];
146 mhpp_vec.emplace_back(mhpp);
constexpr auto abs(T v)
Returns the absolute value of the argument.