File indexing completed on 2024-04-06 11:57:52
0001
0002
0003
0004
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
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
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
0055 TMatacq::~TMatacq() {}
0056
0057 int TMatacq::rawPulseAnalysis(Int_t Nsamp, Double_t *adc)
0058 {
0059 using namespace std;
0060
0061
0062
0063 int k, ithr;
0064 double dsum = 0., dsum2 = 0.;
0065
0066
0067 init();
0068
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
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
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 }