EcalLaserAnalyzer/src/PulseFitWithFunction.cc

0001 /*
0002  *  \class PulseFitWithFunction
0003  *
0004  *  \author: Patrick Jarry - CEA/Saclay
0005  */
0006
0007 // File PulseFitWithFunction.cxx
0008 // ===========================================================
0009 // ==                                                       ==
0010 // ==     Class for a function fit method                   ==
0011 // ==                                                       ==
0012 // ==  Date:   July 17th 2003                               ==
0013 // ==  Author: Patrick Jarry                                ==
0014 // ==                                                       ==
0015 // ==                                                       ==
0016 // ===========================================================
0017
0018 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithFunction.h"
0019
0020 #include <iostream>
0021 #include "TMath.h"
0022
0023 //ClassImp(PulseFitWithFunction)
0024
0025 // Constructor...
0026 PulseFitWithFunction::PulseFitWithFunction() {
0027   fNsamples = 0;
0028   fNum_samp_bef_max = 0;
0029   fNum_samp_after_max = 0;
0030 }
0031
0032 // Destructor
0033 PulseFitWithFunction::~PulseFitWithFunction() {}
0034
0035 // Initialisation
0036
0037 void PulseFitWithFunction::init(int n_samples, int samplb, int sampla, int niter, double alfa, double beta) {
0038   //printf("\n");
0039   //printf(" =========================================================\n");
0040   //printf(" ==     Initialising the function fit method            ==\n");
0041   //printf(" ==          PulseFitWithFunction::init                 ==\n");
0042   //printf(" ==                                                     ==\n");
0043
0044   fNsamples = n_samples;
0045   fAlpha_laser = alfa;
0046   fBeta_laser = beta;
0047   fAlpha_beam = 0.98;
0048   fBeta_beam = 2.04;
0049   fNb_iter = niter;
0050   fNum_samp_bef_max = samplb;
0051   fNum_samp_after_max = sampla;
0052   //printf(" ==  # samples used = %3d                               ==\n",fNsamples);
0053   //printf(" ==  #sample before max= %1d and #sample after maximum= %1d ==\n",fNum_samp_bef_max,fNum_samp_after_max);
0054   //printf(" ==           alpha= %5.4f       beta= %5.4f          ==\n",fAlpha_laser,fBeta_laser);
0055   //printf(" =========================================================\n\n");
0056   return;
0057 }
0058
0059 // Compute the amplitude using as input the Crystaldata
0060 double PulseFitWithFunction::doFit(double *adc) {
0061   double parout[4];  // amp_max ;
0062                      //double amp_parab , tim_parab ;
0063   double chi2;
0064   //int imax ;
0065   //
0066   // first  one has to get starting point first with parabolic fun// //
0067   Fit_parab(&adc[0], 3, fNsamples, parout);
0068   amp_parab = parout[0];
0069   tim_parab = parout[1];
0070   imax = (int)parout[2];
0071   amp_max = parout[3];
0072   //printf("amp_parab= %f tim_parab=%f amp_max=%f imax=%d\n",amp_parab,tim_parab,amp_max,imax);
0073   fNumber_samp_max = imax;
0074   //
0075
0076   if (amp_parab < 1.) {
0077     tim_parab = (double)imax;
0078     amp_parab = amp_max;
0079   }
0080   //
0081   fValue_tim_max = tim_parab;
0082   fFunc_max = amp_parab;
0083   fTim_max = tim_parab;
0084
0085   //  here to fit maximum amplitude and time of arrival ...
0086
0087   chi2 = Fit_electronic(0, &adc[0], 8.);
0088   // adc is an array to be filled with samples
0089   // 0 is for Laser (1 for electron)
0090   // 8 is for sigma of pedestals
0091   // which (can be computed)
0092
0093   //  double amplitude = fAmp_fitted_max ; // amplitude fitted
0094   // double time = fTim_fitted_max ; // time fitted
0095
0096   return chi2;  // return amplitude fitted
0097 }
0098
0099 //-----------------------------------------------------------------------
0100 //----------------------------------------------------------------------
0101
0102 /*************************************************/
0103 double PulseFitWithFunction::Fit_electronic(int data, double *adc_to_fit, double sigmas_sample) {
0104   // fit electronic function from simulation
0105   // parameters fAlpha and fBeta are fixed and fit is providing
0106   // the maximum amplitude ( fAmp_fitted_max ) and the time of
0107   // the maximum amplitude ( fTim_fitted_max)
0108   // initialization of parameters
0109   double chi2 = 0;
0110   double d_alpha, d_beta;
0111   // first initialize parameters fAlpha and fBeta ( depending of beam or laser)
0112   fAlpha = fAlpha_laser;
0113   fBeta = fBeta_laser;
0114   if (data == 1) {
0115     fAlpha = fAlpha_beam;
0116     fBeta = fBeta_beam;
0117   }
0118   //
0119   fAmp_fitted_max = 0.;
0120   fTim_fitted_max = 0.;
0121   double un_sur_sigma = 1.;
0122   double variation_func_max = 0.;
0123   double variation_tim_max = 0.;
0124   //
0125   if (fValue_tim_max > 20. || fValue_tim_max < 3.) {
0126     fValue_tim_max = fNumber_samp_max;
0127   }
0128   int num_fit_min = (int)(fValue_tim_max - fNum_samp_bef_max);
0129   int num_fit_max = (int)(fValue_tim_max + fNum_samp_after_max);
0130   //
0131
0132   if (sigmas_sample > 0.)
0133     un_sur_sigma = 1. / sigmas_sample;
0134
0135   double func, delta;
0136   //          Loop on iterations
0137   for (int iter = 0; iter < fNb_iter; iter++) {
0138     //          initialization inside iteration loop !
0139     chi2 = 0.;
0140     double d11 = 0.;
0141     double d12 = 0.;
0142     double d22 = 0.;
0143     double z1 = 0.;
0144     double z2 = 0.;
0145     fFunc_max += variation_func_max;
0146     fTim_max += variation_tim_max;
0147     int nsamp_used = 0;
0148     //
0149
0150     // Then we loop on samples to be fitted
0151
0152     for (int i = num_fit_min; i < num_fit_max + 1; i++) {
0153       // calculate function to be fitted
0154
0155       func = Electronic_shape((double)i);
0156
0157       // then calculate derivatives of function to be fitted
0158       double dt = (double)i - fTim_max;
0159       double alpha_beta = fAlpha * fBeta;
0160
0161       if (dt > -alpha_beta) {
0162         double dt_sur_beta = dt / fBeta;
0163
0164         double variable = (double)1. + dt / alpha_beta;
0165         double expo = TMath::Exp(-dt_sur_beta);
0166
0167         double puissance = TMath::Power(variable, fAlpha);
0168         d_alpha = un_sur_sigma * puissance * expo;
0169         d_beta = fFunc_max * d_alpha * dt_sur_beta / (alpha_beta * variable);
0170       } else {
0171         continue;
0172       }
0173
0174       nsamp_used++;  // number of samples used inside the fit
0175       // compute matrix elements D (symetric --> d12 = d21 )
0176       d11 += d_alpha * d_alpha;
0177       d12 += d_alpha * d_beta;
0178       d22 += d_beta * d_beta;
0179       // compute delta
0180       delta = (adc_to_fit[i] - func) * un_sur_sigma;
0181       // compute vector elements Z
0182       z1 += delta * d_alpha;
0183       z2 += delta * d_beta;
0184       chi2 += delta * delta;
0185     }  //             end of loop on samples
0186     double denom = d11 * d22 - d12 * d12;
0187     if (denom == 0.) {
0188       //printf( "attention denom = 0 signal pas fitte \n") ;
0189       return 101;
0190     }
0191     if (nsamp_used < 3) {
0192       //printf( "Attention nsamp = %d ---> no function fit provided \n",nsamp_used) ;
0193       return 102;
0194     }
0195     // compute variations of parameters fAmp_max and fTim_max
0196     variation_func_max = (z1 * d22 - z2 * d12) / denom;
0197     variation_tim_max = (-z1 * d12 + z2 * d11) / denom;
0198     chi2 = chi2 / ((double)nsamp_used - 2.);
0199   }  //      end of loop on iterations
0200   // results of the fit are calculated
0201   fAmp_fitted_max = fFunc_max + variation_func_max;
0202   fTim_fitted_max = fTim_max + variation_tim_max;
0203   //
0204   return chi2;
0205 }
0206
0207 //-----------------------------------------------------------------------
0208 //----------------------------------------------------------------------
0209
0210 double PulseFitWithFunction::Electronic_shape(double tim) {
0211   // electronic function (from simulation) to fit ECAL pulse shape
0212   double func_electronic, dtsbeta, variable, puiss;
0213   double albet = fAlpha * fBeta;
0214   if (albet <= 0)
0215     return ((Double_t)0.);
0216   double dt = tim - fTim_max;
0217   if (dt > -albet) {
0218     dtsbeta = dt / fBeta;
0219     variable = 1. + dt / albet;
0220     puiss = TMath::Power(variable, fAlpha);
0221     func_electronic = fFunc_max * puiss * TMath::Exp(-dtsbeta);
0222   } else
0223     func_electronic = 0.;
0224   //
0225   return func_electronic;
0226 }
0227 void PulseFitWithFunction::Fit_parab(Double_t *ampl, Int_t nmin, Int_t nmax, Double_t *parout) {
0228   /* Now we calculate the parabolic adjustement in order to get        */
0229   /*    maximum and time max                                           */
0230
0231   double denom, dt, amp1, amp2, amp3;
0232   double ampmax = 0.;
0233   int imax = 0;
0234   int k;
0235   /*
0236                                                                    */
0237   for (k = nmin; k < nmax; k++) {
0238     //printf("ampl[%d]=%f\n",k,ampl[k]);
0239     if (ampl[k] > ampmax) {
0240       ampmax = ampl[k];
0241       imax = k;
0242     }
0243   }
0244   amp1 = ampl[imax - 1];
0245   amp2 = ampl[imax];
0246   amp3 = ampl[imax + 1];
0247   denom = 2. * amp2 - amp1 - amp3;
0248   /*                                         */
0249   if (denom > 0.) {
0250     dt = 0.5 * (amp3 - amp1) / denom;
0251   } else {
0252     //printf("denom =%f\n",denom)  ;
0253     dt = 0.5;
0254   }
0255   /*                                             */
0256   /* ampmax correspond au maximum d'amplitude parabolique et dt        */
0257   /* decalage en temps par rapport au sample maximum soit k + dt       */
0258
0259   parout[0] = amp2 + (amp3 - amp1) * dt * 0.25;
0260   parout[1] = (double)imax + dt;
0261   parout[2] = (double)imax;
0262   parout[3] = ampmax;
0263   return;
0264 }