EcalLaserAnalyzer/src/TMatacq.cc

0001 /*
0002  *  \class TMatacq
0003  *
0004  *  \author: Patrice Verrecchia - CEA/Saclay
0005  */
0006 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TMatacq.h"
0007
0008 #include <iostream>
0009 #include <cmath>
0010 #include "TVectorD.h"
0011
0012 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TMarkov.h"
0013
0014 using namespace std;
0015 //ClassImp(TMatacq)
0016
0017 void TMatacq::init() {
0018   for (int k = 0; k < NMAXSAMP; k++)
0019     bong[k] = 0.;
0020
0021   for (int k = 0; k <= 100; k++)
0022     bing[k] = 0;
0023
0024   for (int k = 0; k < NSPARAB; k++)
0025     t[k] = (double)k;
0026
0027   return;
0028 }
0029
0030 // Constructor...
0031 TMatacq::TMatacq(
0032     int Ntot, int Nsamp1, int Nsamp2, int cut, int Nbef, int Naft, int niv1, int niv2, int niv3, int nevl, int NSlide) {
0033   fNsamples = Ntot;
0034   presample = Nsamp1;
0035   endsample = Nsamp2;
0036   nsigcut = (double)cut;
0037   fNum_samp_bef_max = Nbef;
0038   fNum_samp_aft_max = Naft;
0039   level1 = ((double)niv1) / 100.;
0040   level2 = ((double)niv2) / 100.;
0041   level3 = ((double)niv3) / 100.;
0042   nevlasers = nevl;
0043   slidingmean = 0.0;
0044   nslide = NSlide;
0045   for (int k = 0; k < nevlasers; k++)
0046     status[k] = 0;
0047   for (int k = 0; k < nevlasers; k++)
0048     status[k + nevlasers] = 0;
0049
0050   nevmtq0 = 0;
0051   nevmtq1 = 0;
0052 }
0053
0054 // Destructor
0055 TMatacq::~TMatacq() {}
0056
0057 int TMatacq::rawPulseAnalysis(Int_t Nsamp, Double_t *adc)  // GHM
0058 {
0059   using namespace std;
0060
0061   //  std::cout << "Entering rawPulseAnalysis" << std::endl;
0062
0063   int k, ithr;
0064   double dsum = 0., dsum2 = 0.;
0065
0066   //  std::cout << "calling init" << std::endl;
0067   init();
0068   //  std::cout << ".......done" << std::endl;
0069
0070   if (Nsamp != fNsamples) {
0071     printf("found different number of samples fNsamples=%d Nsamp=%d\n", fNsamples, Nsamp);
0072     return 100;
0073   }
0074
0075   for (k = 0; k < presample; k++) {
0076     dsum += adc[k];
0077     dsum2 += adc[k] * adc[k];
0078   }
0079   bl = dsum / ((double)presample);
0080   double ss = (dsum2 / ((double)presample) - bl * bl);
0081   if (ss < 0.)
0082     ss = 0.;
0083   sigbl = sqrt(ss);
0084   for (ithr = 0, k = presample; k < endsample; k++) {
0085     if (adc[k] > (bl + nsigcut * sigbl) && ithr == 0) {
0086       ithr = 1;
0087       firstsample = k;
0088     }
0089   }
0090
0091   if (ithr == 0)
0092     return 101;
0093
0094   for (ithr = 0, k = firstsample; k < Nsamp; k++) {
0095     if (adc[k] < (bl + nsigcut * sigbl) && ithr == 0) {
0096       ithr = 1;
0097       lastsample = k;
0098     }
0099   }
0100   if (ithr == 0)
0101     lastsample = Nsamp;
0102
0103   if (lastsample > firstsample + NMAXSAMP)
0104     lastsample = firstsample + NMAXSAMP;
0105
0106   val_max = 0.;
0107   samplemax = 0;
0108   for (Int_t is = firstsample; is < lastsample; is++) {
0109     bong[is - firstsample] = adc[is] - bl;
0110     if (bong[is - firstsample] > val_max) {
0111       val_max = bong[is - firstsample];
0112       samplemax = is;
0113     }
0114   }
0115   if (samplemax == 0)
0116     return 103;
0117   if (samplemax > lastsample)
0118     return 104;
0119   if (samplemax < firstsample)
0120     return 105;
0121
0122   int endslide = samplemax - nslide;
0123   int beginslide = nslide;
0124   int islidingmean = 0;
0125   slidingmean = 0.0;
0126
0127   for (int i = beginslide; i < endslide; i++) {
0128     slidingmean += adc[i];
0129     islidingmean += 1;
0130   }
0131   if (islidingmean != 0)
0132     slidingmean /= double(islidingmean);
0133
0134   return 0;
0135 }
0136 int TMatacq::findPeak() {
0137   int k;
0138   int nbinf = 0;
0139   int jfind = 0;
0140   int nbsup = 0;
0141   double thres = val_max * level1;
0142
0143   for (k = 0, jfind = 0; k < NMAXSAMP; k++) {
0144     if (jfind == 0) {
0145       if (bong[k] > thres) {
0146         nbinf = k;
0147         jfind = 1;
0148       }
0149     }
0150   }
0151   if (jfind == 0)
0152     nbinf = 0;
0153
0154   for (k = NMAXSAMP, jfind = 0; k > nbinf; k--) {
0155     if (jfind == 0) {
0156       if (bong[k] > thres) {
0157         nbsup = k;
0158         jfind = 1;
0159       }
0160     }
0161   }
0162   if (nbsup == 0)
0163     nbsup = nbinf;
0164
0165   pkval = 0.;
0166   sigpkval = 0.5;
0167   if (nbsup == nbinf) {
0168     return 301;
0169   } else {
0170     if (nbinf > 4)
0171       nbinf -= 3;
0172     else
0173       nbinf = 1;
0174     if (nbsup < NMAXSAMP - 4)
0175       nbsup += 3;
0176     else
0177       nbsup = NMAXSAMP;
0178
0179     for (k = 0; k < nbinf; k++)
0180       bong[k] = 0.;
0181     for (k = nbsup; k < NMAXSAMP; k++)
0182       bong[k] = 0.;
0183
0184     for (k = 0, jfind = 0; k < NMAXSAMP; k++) {
0185       if (bong[k] > 0.)
0186         jfind++;
0187     }
0188     if (jfind == 0) {
0189       return 302;
0190     } else if (jfind == 1) {
0191       return 303;
0192     } else {
0193       for (k = 0; k < NMAXSAMP; k++) {
0194         if (k < 100)
0195           bing[k + 1] = (int)bong[k];
0196       }
0197
0198       TMarkov *peak = new TMarkov();
0199
0200       peak->peakFinder(&bing[0]);
0201       pkval = peak->getPeakValue(0);
0202       sigpkval = peak->getPeakValue(1);
0203
0204       delete peak;
0205
0206       pkval += (firstsample - 1);
0207     }
0208   }
0209
0210   return 0;
0211 }
0212
0213 int TMatacq::doFit() {
0214   int testneg = 0;
0215   ampl = 0.;
0216   timeatmax = 0.;
0217   int heresamplemax = samplemax - firstsample;
0218
0219   int beginfit = heresamplemax - fNum_samp_bef_max;
0220   int endfit = heresamplemax + fNum_samp_aft_max;
0221
0222   int nval = endfit - beginfit + 1;
0223   if (nval != NSPARAB)
0224     return 201;
0225   for (int kn = beginfit; kn <= endfit; kn++) {
0226     if (bong[kn] <= 0)
0227       testneg = 1;
0228   }
0229   if (testneg == 1)
0230     return 202;
0231
0232   for (int i = 0; i < nval; i++) {
0233     val[i] = bong[beginfit + i];
0234     fv1[i] = 1.;
0235     fv2[i] = (double)(i);
0236     fv3[i] = ((double)(i)) * ((double)(i));
0237   }
0238
0239   TVectorD y(nval, val);
0240   TVectorD f1(nval, fv1);
0241   TVectorD f2(nval, fv2);
0242   TVectorD f3(nval, fv3);
0243
0244   double bj[3];
0245   bj[0] = f1 * y;
0246   bj[1] = f2 * y;
0247   bj[2] = f3 * y;
0248   TVectorD b(3, bj);
0249
0250   double aij[9];
0251   aij[0] = f1 * f1;
0252   aij[1] = f1 * f2;
0253   aij[2] = f1 * f3;
0254   aij[3] = f2 * f1;
0255   aij[4] = f2 * f2;
0256   aij[5] = f2 * f3;
0257   aij[6] = f3 * f1;
0258   aij[7] = f3 * f2;
0259   aij[8] = f3 * f3;
0260   TMatrixD a(3, 3, aij);
0261
0262   double det1;
0263   a.InvertFast(&det1);
0264
0265   TVectorD c = a * b;
0266
0267   double *par = c.GetMatrixArray();
0268   if (par[2] != 0.) {
0269     timeatmax = -par[1] / (2. * par[2]);
0270     ampl = par[0] - par[2] * (timeatmax * timeatmax);
0271   }
0272
0273   if (ampl <= 0.) {
0274     ampl = bong[heresamplemax];
0275     return 203;
0276   }
0277
0278   if ((int)timeatmax > NSPARAB)
0279     return 204;
0280   if ((int)timeatmax < 0)
0281     return 205;
0282
0283   timeatmax += (beginfit + firstsample);
0284
0285   int tplus[3], tminus[3];
0286   double xampl[3];
0287   xampl[0] = 0.2 * ampl;
0288   xampl[1] = 0.5 * ampl;
0289   xampl[2] = 0.8 * ampl;
0290
0291   int hitplus[3];
0292   int hitminus[3];
0293   double width[3];
0294
0295   for (int i = 0; i < 3; i++) {
0296     hitplus[i] = 0;
0297     hitminus[i] = 0;
0298     width[i] = 0.0;
0299     tplus[i] = firstsample + lastsample;
0300     tminus[i] = 0;
0301   }
0302
0303   // calculate first estimate of half width
0304   int im = heresamplemax;
0305   int iplusbound = firstsample + lastsample - im;
0306
0307   for (int j = 0; j < 3; j++) {
0308     for (int i = 1; i < im; i++) {
0309       if (bong[im - i] < xampl[j] && bong[im - i + 1] > xampl[j]) {
0310         tminus[j] = im - i;
0311         hitminus[j]++;
0312         i = im;
0313       }
0314     }
0315
0316     for (int i = 0; i < iplusbound; i++) {
0317       if (bong[im + i] > xampl[j] && bong[im + i + 1] < xampl[j]) {
0318         tplus[j] = im + i;
0319         hitplus[j]++;
0320         i = iplusbound;
0321       }
0322     }
0323
0324     // interpolate to get a better estimate
0325
0326     double slopeplus = double(bong[tplus[j] + 1] - bong[tplus[j]]);
0327     double slopeminus = double(bong[tminus[j] + 1] - bong[tminus[j]]);
0328
0329     double timeminus = double(tminus[j]) + 0.5;
0330     if (slopeminus != 0)
0331       timeminus = tminus[j] + (xampl[j] - double(bong[tminus[j]])) / slopeminus;
0332
0333     double timeplus = double(tplus[j]) + 0.5;
0334     if (slopeplus != 0)
0335       timeplus = tplus[j] + (xampl[j] - double(bong[tplus[j]])) / slopeplus;
0336
0337     width[j] = timeplus - timeminus;
0338   }
0339
0340   width20 = width[0];
0341   width50 = width[1];
0342   width80 = width[2];
0343
0344   return 0;
0345 }
0346
0347 int TMatacq::compute_trise() {
0348   int error;
0349   trise = 0.;
0350
0351   double t20 = interpolate(ampl * level2);
0352   if (t20 < 0.) {
0353     error = (int)-t20;
0354     return error;
0355   }
0356   double t80 = interpolate(ampl * level3);
0357   if (t80 < 0.) {
0358     error = (int)-t80;
0359     return error;
0360   }
0361   trise = t80 - t20;
0362
0363   return 0;
0364 }
0365
0366 double TMatacq::interpolate(double amplx) {
0367   double T;
0368   int kmax = (int)pkval - firstsample;
0369
0370   int bin_low = 0;
0371   for (Int_t k = 0; k < kmax; k++)
0372     if (0. < bong[k] && bong[k] < amplx) {
0373       bin_low = k;
0374     }
0375   if (bin_low == 0)
0376     return -301.;
0377   int bin_high = 0;
0378   for (Int_t k = kmax; k >= 0; k--)
0379     if (bong[k] > amplx) {
0380       bin_high = k;
0381     }
0382   if (bin_high == 0)
0383     return -302.;
0384   if (bin_high < bin_low)
0385     return -303.;
0386
0387   if (bin_low == bin_high) {
0388     T = (double)bin_high;
0389   } else {
0390     double slope = (bong[bin_high] - bong[bin_low]) / ((double)(bin_high - bin_low));
0391     T = (double)bin_high - (bong[bin_high] - amplx) / slope;
0392   }
0393
0394   return T;
0395 }
0396
0397 void TMatacq::enterdata(Int_t anevt) {
0398   if (anevt < 2 * nevlasers) {
0399     if (anevt < nevlasers) {
0400       nevmtq0++;
0401       status[nevmtq0 - 1] = anevt;
0402       comp_trise[nevmtq0 - 1] = trise;
0403       comp_peak[nevmtq0 - 1] = pkval;
0404     } else {
0405       nevmtq1++;
0406       status[nevmtq0 + nevmtq1 - 1] = anevt;
0407       comp_trise[nevmtq0 + nevmtq1 - 1] = trise;
0408       comp_peak[nevmtq0 + nevmtq1 - 1] = pkval;
0409     }
0410   } else {
0411     if (anevt < 3 * nevlasers) {
0412       nevmtq0++;
0413       status[nevmtq0 - 1] = anevt;
0414       comp_trise[nevmtq0 - 1] = trise;
0415       comp_peak[nevmtq0 - 1] = pkval;
0416     } else {
0417       nevmtq1++;
0418       status[nevmtq0 + nevmtq1 - 1] = anevt;
0419       comp_trise[nevmtq0 + nevmtq1 - 1] = trise;
0420       comp_peak[nevmtq0 + nevmtq1 - 1] = pkval;
0421     }
0422   }
0423 }
0424
0425 void TMatacq::printmatacqData(int gRunNumber, int color, int timestart) {
0426   FILE *fmatacq;
0427   char filename[80];
0428   int i;
0429   double ss;
0430   sprintf(filename, "runMatacq%d.pedestal", gRunNumber);
0431   fmatacq = fopen(filename, "w");
0432   if (fmatacq == nullptr)
0433     printf("Error while opening file : %s\n", filename);
0434
0435   double sumtrise = 0.;
0436   double sumtrise2 = 0.;
0437   int timestop = timestart + 3;
0438   double mintrise = 10000.;
0439   double maxtrise = 0.;
0440   for (i = 0; i < nevmtq0; i++) {
0441     if (comp_trise[i] > maxtrise) {
0442       maxtrise = comp_trise[i];
0443     }
0444     if (comp_trise[i] < mintrise) {
0445       mintrise = comp_trise[i];
0446     }
0447     sumtrise += comp_trise[i];
0448     sumtrise2 += comp_trise[i] * comp_trise[i];
0449   }
0450   meantrise = sumtrise / ((double)nevmtq0);
0451   ss = (sumtrise2 / ((double)nevmtq0) - meantrise * meantrise);
0452   if (ss < 0.)
0453     ss = 0.;
0454   sigtrise = sqrt(ss);
0455   fprintf(
0456       fmatacq, "%d %d %d %d %f %f %f %f\n", nevmtq0, color, timestart, timestop, meantrise, sigtrise, mintrise, maxtrise);
0457
0458   sumtrise = 0.;
0459   sumtrise2 = 0.;
0460   mintrise = 10000.;
0461   maxtrise = 0.;
0462   for (i = nevmtq0; i < nevmtq0 + nevmtq1; i++) {
0463     if (comp_trise[i] > maxtrise) {
0464       maxtrise = comp_trise[i];
0465     }
0466     if (comp_trise[i] < mintrise) {
0467       mintrise = comp_trise[i];
0468     }
0469     sumtrise += comp_trise[i];
0470     sumtrise2 += comp_trise[i] * comp_trise[i];
0471   }
0472   meantrise = sumtrise / ((double)nevmtq1);
0473   ss = (sumtrise2 / ((double)nevmtq1) - meantrise * meantrise);
0474   if (ss < 0.)
0475     ss = 0.;
0476   sigtrise = sqrt(ss);
0477   fprintf(
0478       fmatacq, "%d %d %d %d %f %f %f %f\n", nevmtq1, color, timestart, timestop, meantrise, sigtrise, mintrise, maxtrise);
0479
0480   int iret = fclose(fmatacq);
0481   printf(" Closing file : %d\n", iret);
0482 }
0483
0484 int TMatacq::countBadPulses(int gRunNumber) {
0485   FILE *fmatacq;
0486   char filename[80];
0487   sprintf(filename, "badevtsMatacq%d.dat", gRunNumber);
0488   fmatacq = fopen(filename, "w");
0489   if (fmatacq == nullptr)
0490     printf("Error while opening file : %s\n", filename);
0491
0492   int nevbad = 0;
0493   for (Int_t i = 0; i < nevmtq0 + nevmtq1; i++) {
0494     if (comp_trise[i] < meantrise - 3. * sigtrise) {
0495       nevbad++;
0496       fprintf(fmatacq, "%d \n", status[i]);
0497       status[i] = -1;
0498     }
0499     if (comp_trise[i] > meantrise + 3. * sigtrise) {
0500       nevbad++;
0501       fprintf(fmatacq, "%d \n", status[i]);
0502       status[i] = -1;
0503     }
0504   }
0505
0506   int iret = fclose(fmatacq);
0507   printf(" Closing file : %d\n", iret);
0508   return nevbad;
0509 }
0510
0511 void TMatacq::printitermatacqData(int gRunNumber, int color, int timestart) {
0512   FILE *fmatacq;
0513   char filename[80];
0514   int i;
0515   double ss;
0516   sprintf(filename, "runiterMatacq%d.pedestal", gRunNumber);
0517   fmatacq = fopen(filename, "w");
0518   if (fmatacq == nullptr)
0519     printf("Error while opening file : %s\n", filename);
0520
0521   int nevmtqgood = 0;
0522   double sumtrise = 0.;
0523   double sumtrise2 = 0.;
0524   int timestop = timestart + 3;
0525   double mintrise = 10000.;
0526   double maxtrise = 0.;
0527   for (i = 0; i < nevmtq0; i++) {
0528     if (status[i] >= 0) {
0529       nevmtqgood++;
0530       if (comp_trise[i] > maxtrise) {
0531         maxtrise = comp_trise[i];
0532       }
0533       if (comp_trise[i] < mintrise) {
0534         mintrise = comp_trise[i];
0535       }
0536       sumtrise += comp_trise[i];
0537       sumtrise2 += comp_trise[i] * comp_trise[i];
0538     }
0539   }
0540   meantrise = sumtrise / ((double)nevmtqgood);
0541   ss = (sumtrise2 / ((double)nevmtqgood) - meantrise * meantrise);
0542   if (ss < 0.)
0543     ss = 0.;
0544   sigtrise = sqrt(ss);
0545   fprintf(fmatacq,
0546           "%d %d %d %d %f %f %f %f\n",
0547           nevmtqgood,
0548           color,
0549           timestart,
0550           timestop,
0551           meantrise,
0552           sigtrise,
0553           mintrise,
0554           maxtrise);
0555
0556   nevmtqgood = 0;
0557   sumtrise = 0.;
0558   sumtrise2 = 0.;
0559   mintrise = 10000.;
0560   maxtrise = 0.;
0561   for (i = nevmtq0; i < nevmtq0 + nevmtq1; i++) {
0562     if (status[i] >= 0) {
0563       nevmtqgood++;
0564       if (comp_trise[i] > maxtrise) {
0565         maxtrise = comp_trise[i];
0566       }
0567       if (comp_trise[i] < mintrise) {
0568         mintrise = comp_trise[i];
0569       }
0570       sumtrise += comp_trise[i];
0571       sumtrise2 += comp_trise[i] * comp_trise[i];
0572     }
0573   }
0574   meantrise = sumtrise / ((double)nevmtqgood);
0575   ss = (sumtrise2 / ((double)nevmtqgood) - meantrise * meantrise);
0576   if (ss < 0.)
0577     ss = 0.;
0578   sigtrise = sqrt(ss);
0579   fprintf(fmatacq,
0580           "%d %d %d %d %f %f %f %f\n",
0581           nevmtqgood,
0582           color,
0583           timestart,
0584           timestop,
0585           meantrise,
0586           sigtrise,
0587           mintrise,
0588           maxtrise);
0589
0590   int iret = fclose(fmatacq);
0591   printf(" Closing file : %d\n", iret);
0592 }