EcalLaserAnalyzer/src/TPNFit.cc

0001 /*
0002  *  \class TPNFit
0003  *
0004  *  \author: Patrice Verrecchia - CEA/Saclay
0005  */
0006
0007 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h>
0008
0009 #include <iostream>
0010 #include "TVectorD.h"
0011
0012 //ClassImp(TPNFit)
0013
0014 // Constructor...
0015 TPNFit::TPNFit() {
0016   fNsamples = 0;
0017   fNum_samp_bef_max = 0;
0018   fNum_samp_after_max = 0;
0019 }
0020
0021 // Destructor
0022 TPNFit::~TPNFit() {}
0023
0024 void TPNFit::init(int nsamples, int firstsample, int lastsample) {
0025   fNsamples = nsamples;
0026   fNum_samp_bef_max = firstsample;
0027   fNum_samp_after_max = lastsample;
0028   //printf("nsamples=%d firstsample=%d lastsample=%d\n",nsamples,firstsample,lastsample);
0029
0030   if (fNsamples > NMAXSAMP2)
0031     printf("number of PN samples exceed maximum\n");
0032
0033   for (int k = 0; k < NMAXSAMP2; k++)
0034     t[k] = (double)k;
0035
0036   return;
0037 }
0038
0039 double TPNFit::doFit(int maxsample, double *adc) {
0040   double chi2 = 0.;
0041   ampl = 0.;
0042   timeatmax = 0.;
0043
0044   //printf("maxsample=%d\n",maxsample);
0045   firstsample = maxsample - fNum_samp_bef_max;
0046   lastsample = maxsample + fNum_samp_after_max;
0047
0048   if (firstsample <= 0)
0049     return 101;
0050   if (lastsample >= fNsamples)
0051     lastsample = fNsamples - 1;
0052   if (lastsample - firstsample < 1)
0053     return 102;
0054   int nval = lastsample - firstsample + 1;
0055   //printf("firstsample=%d lastsample=%d nval=%d\n",
0056   //                        firstsample,lastsample,nval);
0057   int testneg = 0;
0058   for (int kn = firstsample; kn <= lastsample; kn++) {
0059     //printf("adc[%d]=%f\n",kn,adc[kn]);
0060     if (adc[kn] < 0.)
0061       testneg = 1;
0062   }
0063   if (testneg == 1)
0064     return 103;
0065
0066   for (int i = firstsample; i <= lastsample; i++) {
0067     val[i - firstsample] = adc[i];
0068     fv1[i - firstsample] = 1.;
0069     fv2[i - firstsample] = (double)(i);
0070     fv3[i - firstsample] = ((double)(i)) * ((double)(i));
0071   }
0072
0073   TVectorD y(nval, val);
0074   //y.Print();
0075   TVectorD f1(nval, fv1);
0076   //f1.Print();
0077   TVectorD f2(nval, fv2);
0078   //f2.Print();
0079   TVectorD f3(nval, fv3);
0080   //f3.Print();
0081
0082   double bj[3];
0083   bj[0] = f1 * y;
0084   bj[1] = f2 * y;
0085   bj[2] = f3 * y;
0086   TVectorD b(3, bj);
0087   //b.Print();
0088
0089   double aij[9];
0090   aij[0] = f1 * f1;
0091   aij[1] = f1 * f2;
0092   aij[2] = f1 * f3;
0093   aij[3] = f2 * f1;
0094   aij[4] = f2 * f2;
0095   aij[5] = f2 * f3;
0096   aij[6] = f3 * f1;
0097   aij[7] = f3 * f2;
0098   aij[8] = f3 * f3;
0099   TMatrixD a(3, 3, aij);
0100   //a.Print();
0101
0102   double det1;
0103   a.InvertFast(&det1);
0104   //a.Print();
0105
0106   TVectorD c = a * b;
0107   //c.Print();
0108
0109   double *par = c.GetMatrixArray();
0110   //printf("par0=%f par1=%f par2=%f\n",par[0],par[1],par[2]);
0111   if (par[2] != 0.) {
0112     timeatmax = -par[1] / (2. * par[2]);
0113     ampl = par[0] - par[2] * (timeatmax * timeatmax);
0114   }
0115   //printf("amp=%f time=%f\n",amp_max,timeatmax);
0116
0117   if (ampl <= 0.) {
0118     ampl = adc[maxsample];
0119     return 1.;
0120   }
0121
0122   if ((int)timeatmax > lastsample)
0123     return 103;
0124   if ((int)timeatmax < firstsample)
0125     return 103;
0126
0127   return chi2;
0128 }