EcalLaserAnalyzer/src/PulseFitWithShape.cc

0001 /*
0002  *  \class PulseFitWithShape
0003  *
0004  *  \author: Julie Malcles - CEA/Saclay
0005  */
0006
0007 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithShape.h"
0008
0009 #include <iostream>
0010 #include "TMath.h"
0011 #include <cmath>
0012
0013 //ClassImp(PulseFitWithShape)
0014
0015 // Constructor...
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 // Destructor
0025 PulseFitWithShape::~PulseFitWithShape() {}
0026
0027 // Initialisation
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 // Compute the amplitude using as input the Crystaldata
0057
0058 double PulseFitWithShape::doFit(double *adc, double *cova) {
0059   // xpar = fit paramaters
0060   //     [0] = signal amplitude
0061   //     [1] = residual pedestal
0062   //     [2] = clock phase
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   // for now don't fit pedestal
0075
0076   xpar[1] = 0.0;
0077
0078   // Sample noise. If the cova matrix is defined, use it :
0079
0080   double noise = fNoise;
0081   //if(cova[0] > 0.) noise=1./sqrt(cova[0]);
0082
0083   xpar[0] = 0.;
0084   xpar[2] = 0.;
0085
0086   // first locate max:
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   // Shift pulse shape and calculate its derivative:
0098
0099   double qms = 0.;
0100   int ims = 0;
0101
0102   // 1) search for maximum
0103
0104   for (int is = 0; is < fNsamplesShape; is++) {
0105     if (pshape[is] > qms) {
0106       qms = pshape[is];
0107       ims = is;
0108     }
0109
0110     // 2) compute shape derivative :
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   // 3) compute pol2 max
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   // Starting point of the fit : qmax and tmax given by a
0132   // P2 fit on 3 max samples.
0133
0134   double qm = 0.;
0135   int im = 0;
0136
0137   int nsamplef = fNum_samp_bef_max + fNum_samp_after_max + 1;  // number of samples used in the fit
0138   int nsampleo = imax - fNum_samp_bef_max;                     // first sample number = sample max-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;  // Just one iteration for very low signal
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) {  // Interpolate shape values to get the right number :
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 {  // Or extraoplate if we are outside the array :
0210
0211         amp1 = pshape[fNsamplesShape - 1] + dshape[fNsamplesShape - 1] * (xbin - double(fNsamplesShape - 1)) / 25.;
0212         der1 = dshape[fNsamplesShape - 1];
0213       }
0214
0215       if (useCova) {  // Use covariance matrix:
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 {  // Don't use covariance matrix:
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 }