EcalLaserAnalyzer/src/TSFit.cc

0001 /*
0002  *  \class TSFit
0003  *
0004  *  \author: Jean-Pierre Pansart - CEA/Saclay
0005  */
0006
0007 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TSFit.h>
0008
0009 #include <cstdio>
0010 #include <cmath>
0011 #include <cstring>
0012
0013 //ClassImp(TSFit)
0014
0015 //------------------------------------------------------------------------
0016 // TSFit provide fitting methods of pulses with third degree polynomial
0017 //
0018
0019 TSFit::TSFit(int size, int size_sh) {
0020   sdim = size;
0021   plshdim = size_sh;
0022   // sample_flag = new int[size];
0023   //   t = new double[sdim];
0024   //   z = new double[sdim];
0025   //   f = new double[sdim];
0026   //   acc = new double[sdim];
0027   //   adfmx = new double[sdim];
0028   //   adcp = new double[sdim];
0029   //   maskp3 = new double[sdim];
0030   //   corel = new double[sdim];
0031   //   nbcor = new double[sdim];
0032   //   int k;
0033   //   ff = new double *[sdim];
0034   //   for(k=0;k<sdim;k++)ff[k] = new double[4];
0035   //   der = new double *[sdim];
0036   //   for(k=0;k<sdim;k++)der[k] = new double[5];
0037 }
0038
0039 void TSFit::set_params(int n_samples,
0040                        int niter,
0041                        int n_presmpl,
0042                        int sample_min,
0043                        int sample_max,
0044                        double time_of_max,
0045                        double chi2_max,
0046                        int nsbm,
0047                        int nsam) {
0048   // parameters initialisation of the fit package
0049   // nsbm  n_samples_bef_max, nsam n_samples_aft_max
0050
0051   nbs = n_samples;
0052   nbr_iter_fit = niter;
0053   n_presamples = n_presmpl;
0054   iinf = sample_min;
0055   isup = sample_max;
0056   avtm = time_of_max;
0057   xki2_max = chi2_max;
0058   n_samples_bef_max = nsbm;
0059   n_samples_aft_max = nsam;
0060
0061   norme = 0.;
0062   alpha_th = 2.20;
0063   beta_th = 1.11;
0064
0065   int k;
0066
0067   for (k = 0; k <= nbs; k++) {
0068     sample_flag[k] = 0;
0069   }
0070
0071   for (k = sample_min; k <= sample_max; k++) {
0072     sample_flag[k] = 1;
0073   }
0074   /*
0075   int lim1 = ( iinf > n_presamples ) ? n_presamples : iinf;
0076   for(int k=lim1;k<=sample_max;k++){
0077     sample_flag[k] = 2;
0078   }
0079   */
0080   //  printf( "sample_fag : " );
0081   //   for(k=0;k<=nbs;k++){
0082   //     printf( "%1d ", sample_flag[k] );
0083   //   }
0084   //   printf( "\n" );
0085 }
0086
0087 void TSFit::init_errmat(double noise_initialvalue) {
0088   //  Input:  noise_initial value  noise (in adc channels) as read in the
0089   //  data base.
0090   /*------------------------------------------------------------------------*/
0091
0092   int i, j;
0093   //double one_over_noisesq;
0094
0095   //one_over_noisesq = 1. / ( noise_initialvalue * noise_initialvalue );
0096   for (i = 0; i < sdim; i++) {
0097     for (j = 0; j < sdim; j++) {
0098       errmat[i][j] = noise_initialvalue;
0099       //errmat[i][j] = 0.;
0100     }
0101     //errmat[i][i] = one_over_noisesq;
0102   }
0103 }
0104
0105 double TSFit::fpol3dg(int nmxul, double *parom, double *mask, double *adc) {
0106   // fit third degree polynomial
0107   // nmxul   array adc[] length
0108   // parom   return parameters (a0,a1,a2,a3,pos max,height)
0109   // fplo3dg uses only the diagonal terms of errmat[][]
0110   // errmat  inverse of the error matrix
0111
0112   int i, k, l;
0113   double h, t2, tm, delta, tmp;
0114   double xki2, dif, difmx, deglib;
0115   double bv[4], s;
0116
0117   deglib = (double)nmxul - 4.;
0118   for (i = 0; i < nmxul; i++) {
0119     t[i] = i;
0120     ff[i][0] = 1.;
0121     ff[i][1] = t[i];
0122     ff[i][2] = t[i] * t[i];
0123     ff[i][3] = ff[i][2] * t[i];
0124   }
0125   /*   computation of covariance matrix     */
0126   for (k = 0; k < 4; k++) {
0127     for (l = 0; l < 4; l++) {
0128       s = 0.;
0129       for (i = 0; i < nmxul; i++) {
0130         s = s + ff[i][k] * ff[i][l] * errmat[i][i] * mask[i];
0131       }
0132       cov[k][l] = s;
0133     }
0134     s = 0.;
0135     for (i = 0; i < nmxul; i++) {
0136       s = s + ff[i][k] * adc[i] * errmat[i][i] * mask[i];
0137     }
0138     bv[k] = s;
0139   }
0140   /*     parameters                          */
0141   inverms(4, cov, invcov);
0142   for (k = 0; k < 4; k++) {
0143     s = 0.;
0144     for (l = 0; l < 4; l++) {
0145       s = s + bv[l] * invcov[l][k];
0146     }
0147     parom[k] = s;
0148   }
0149
0150   if (parom[3] == 0.) {
0151     parom[4] = -1000.;
0152     parom[5] = -1000.;
0153     parom[6] = -1000.;
0154     return 1000000.;
0155   }
0156   /*    worst hit and ki2                    */
0157   xki2 = 0.;
0158   difmx = 0.;
0159   for (i = 0; i < nmxul; i++) {
0160     t2 = t[i] * t[i];
0161     h = parom[0] + parom[1] * t[i] + parom[2] * t2 + parom[3] * t2 * t[i];
0162     dif = (adc[i] - h) * mask[i];
0163     xki2 = xki2 + dif * dif * errmat[i][i];
0164     if (dif > difmx) {
0165       difmx = dif;
0166     }
0167   }
0168   if (deglib > 0.5)
0169     xki2 = xki2 / deglib;
0170   /*     amplitude and maximum position                    */
0171   delta = parom[2] * parom[2] - 3. * parom[3] * parom[1];
0172   if (delta > 0.) {
0173     delta = sqrt(delta);
0174     tm = -(delta + parom[2]) / (3. * parom[3]);
0175     tmp = (delta - parom[2]) / (3. * parom[3]);
0176   } else {
0177     parom[4] = -1000.;
0178     parom[5] = -1000.;
0179     parom[6] = -1000.;
0180     return xki2;
0181   }
0182   parom[4] = tm;
0183   parom[5] = parom[0] + parom[1] * tm + parom[2] * tm * tm + parom[3] * tm * tm * tm;
0184   parom[6] = tmp;
0185   return xki2;
0186 }
0187 double TSFit::inverms(int n, double g[matdim][matdim], double ginv[matdim][matdim]) {
0188   // inversion of a positive definite symetric matrix of size n
0189
0190   int i, j, k, jj;
0191   double r, s;
0192   double deter = 0;
0193
0194   /*   initialisation  */
0195
0196   if (n > matdim) {
0197     printf("ERROR : trying to use TSFit::inverms with size %d( max allowed %d\n", n, matdim);
0198     return -999.;
0199   }
0200
0201   int zero = 0;
0202   memset((char *)al, zero, 8 * n * n);
0203   memset((char *)be, zero, 8 * n * n);
0204   /*
0205   for(i=0;i<n;i++){
0206     for(j=0;j<n;j++){
0207       al[i][j] = 0.;
0208       be[i][j] = 0.;
0209     }
0210   }
0211   */
0212   /*  decomposition en vecteurs sur une base orthonormee  */
0213   al[0][0] = sqrt(g[0][0]);
0214   for (i = 1; i < n; i++) {
0215     al[i][0] = g[0][i] / al[0][0];
0216     for (j = 1; j <= i; j++) {
0217       s = 0.;
0218       for (k = 0; k <= j - 1; k++) {
0219         s = s + al[i][k] * al[j][k];
0220       }
0221       r = g[i][j] - s;
0222       if (j < i)
0223         al[i][j] = r / al[j][j];
0224       if (j == i)
0225         al[i][j] = sqrt(r);
0226     }
0227   }
0228   /*  inversion de la matrice al   */
0229   be[0][0] = 1. / al[0][0];
0230   for (i = 1; i < n; i++) {
0231     be[i][i] = 1. / al[i][i];
0232     for (j = 0; j < i; j++) {
0233       jj = i - j - 1;
0234       s = 0.;
0235       for (k = jj + 1; k <= i; k++) {
0236         s = s + be[i][k] * al[k][jj];
0237       }
0238       be[i][jj] = -s / al[jj][jj];
0239     }
0240   }
0241   /*   calcul de la matrice ginv   */
0242   for (i = 0; i < n; i++) {
0243     for (j = 0; j < n; j++) {
0244       s = 0.;
0245       for (k = 0; k < n; k++) {
0246         s = s + be[k][i] * be[k][j];
0247       }
0248       ginv[i][j] = s;
0249     }
0250   }
0251
0252   return deter;
0253 }
0254
0255 double TSFit::fit_third_degree_polynomial(double *bdc, double *ret_dat) {
0256   //  third degree polynomial fit of the pulse summit.
0257   //  samples are contained in array bdc and must be pedestal
0258   //  substracted.
0259   //  only samples having sample_flag >= 1 are used.
0260   //  the unit of time is one clock unit, that is to say 25 ns.
0261   //  output: ret_dat[0] = pulse height
0262   //          ret_dat[1]   position of maximum in the sample frame in clock units
0263   //          ret_dat[2]   adc value of the highest sample
0264   //          ret_dat[3]   number of the highest sample
0265   //          ret_dat[4]   lower sample number used for fitting
0266   //          ret_dat[5]   upper sample number used for fitting
0267   // errmat  inverse of the error matrix
0268
0269   int i;
0270   int nus;
0271   double xki2;
0272   double tm, tmp, amp;
0273
0274   ret_dat[0] = -999.;
0275   ret_dat[1] = -999.;
0276
0277   //    search the maximum
0278   double val_max = 0.;
0279   int imax = 0;
0280   for (i = 0; i < nbs; i++) {
0281     if (sample_flag[i] == 0)
0282       continue;
0283     if (bdc[i] > val_max) {
0284       val_max = bdc[i];
0285       imax = i;
0286     }
0287   }
0288
0289   if ((val_max * val_max) * errmat[imax][imax] < 16.)
0290     return -118;
0291
0292   //  if( imax != 9 )printf( "imax : %d !!!!!!!!!!!!!!!!!!!!!!!!!!!\n", imax );
0293
0294   if (norme == 0.)
0295     norme = val_max;
0296
0297   // look for samples above 1/3 of maximum before and 1/2 after
0298   double val2 = val_max / 2.;
0299   double val3 = val_max / 2.;
0300   int ilow = iinf;
0301   int ihig = 0;
0302
0303   for (i = iinf; i <= isup; i++) {
0304     if (sample_flag[i] >= 1) {
0305       if ((bdc[i] < val3) && (i < imax))
0306         ilow = i;
0307       if (bdc[i] > val2)
0308         ihig = i;
0309     }
0310   }
0311
0312   ilow++;
0313
0314   //ilow = imax - 1;
0315
0316   /*  le test suivant, apparemment idiot, est mis a cause des sequences 0. 2048. qui apparaissent dans certains mauvais evts     JPP 11/09/00 */
0317
0318   if (ihig == ilow)
0319     return -105;
0320   if (ilow == imax)
0321     ilow = ilow - 1;
0322   //if( ihig - ilow < 3 )ihig = ilow + 3;
0323   ihig = ilow + 3;
0324
0325   /*   printf("  third degree:   ilow %d ihig %d \n",ilow,ihig);  */
0326   nus = 0;
0327   int number_of_good_samples = 0;
0328   for (i = ilow; i <= ihig; i++) {
0329     maskp3[nus] = 0;
0330     adfmx[nus] = 0.;
0331     /*    printf(" adc %f sample_flag %d number_of good_samples %d \n",bdc[i],sample_flag[i],number_of_good_samples);  */
0332     if (sample_flag[i] >= 1) {
0333       adfmx[nus] = bdc[i];
0334       maskp3[nus] = 1.;
0335       number_of_good_samples++;
0336     }
0337     nus++;
0338   }
0339
0340   if (number_of_good_samples < 4) {
0341     return (-106);
0342   }
0343
0344   xki2 = fpol3dg(nus, &parfp3[0], &maskp3[0], &adfmx[0]);
0345
0346   /* printf( "fpol3dg-----------------------------------> %f %f %f %f %f\n",
0347       parfp3[0], parfp3[1], parfp3[2], parfp3[3], parfp3[4] );  */
0348
0349   tm = parfp3[4] + (float)ilow;
0350   amp = parfp3[5];
0351
0352   if (amp * amp * errmat[0][0] < 2.)
0353     return -101.;
0354   tmp = parfp3[6] + (float)ilow;
0355
0356   /*
0357     validation of fit quality.  Most of the time the fit is done with
0358     four samples, therefore there is no possible ki2 check. When more than
0359     4 samples are used the ki2 is often bad. So, in order to suppress some
0360     events with bad samples, a consistency check on the position of the
0361     maximum and minimum of the 3rd degree polynomial is used.
0362   */
0363
0364   if (xki2 > xki2_max) {
0365     return -102.;
0366   }
0367   if ((tm < (double)ilow) || (tm > (double)ihig)) {
0368     return -103.;
0369   }
0370
0371   if ((tmp > (double)ilow) && (tmp < (double)ihig - 1.)) {
0372     return -104.;
0373   }
0374
0375   ret_dat[0] = amp;
0376   ret_dat[1] = tm;
0377   ret_dat[2] = val_max;
0378   ret_dat[3] = (double)imax;
0379   ret_dat[4] = (double)ilow;
0380   ret_dat[5] = (double)ihig;
0381   ret_dat[6] = (double)tmp;
0382   int k;
0383   for (i = 0; i < 4; i++) {
0384     k = i + 7;
0385     ret_dat[k] = parfp3[i];
0386   }
0387
0388   return xki2;
0389 }