34 const std::vector<double> &shape,
43 std::cout <<
"PulseFitWithShape::init: Error! Configuration samples in fit greater than total number of samples!" 120 double sq1 =
pshape[ims - 1];
122 double sq3 =
pshape[ims + 1];
124 double a2 = (sq3 + sq1) / 2.0 - sq2;
125 double a1 = sq2 - sq1 +
a2 * (1 - 2 * ims);
129 t_ims = -a1 / (2.0 *
a2);
140 for (
int is = 0; is < nsamplef; is++) {
141 if (
adc[is + nsampleo] > qm) {
142 qm =
adc[is + nsampleo];
148 if (qm > 5. *
noise) {
149 if (im >= nsamplef + nsampleo)
150 im = nsampleo + nsamplef - 1;
151 double q1 =
adc[im - 1];
153 double q3 =
adc[im + 1];
154 double y2 = (q1 + q3) / 2. - q2;
155 double y1 = q2 - q1 +
y2 * (1 - 2 * im);
156 double y0 = q2 -
y1 * (double)im -
y2 * (
double)(im * im);
158 xpar[0] = y0 +
y1 * tm +
y2 * tm * tm;
159 xpar[2] = (double)ims / 25. - tm;
162 double chi2old = 999999.;
164 int nsfit = nsamplef;
170 std::vector<double> resi(
fNsamples, 0.0);
172 while (std::fabs(chi2old -
chi2) > 0.1 && iloop < nloop) {
186 for (
int is = 0; is < nsfit; is++) {
187 const int iis = is + nsampleo;
189 double t1 = (double)iis + xpar[2];
190 double xbin =
t1 * 25.;
191 int ibin1 = (
int)xbin;
196 double amp1, amp11, amp12, der1, der11, der12;
200 int ibin2 = ibin1 + 1;
201 double xfrac = xbin - ibin1;
204 amp1 = (1. - xfrac) * amp11 + xfrac * amp12;
207 der1 = (1. - xfrac) * der11 + xfrac * der12;
216 for (
int js = 0; js < nsfit; js++) {
220 t1 = (double)jjs + xpar[2];
227 int ibin2 = ibin1 + 1;
228 double xfrac = xbin - ibin1;
231 double amp2 = (1. - xfrac) * amp11 + xfrac * amp12;
234 double der2 = (1. - xfrac) * der11 + xfrac * der12;
237 s1 = s1 + amp1 * cova[js *
fNsamples + is];
238 s2 = s2 + amp1 * amp2 * cova[js *
fNsamples + is];
239 ys1 = ys1 +
adc[iis] * amp2 * cova[js *
fNsamples + is];
240 sp1 = sp1 + der1 * cova[is *
fNsamples + js];
241 sp2 = sp2 + der1 * der2 * cova[js *
fNsamples + is];
242 ssp = ssp + amp1 * der2 * cova[js *
fNsamples + is];
243 ysp = ysp +
adc[iis] * der2 * cova[js *
fNsamples + is];
249 s2 = s2 + amp1 * amp1;
250 ys1 = ys1 + amp1 *
adc[iis];
252 sp2 = sp2 + der1 * der1;
253 ssp = ssp + der1 * amp1;
254 ysp = ysp + der1 *
adc[iis];
257 xpar[0] = (ysp * ssp - ys1 * sp2) / (ssp * ssp - s2 * sp2);
258 xpar[2] += (ysp / xpar[0] / sp2 - ssp / sp2);
260 for (
int is = 0; is < nsfit; is++) {
261 const int iis = is + nsampleo;
263 double t1 = (double)iis + xpar[2];
264 double xbin =
t1 * 25.;
265 int ibin1 = (
int)xbin;
273 int ibin2 = ibin1 + 1;
274 double xfrac = xbin - ibin1;
275 double amp11 = xpar[0] *
pshape[ibin1];
276 double amp12 = xpar[0] *
pshape[ibin2];
277 double amp1 = xpar[1] + (1. - xfrac) * amp11 + xfrac * amp12;
278 resi[iis] =
adc[iis] - amp1;
282 for (
int is = 0; is < nsfit; is++) {
283 const int iis = is + nsampleo;
286 for (
int js = 0; js < nsfit; js++) {
293 chi2 += resi[iis] * resi[iis];
virtual double doFit(double *, double *cova=nullptr)
virtual void init(int, int, int, int, int, const std::vector< double > &, double)
std::vector< double > dshape
std::vector< double > pshape
~PulseFitWithShape() override
uint16_t *__restrict__ uint16_t const *__restrict__ adc