File indexing completed on 2023-03-17 10:41:25
0001
0002
0003
0004
0005
0006
0007 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithShape.h"
0008
0009 #include <iostream>
0010 #include "TMath.h"
0011 #include <cmath>
0012
0013
0014
0015
0016 PulseFitWithShape::PulseFitWithShape() {
0017 fNsamples = 0;
0018 fNsamplesShape = 0;
0019 fNum_samp_bef_max = 0;
0020 fNum_samp_after_max = 0;
0021 fNoise = 0.0;
0022 }
0023
0024
0025 PulseFitWithShape::~PulseFitWithShape() {}
0026
0027
0028
0029 void PulseFitWithShape::init(int n_samples,
0030 int samplb,
0031 int sampla,
0032 int niter,
0033 int n_samplesShape,
0034 const std::vector<double> &shape,
0035 double nois) {
0036 fNsamples = n_samples;
0037 fNsamplesShape = n_samplesShape;
0038 fNb_iter = niter;
0039 fNum_samp_bef_max = samplb;
0040 fNum_samp_after_max = sampla;
0041
0042 if (fNsamplesShape < fNum_samp_bef_max + fNum_samp_after_max + 1) {
0043 std::cout << "PulseFitWithShape::init: Error! Configuration samples in fit greater than total number of samples!"
0044 << std::endl;
0045 }
0046
0047 for (int i = 0; i < fNsamplesShape; i++) {
0048 pshape.push_back(shape[i]);
0049 dshape.push_back(0.0);
0050 }
0051
0052 fNoise = nois;
0053 return;
0054 }
0055
0056
0057
0058 double PulseFitWithShape::doFit(double *adc, double *cova) {
0059
0060
0061
0062
0063
0064 bool useCova = true;
0065 if (cova == nullptr)
0066 useCova = false;
0067
0068 double xpar[3];
0069 double chi2;
0070
0071 fAmp_fitted_max = 0.;
0072 fTim_fitted_max = 0.;
0073
0074
0075
0076 xpar[1] = 0.0;
0077
0078
0079
0080 double noise = fNoise;
0081
0082
0083 xpar[0] = 0.;
0084 xpar[2] = 0.;
0085
0086
0087
0088 int imax = 0;
0089 double amax = 0.0;
0090 for (int i = 0; i < fNsamples; i++) {
0091 if (adc[i] > amax) {
0092 amax = adc[i];
0093 imax = i;
0094 }
0095 }
0096
0097
0098
0099 double qms = 0.;
0100 int ims = 0;
0101
0102
0103
0104 for (int is = 0; is < fNsamplesShape; is++) {
0105 if (pshape[is] > qms) {
0106 qms = pshape[is];
0107 ims = is;
0108 }
0109
0110
0111
0112 if (is < fNsamplesShape - 2)
0113 dshape[is] = (pshape[is + 2] - pshape[is]) * 12.5;
0114 else
0115 dshape[is] = dshape[is - 1];
0116 }
0117
0118
0119
0120 double sq1 = pshape[ims - 1];
0121 double sq2 = pshape[ims];
0122 double sq3 = pshape[ims + 1];
0123
0124 double a2 = (sq3 + sq1) / 2.0 - sq2;
0125 double a1 = sq2 - sq1 + a2 * (1 - 2 * ims);
0126
0127 double t_ims = 0;
0128 if (a2 != 0)
0129 t_ims = -a1 / (2.0 * a2);
0130
0131
0132
0133
0134 double qm = 0.;
0135 int im = 0;
0136
0137 int nsamplef = fNum_samp_bef_max + fNum_samp_after_max + 1;
0138 int nsampleo = imax - fNum_samp_bef_max;
0139
0140 for (int is = 0; is < nsamplef; is++) {
0141 if (adc[is + nsampleo] > qm) {
0142 qm = adc[is + nsampleo];
0143 im = nsampleo + is;
0144 }
0145 }
0146
0147 double tm;
0148 if (qm > 5. * noise) {
0149 if (im >= nsamplef + nsampleo)
0150 im = nsampleo + nsamplef - 1;
0151 double q1 = adc[im - 1];
0152 double q2 = adc[im];
0153 double q3 = adc[im + 1];
0154 double y2 = (q1 + q3) / 2. - q2;
0155 double y1 = q2 - q1 + y2 * (1 - 2 * im);
0156 double y0 = q2 - y1 * (double)im - y2 * (double)(im * im);
0157 tm = -y1 / 2. / y2;
0158 xpar[0] = y0 + y1 * tm + y2 * tm * tm;
0159 xpar[2] = (double)ims / 25. - tm;
0160 }
0161
0162 double chi2old = 999999.;
0163 chi2 = 99999.;
0164 int nsfit = nsamplef;
0165 int iloop = 0;
0166 int nloop = fNb_iter;
0167 if (qm <= 5 * noise)
0168 nloop = 1;
0169
0170 std::vector<double> resi(fNsamples, 0.0);
0171
0172 while (std::fabs(chi2old - chi2) > 0.1 && iloop < nloop) {
0173 iloop++;
0174 chi2old = chi2;
0175
0176 double c = 0.;
0177 double y1 = 0.;
0178 double s1 = 0.;
0179 double s2 = 0.;
0180 double ys1 = 0.;
0181 double sp1 = 0.;
0182 double sp2 = 0.;
0183 double ssp = 0.;
0184 double ysp = 0.;
0185
0186 for (int is = 0; is < nsfit; is++) {
0187 const int iis = is + nsampleo;
0188
0189 double t1 = (double)iis + xpar[2];
0190 double xbin = t1 * 25.;
0191 int ibin1 = (int)xbin;
0192
0193 if (ibin1 < 0)
0194 ibin1 = 0;
0195
0196 double amp1, amp11, amp12, der1, der11, der12;
0197
0198 if (ibin1 <= fNsamplesShape - 2) {
0199
0200 int ibin2 = ibin1 + 1;
0201 double xfrac = xbin - ibin1;
0202 amp11 = pshape[ibin1];
0203 amp12 = pshape[ibin2];
0204 amp1 = (1. - xfrac) * amp11 + xfrac * amp12;
0205 der11 = dshape[ibin1];
0206 der12 = dshape[ibin2];
0207 der1 = (1. - xfrac) * der11 + xfrac * der12;
0208
0209 } else {
0210
0211 amp1 = pshape[fNsamplesShape - 1] + dshape[fNsamplesShape - 1] * (xbin - double(fNsamplesShape - 1)) / 25.;
0212 der1 = dshape[fNsamplesShape - 1];
0213 }
0214
0215 if (useCova) {
0216 for (int js = 0; js < nsfit; js++) {
0217 int jjs = js;
0218 jjs += nsampleo;
0219
0220 t1 = (double)jjs + xpar[2];
0221 xbin = t1 * 25.;
0222 ibin1 = (int)xbin;
0223 if (ibin1 < 0)
0224 ibin1 = 0;
0225 if (ibin1 > fNsamplesShape - 2)
0226 ibin1 = fNsamplesShape - 2;
0227 int ibin2 = ibin1 + 1;
0228 double xfrac = xbin - ibin1;
0229 amp11 = pshape[ibin1];
0230 amp12 = pshape[ibin2];
0231 double amp2 = (1. - xfrac) * amp11 + xfrac * amp12;
0232 der11 = dshape[ibin1];
0233 der12 = dshape[ibin2];
0234 double der2 = (1. - xfrac) * der11 + xfrac * der12;
0235 c = c + cova[js * fNsamples + is];
0236 y1 = y1 + adc[iis] * cova[js * fNsamples + is];
0237 s1 = s1 + amp1 * cova[js * fNsamples + is];
0238 s2 = s2 + amp1 * amp2 * cova[js * fNsamples + is];
0239 ys1 = ys1 + adc[iis] * amp2 * cova[js * fNsamples + is];
0240 sp1 = sp1 + der1 * cova[is * fNsamples + js];
0241 sp2 = sp2 + der1 * der2 * cova[js * fNsamples + is];
0242 ssp = ssp + amp1 * der2 * cova[js * fNsamples + is];
0243 ysp = ysp + adc[iis] * der2 * cova[js * fNsamples + is];
0244 }
0245 } else {
0246 c++;
0247 y1 = y1 + adc[iis];
0248 s1 = s1 + amp1;
0249 s2 = s2 + amp1 * amp1;
0250 ys1 = ys1 + amp1 * adc[iis];
0251 sp1 = sp1 + der1;
0252 sp2 = sp2 + der1 * der1;
0253 ssp = ssp + der1 * amp1;
0254 ysp = ysp + der1 * adc[iis];
0255 }
0256 }
0257 xpar[0] = (ysp * ssp - ys1 * sp2) / (ssp * ssp - s2 * sp2);
0258 xpar[2] += (ysp / xpar[0] / sp2 - ssp / sp2);
0259
0260 for (int is = 0; is < nsfit; is++) {
0261 const int iis = is + nsampleo;
0262
0263 double t1 = (double)iis + xpar[2];
0264 double xbin = t1 * 25.;
0265 int ibin1 = (int)xbin;
0266 if (ibin1 < 0)
0267 ibin1 = 0;
0268
0269 if (ibin1 < 0)
0270 ibin1 = 0;
0271 if (ibin1 > fNsamplesShape - 2)
0272 ibin1 = fNsamplesShape - 2;
0273 int ibin2 = ibin1 + 1;
0274 double xfrac = xbin - ibin1;
0275 double amp11 = xpar[0] * pshape[ibin1];
0276 double amp12 = xpar[0] * pshape[ibin2];
0277 double amp1 = xpar[1] + (1. - xfrac) * amp11 + xfrac * amp12;
0278 resi[iis] = adc[iis] - amp1;
0279 }
0280
0281 chi2 = 0.;
0282 for (int is = 0; is < nsfit; is++) {
0283 const int iis = is + nsampleo;
0284
0285 if (useCova) {
0286 for (int js = 0; js < nsfit; js++) {
0287 int jjs = js;
0288 jjs += nsampleo;
0289 chi2 += resi[iis] * resi[jjs] * cova[js * fNsamples + is];
0290 }
0291
0292 } else
0293 chi2 += resi[iis] * resi[iis];
0294 }
0295 }
0296
0297 fAmp_fitted_max = xpar[0];
0298 fTim_fitted_max = (double)t_ims / 25. - xpar[2];
0299
0300 return chi2;
0301 }