Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-08 08:15:47

0001 /** This routine is used to fit the signal line
0002       shape of CMS barrel calorimeter
0003 
0004   The method used is the one described in note LPC 84-30 (Billoir 1984) :
0005     "Methode d'ajustement dans un probleme a parametrisation hierarchisee"
0006   First we read the data file which contains a least 1000 events with sample
0007    data , then we adjust the best function shape to the data in order to get
0008    the general parameters which will be used later to get the maximum and the
0009    time of any signal
0010   This function is only used to get the 2 parameters of the general function */
0011 
0012 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TFParams.h>
0013 #include "TMatrixD.h"
0014 #include "TMath.h"
0015 
0016 #include <iostream>
0017 #include <ctime>
0018 
0019 //ClassImp(TFParams)
0020 
0021 using namespace std;
0022 
0023 TFParams::TFParams(int size, int size_sh) {
0024   //int  sdim = size;
0025   //int plshdim = size_sh;
0026 
0027   for (int i = 0; i < 10; i++) {
0028     for (int j = 0; j < 10; j++) {
0029       weight_matrix[i][j] = 8.;
0030     }
0031   }
0032 }
0033 
0034 double TFParams::fitpj(double **adcval, double *parout, double **db_i, double noise_val, int debug) {
0035 #define dimn 10
0036 #define dimin 10
0037 #define plshdim 300
0038 #define nsamp 10
0039 #define ntrack 500
0040   //#define debug debug1
0041 
0042   // ******************************************************************
0043   // *  Definitions of variables used in the routine
0044   // ******************************************************************
0045 
0046   double a1, a2, a3, b1, b2;
0047   int iter, nevt;
0048   //double errpj[dimmat][dimmat]  ;
0049   double bi[ntrack][2], dbi[ntrack][2];
0050   double zi[ntrack][2];
0051   double par3degre[3];
0052   int ioktk[ntrack], nk, nborn_min = 0, nborn_max = 0;
0053   double cti[ntrack][6], dm1i[ntrack][4];
0054   double par[4], tsig[1];
0055   double amp, delta[nsamp], delta2, fun;
0056   double num_fit_min[ntrack], num_fit_max[ntrack];
0057   int i, j, k, imax[ntrack];
0058 
0059   double ampmax[ntrack], dt, t;
0060   double chi2, chi2s, da1[nsamp], da2[nsamp], db1[nsamp], db2[nsamp];
0061   double chi2tot;
0062   double fact2;
0063   double albet, dtsbeta, variab, alpha, beta;
0064   double unsurs1 /*,unsurs2*/;
0065   //  double fit3 ;
0066   int numb_a, numb_b, numb_x;
0067 
0068   fun = 0;
0069   chi2s = 0;
0070   chi2tot = 0;
0071   matrice DA, DAT, BK, DB, DBT, C, CT, D, DM1, CDM1, CDM1CT, Z, CDM1Z, YK, Y, B, X, XINV, RES2;
0072   matrice A_CROISS, ZMCT;
0073 
0074   double *amplu;
0075   amplu = new double[nsamp];
0076 
0077   parout[0] = 0.;
0078   parout[1] = 0.;
0079   parout[2] = 0.;
0080 
0081   //
0082   //  Initialisation of fit parameters
0083   //
0084 
0085   a1 = a1ini;
0086   a2 = a2ini;
0087   a3 = a3ini;
0088   if (METHODE == 2) {
0089     a2 = a3ini;  // for lastshape BETA is the third parameter ( ... ! )
0090   }
0091   if (debug == 1) {
0092     printf(" ------> __> valeurs de a1 %f a2 %f a3 %f\n", a1, a2, a3);
0093   }
0094   for (i = 0; i < ntrack; i++) {
0095     for (j = 0; j < 2; j++) {
0096       bi[i][j] = (double)0.;
0097       dbi[i][j] = (double)0.;
0098       zi[i][j] = (double)0.;
0099       cti[i][j] = (double)0.;
0100       dm1i[i][j] = (double)0.;
0101     }
0102   }
0103 
0104   numb_a = 2;
0105 
0106   //
0107   //  Matrices initialisation
0108   //
0109 
0110   numb_x = 1;
0111   numb_b = 2;
0112   DA = cree_mat(numb_a, numb_x);
0113   DAT = cree_mat(numb_x, numb_a);
0114   BK = cree_mat_prod(DA, DAT);
0115   DB = cree_mat(numb_b, numb_x);
0116   DBT = cree_mat(numb_x, numb_b);
0117   C = cree_mat(numb_a, numb_b);
0118   CT = cree_mat(numb_b, numb_a);
0119   D = cree_mat_prod(DB, DBT);
0120   DM1 = cree_mat_prod(DB, DBT);
0121   CDM1 = cree_mat_prod(C, DM1);
0122   CDM1CT = cree_mat_prod(CDM1, CT);
0123   Z = cree_mat(numb_b, numb_x);
0124   CDM1Z = cree_mat_prod(CDM1, Z);
0125   YK = cree_mat(numb_a, numb_x);
0126   Y = cree_mat(numb_a, numb_x);
0127   B = cree_mat_prod(DA, DAT);
0128   X = cree_mat_prod(DA, DAT);
0129   XINV = cree_mat_prod(DA, DAT);
0130   RES2 = cree_mat(numb_a, numb_x);
0131   A_CROISS = cree_mat(numb_a, numb_x);
0132   ZMCT = cree_mat(numb_b, numb_x);
0133 
0134   //////////////////////////////
0135   // First Loop on iterations //
0136   /////////////////////////////
0137 
0138   for (iter = 0; iter < 6; iter++) {
0139     chi2tot = 0;
0140 
0141     //
0142     //    Set zeros for general matrices
0143     //
0144 
0145     if (debug == 1) {
0146       printf(" Debut de l'iteration numero %d \n", iter);
0147     }
0148     zero_mat(CDM1Z);
0149     zero_mat(Y);
0150     zero_mat(CDM1CT);
0151     zero_mat(B);
0152     zero_mat(X);
0153     zero_mat(CDM1);
0154 
0155     nk = -1;
0156     if (debug == 1) {
0157       printf(" resultats injectes a iterations %d \n", iter);
0158       printf(" parametre a1 = %f \n", a1);
0159       printf(" parametre a2 = %f \n", a2);
0160       printf(" chi2 du fit chi2s = %f \n", chi2s);
0161 
0162       printf(" value de nevtmax _______________ %d \n", nevtmax);
0163     }
0164 
0165     /////////////////////
0166     //  Loop on events //
0167     /////////////////////
0168 
0169     for (nevt = 0; nevt < nevtmax; nevt++) {
0170       //       B1 = BI[nk,1] est la normalisation du signal
0171       //       B2 = BI[nk,2] ewst le dephasage par rapport a une
0172       //                 fonction centree en zero
0173       //        Nous choisissons ici de demarrer avec les resultats
0174       //            de l'ajustement parabolique mais il faudra bien
0175       //            entendu verifier que cela ne biaise pas le resultat !
0176       //            mise a zero des matrices utilisees dans la boucle
0177 
0178       zero_mat(Z);
0179       zero_mat(YK);
0180       zero_mat(BK);
0181       zero_mat(C);
0182       zero_mat(D);
0183 
0184       nk = nevt;
0185       ampmax[nk] = 0.;
0186       imax[nk] = 0;
0187       for (k = 0; k < 10; k++) {
0188         amplu[k] = adcval[nevt][k];
0189         if (amplu[k] > ampmax[nk]) {
0190           ampmax[nk] = amplu[k];
0191           imax[nk] = k;
0192         }
0193       }
0194 
0195       if (iter == 0) {
0196         // start with degree 3 polynomial ....
0197         //fit3 =polfit(ns ,imax[nk] ,par3degre ,errpj ,amplu ) ;
0198         //  std::cout << "Poly Fit Param  :"<< par3degre[0] <<" "<< par3degre[1]<< std::endl;
0199 
0200         // start with parabol
0201         //fit3 = parab(amplu,4,12,par3degre) ;
0202         /*fit3 =*/parab(amplu, 2, 9, par3degre);
0203         //std::cout << "Parab Fit Param :"<< par3degre[0] <<" "<< par3degre[1]<< std::endl;
0204 
0205         // start with basic initial values
0206         //par3degre[0]= ampmax+ampmax/10. ;
0207         //par3degre[1]= (double)imax[nk]+0.1 ;
0208         //bi[nk][0] = ampmax[nk] ;
0209         //bi[nk][1] = (double)imax[nk] ;
0210 
0211         num_fit_min[nevt] = (double)imax[nk] - (double)nsmin;
0212         num_fit_max[nevt] = (double)imax[nk] + (double)nsmax;
0213 
0214         bi[nk][0] = par3degre[0];
0215         bi[nk][1] = par3degre[1];
0216 
0217         if (debug == 1) {
0218           printf("---------> depart ampmax[%d]=%f   maximum %f tim %f \n", nk, ampmax[nk], bi[nk][0], bi[nk][1]);
0219         }
0220 
0221       } else {
0222         //  in other iterations  :
0223         //   increment bi[][] parameters with bdi[][]
0224         //   calculated in previous
0225         //   iteration
0226 
0227         bi[nk][0] += dbi[nk][0];
0228         bi[nk][1] += dbi[nk][1];
0229 
0230         if (debug == 1) {
0231           printf("iter %d valeur de max %f et norma %f poly 3 \n", iter, bi[nk][1], bi[nk][0]);
0232         }
0233       }
0234 
0235       b1 = bi[nk][0];
0236       b2 = bi[nk][1];
0237 
0238       /////////////////////
0239       // Loop on samples //
0240       /////////////////////
0241 
0242       chi2 = 0.;
0243       ioktk[nk] = 0;
0244       ns = nborn_max - nborn_min + 1;
0245 
0246       for (k = 0; k < 10; k++) {
0247         //
0248         //  calculation of fonction used to fit
0249         //
0250 
0251         dt = (double)k - b2;
0252         t = (double)k;
0253         amp = amplu[k];
0254         if (debug == 1) {
0255           printf(" CHECK sample %f ampli %f \n", t, amp);
0256         }
0257         //unsurs1 = 1./sig_err ;
0258         //unsurs2 = 1./(sig_err*sig_err) ;
0259         //unsurs1 = 0.1 ;
0260         //unsurs2 = 0.01 ;
0261 
0262         unsurs1 = 1. / noise_val;
0263         //unsurs2=(1./noise_val)*(1./noise_val);
0264 
0265         //
0266         // Pulse shape function used: pulseShapepj
0267         //
0268 
0269         nborn_min = (int)num_fit_min[nevt];
0270         nborn_max = (int)num_fit_max[nevt];
0271         if (k < nborn_min || k > nborn_max)
0272           continue;
0273         tsig[0] = (double)k;
0274 
0275         if (METHODE == 2) {
0276           par[0] = b1;
0277           par[1] = b2;
0278           par[2] = a1;
0279           par[3] = a2;
0280           fun = pulseShapepj(tsig, par);
0281         }
0282         if (debug == 1) {
0283           printf(" valeur ampli %f et function %f min %d max %d \n", amp, fun, nsmin, nsmax);
0284           printf("min %f max %f \n", num_fit_min[nevt], num_fit_max[nevt]);
0285         }
0286 
0287         //   we need to determine a1,a2 which are global parameters
0288         //    and  b1, b2 which are parameters for each individual signal:
0289         //        b1, b2 = amplitude and time event by event
0290         //        a1, a2 = alpha and beta global
0291         //    we first begin to calculate the derivatives used in the following calculation
0292 
0293         if (METHODE == 2) {
0294           alpha = a1;
0295           beta = a2;
0296           albet = alpha * beta;
0297           if (dt > -albet) {
0298             variab = (double)1. + dt / albet;
0299             dtsbeta = dt / beta;
0300             db1[k] = unsurs1 * fun / b1;
0301             fact2 = fun;
0302             db2[k] = unsurs1 * fact2 * dtsbeta / (albet * variab);
0303             da1[k] = unsurs1 * fact2 * (log(variab) - dtsbeta / (alpha * variab));
0304             da2[k] = unsurs1 * fact2 * dtsbeta * dtsbeta / (albet * variab);
0305           }
0306         }
0307         delta[k] = (amp - fun) * unsurs1;
0308         if (debug == 1) {
0309           printf(" ------->iter %d valeur de k %d amp %f fun %f delta %f \n", iter, k, amp, fun, delta[k]);
0310           printf(" -----> valeur de k %d delta %f da1 %f da2 %f  \n", k, delta[k], da1[k], da2[k]);
0311         }
0312 
0313         chi2 = chi2 + delta[k] * delta[k];
0314 
0315         if (debug == 1) {
0316           printf(" CHECK chi2 %f deltachi2 %f sample %d iter %d \n", chi2, delta[k] * delta[k], k, iter);
0317         }
0318       }
0319 
0320       /////////////////////////
0321       // End Loop on samples //
0322       /////////////////////////
0323 
0324       double wk1wk2;
0325 
0326       ///////////////////////////
0327       // Start Loop on samples //
0328       ///////////////////////////
0329 
0330       for (int k1 = nborn_min; k1 < nborn_max + 1; k1++) {
0331         wk1wk2 = 1.;
0332         int k2 = k1;
0333 
0334         DA.coeff[0][0] = da1[k1] * wk1wk2;
0335         DA.coeff[1][0] = da2[k1] * wk1wk2;
0336         DAT.coeff[0][0] = da1[k2];
0337         DAT.coeff[0][1] = da2[k2];
0338         DB.coeff[0][0] = db1[k1] * wk1wk2;
0339         DB.coeff[1][0] = db2[k1] * wk1wk2;
0340         DBT.coeff[0][0] = db1[k2];
0341         DBT.coeff[0][1] = db2[k2];
0342 
0343         //  Compute derivative matrix : matrix b[2][2]
0344 
0345         produit_mat_int(DA, DAT, BK);
0346 
0347         //  Compute matrix c[2][2]
0348 
0349         produit_mat_int(DA, DBT, C);
0350 
0351         //  Compute matrix d[2][2]
0352 
0353         produit_mat_int(DB, DBT, D);
0354 
0355         //  Compute matrix y[3] and z[2] depending of delta (amp-fun)
0356 
0357         delta2 = delta[k2];
0358 
0359         somme_mat_int_scale(DA, YK, delta2);
0360         somme_mat_int_scale(DB, Z, delta2);
0361 
0362         ioktk[nk]++;
0363       }
0364 
0365       /////////////////////////
0366       // End Loop on samples //
0367       /////////////////////////
0368 
0369       //  Remove events with a bad shape
0370 
0371       if (ioktk[nk] < 4) {
0372         printf(" event rejected because npamp_used = %d \n", ioktk[nk]);
0373         continue;
0374       }
0375       chi2s = chi2 / (2. + (double)ns + 2.);
0376       chi2tot += chi2s;
0377 
0378       if (debug == 1) {
0379         if (nevt == 198 || nevt == 199) {
0380           std::cout << "adc123 pour l'evt " << nevt << " = " << adcval[nevt][nborn_min] << " = "
0381                     << adcval[nevt][imax[nevt]] << " = " << adcval[nevt][nborn_max] << std::endl;
0382           std::cout << "chi2s  pour l'evt " << nevt << " = " << chi2s << " " << chi2 << " " << ns << "  " << iter
0383                     << std::endl;
0384           std::cout << "chi2tot           " << nevt << " = " << chi2tot << "  " << iter << std::endl;
0385         }
0386       }
0387 
0388       //  Transpose matrix C ---> CT
0389 
0390       transpose_mat(C, CT);
0391 
0392       //  Calculate DM1 (inverse of D matrix 2x2)
0393 
0394       inverse_mat(D, DM1);
0395 
0396       //  Set matrix product c*d in memory in order to compute variations
0397       //   of parameters B at the end of the iteration loop
0398       //   the variations of parameters b are dependant of the variations of
0399       //   parameters da[1],da[2]
0400 
0401       cti[nk][0] = CT.coeff[0][0];
0402       cti[nk][1] = CT.coeff[0][1];
0403       cti[nk][2] = CT.coeff[1][0];
0404       cti[nk][3] = CT.coeff[1][1];
0405 
0406       dm1i[nk][0] = DM1.coeff[0][0];
0407       dm1i[nk][1] = DM1.coeff[0][1];
0408       dm1i[nk][2] = DM1.coeff[1][0];
0409       dm1i[nk][3] = DM1.coeff[1][1];
0410 
0411       zi[nk][0] = Z.coeff[0][0];
0412       zi[nk][1] = Z.coeff[1][0];
0413 
0414       //   Sum the matrix b and y after every event
0415 
0416       for (k = 0; k < numb_a; k++) {
0417         Y.coeff[k][0] += YK.coeff[k][0];
0418       }
0419       somme_mat_int(BK, B);
0420 
0421       //   Calculate c(d-1)
0422 
0423       produit_mat(C, DM1, CDM1);
0424 
0425       // Compute c(d-1)ct
0426 
0427       produit_mat_int(CDM1, CT, CDM1CT);
0428 
0429       // Compute c(d-1)z
0430 
0431       produit_mat_int(CDM1, Z, CDM1Z);
0432     }
0433     /////////////////////////
0434     // End Loop on events //
0435     /////////////////////////
0436 
0437     //  Compute b-cdm1ct
0438 
0439     diff_mat(B, CDM1CT, X);
0440     inverse_mat(X, XINV);
0441     diff_mat(Y, CDM1Z, RES2);
0442 
0443     // Calculation is now easy for da[0] da[1]
0444 
0445     produit_mat(XINV, RES2, A_CROISS);
0446 
0447     //  A la fin, on peut iterer en mesurant l'accroissement a apporter
0448     //    des parametres globaux par la formule db[i] = dm1(z-ct*da[i])
0449 
0450     for (k = 0; k < nk + 1; k++) {
0451       if (METHODE == 2) {
0452         ZMCT.coeff[0][0] = zi[k][0] - (cti[k][0] * A_CROISS.coeff[0][0] + cti[k][1] * A_CROISS.coeff[1][0]);
0453         ZMCT.coeff[1][0] = zi[k][1] - (cti[k][2] * A_CROISS.coeff[0][0] + cti[k][3] * A_CROISS.coeff[1][0]);
0454       }
0455 
0456       dbi[k][0] = dm1i[k][0] * ZMCT.coeff[0][0] + dm1i[k][1] * ZMCT.coeff[1][0];
0457       dbi[k][1] = dm1i[k][2] * ZMCT.coeff[0][0] + dm1i[k][3] * ZMCT.coeff[1][0];
0458       if (debug == 1) {
0459         if (k < 100) {
0460           printf(" variations de b1= %f et b2= %f  \n", dbi[k][0], dbi[k][1]);
0461         }
0462       }
0463       db_i[k][0] = bi[k][0] + dbi[k][0];
0464       db_i[k][1] = bi[k][1] + dbi[k][1];
0465     }
0466 
0467     //   dbi[0] et dbi[1] mesurent les variations a apporter aux
0468     //   parametres du signal
0469 
0470     a1 += A_CROISS.coeff[0][0];
0471     a2 += A_CROISS.coeff[1][0];
0472 
0473     if (debug == 1) {
0474       printf(" CHECK croiss coef0: %f  croiss coef1: %f iter %d \n",
0475              fabs(A_CROISS.coeff[0][0]),
0476              fabs(A_CROISS.coeff[1][0]),
0477              iter);
0478     }
0479     if (fabs(A_CROISS.coeff[0][0]) < 0.001 && fabs(A_CROISS.coeff[1][0]) < 0.001)
0480       break;
0481   }
0482 
0483   /////////////////////////////
0484   //  End Loop on iterations //
0485   /////////////////////////////
0486 
0487   parout[0] = a1;
0488   parout[1] = a2;
0489   parout[2] = a3;
0490   if (debug == 1) {
0491     printf(" resultats trouves au bout de %d iterations \n", iter);
0492     printf(" parametre a1 = %f \n", a1);
0493     printf(" parametre a2 = %f \n", a2);
0494   }
0495 
0496   if (debug == 1) {
0497     std::cout << " Final chi2 / NDOF  :  " << chi2tot / nevtmax << std::endl;
0498     std::cout << " Final (alpha,beta) : (" << a1 << "," << a2 << ")" << std::endl;
0499   }
0500 
0501   return chi2tot / nevtmax;
0502 }
0503 
0504 ////////////////
0505 // End  Fitpj //
0506 ////////////////
0507 
0508 /**************************************************************************/
0509 void TFParams::set_const(int n_samples, int sample_min, int sample_max, double alpha, double beta, int nevtmaximum) {
0510   /*------------------------------------------------------------------------*/
0511   ns = n_samples;
0512   nsmin = sample_min;
0513   nsmax = sample_max;
0514   nevtmax = nevtmaximum;
0515   a1ini = alpha;
0516   a2ini = 0.0;
0517   a3ini = beta;
0518   step_shape = .04;
0519   METHODE = 2;
0520   if (ns > SDIM2)
0521     printf("warning: NbOfsamples exceed maximum\n");
0522 }
0523 void TFParams::produit_mat(matrice A, matrice B, matrice M) {
0524   int i, j, k;
0525   //  resultat du produit A*B = M
0526   if (A.nb_colonnes != B.nb_lignes) {
0527     printf(" Erreur : produit de matrices de tailles incompatibles \n ");
0528     M.coeff = nullptr;
0529     return;
0530   }
0531   M.nb_lignes = A.nb_lignes;
0532   M.nb_colonnes = B.nb_colonnes;
0533   zero_mat(M);
0534   for (i = 0; i < M.nb_lignes; i++) {
0535     for (j = 0; j < M.nb_colonnes; j++) {
0536       for (k = 0; k < A.nb_colonnes; k++) {
0537         M.coeff[i][j] += A.coeff[i][k] * B.coeff[k][j];
0538       }
0539     }
0540   }
0541   return;
0542 }
0543 
0544 void TFParams::produit_mat_int(matrice A, matrice B, matrice M) {
0545   int i, j, k;
0546   if (A.nb_colonnes != B.nb_lignes) {
0547     printf(" Erreur : produit de matrices de tailles incompatibles \n ");
0548     M.coeff = nullptr;
0549     return;
0550   }
0551   M.nb_lignes = A.nb_lignes;
0552   M.nb_colonnes = B.nb_colonnes;
0553   for (i = 0; i < M.nb_lignes; i++) {
0554     for (j = 0; j < M.nb_colonnes; j++) {
0555       for (k = 0; k < A.nb_colonnes; k++) {
0556         M.coeff[i][j] += A.coeff[i][k] * B.coeff[k][j];
0557       }
0558     }
0559   }
0560   return;
0561 }
0562 void TFParams::diff_mat(matrice A, matrice B, matrice M) {
0563   int i, j;
0564   //resultat de la difference A-B = M
0565   if (A.nb_lignes != B.nb_lignes) {
0566     printf(" Erreur : difference de matrices de tailles incompatibles \n ");
0567     M.coeff = nullptr;
0568     return;
0569   }
0570   M.nb_lignes = A.nb_lignes;
0571   M.nb_colonnes = A.nb_colonnes;
0572   for (i = 0; i < M.nb_lignes; i++) {
0573     for (j = 0; j < M.nb_colonnes; j++) {
0574       M.coeff[i][j] = A.coeff[i][j] - B.coeff[i][j];
0575     }
0576   }
0577   return;
0578 }
0579 void TFParams::copie_colonne_mat(matrice A, matrice M, int nk) {
0580   int i, j;
0581   int k;
0582   /* resultat de la copie de A dans un vecteur colonne M */
0583   k = 0;
0584   for (i = 0; i < A.nb_lignes; i++) {
0585     for (j = 0; j < A.nb_colonnes; j++) {
0586       M.coeff[nk][k] = A.coeff[i][j];
0587       printf(" copie nk %d  i %d j %d k %d A %e M %e \n ", nk, i, j, k, A.coeff[i][j], M.coeff[nk][k]);
0588       k++;
0589     }
0590   }
0591   return;
0592 }
0593 
0594 void TFParams::somme_mat_int(matrice A, matrice M) {
0595   int i, j;
0596   /* resultat de la somme integree M += A */
0597   if (A.nb_lignes != M.nb_lignes) {
0598     printf(" Erreur : somme de matrices de tailles incompatibles \n ");
0599     M.coeff = nullptr;
0600     return;
0601   }
0602   M.nb_lignes = A.nb_lignes;
0603   M.nb_colonnes = A.nb_colonnes;
0604   for (i = 0; i < M.nb_lignes; i++) {
0605     for (j = 0; j < M.nb_colonnes; j++)
0606       M.coeff[i][j] += A.coeff[i][j];
0607   }
0608   return;
0609 }
0610 void TFParams::somme_mat_int_scale(matrice A, matrice M, double delta) {
0611   int i, j;
0612   M.nb_lignes = A.nb_lignes;
0613   M.nb_colonnes = A.nb_colonnes;
0614   for (i = 0; i < M.nb_lignes; i++) {
0615     for (j = 0; j < M.nb_colonnes; j++)
0616       M.coeff[i][j] += A.coeff[i][j] * delta;
0617   }
0618   return;
0619 }
0620 void TFParams::transpose_mat(matrice A, matrice M) {
0621   int i, j;
0622   // resultat de la transposition = matrice M
0623   for (i = 0; i < A.nb_lignes; i++) {
0624     for (j = 0; j < A.nb_colonnes; j++) {
0625       M.coeff[j][i] = A.coeff[i][j];
0626     }
0627   }
0628   return;
0629 }
0630 matrice cree_mat_prod(matrice A, matrice B) {
0631   int i, j;
0632   matrice M; /* resultat de la creation */
0633 
0634   M.nb_lignes = A.nb_lignes;
0635   M.nb_colonnes = B.nb_colonnes;
0636   M.coeff = (double **)malloc(M.nb_lignes * sizeof(double *));
0637   for (i = 0; i < M.nb_lignes; i++)
0638     M.coeff[i] = (double *)calloc(M.nb_colonnes, sizeof(double));
0639   for (i = 0; i < M.nb_lignes; i++) {
0640     for (j = 0; j < M.nb_colonnes; j++) {
0641       M.coeff[i][j] = 0.;
0642     }
0643   }
0644   //printf(" creation de matrice ---->  nlignes %d ncolonnes %d  \n",
0645   //     M.nb_lignes,M.nb_colonnes) ;
0646   return (M);
0647 }
0648 matrice cree_mat(int n_lignes, int n_colonnes) {
0649   int i, j;
0650   matrice M; /* resultat de la creation */
0651 
0652   M.nb_lignes = n_lignes;
0653   M.nb_colonnes = n_colonnes;
0654   M.coeff = (double **)malloc(M.nb_lignes * sizeof(double *));
0655   for (i = 0; i < M.nb_lignes; i++)
0656     M.coeff[i] = (double *)calloc(M.nb_colonnes, sizeof(double));
0657   for (i = 0; i < M.nb_lignes; i++) {
0658     for (j = 0; j < M.nb_colonnes; j++) {
0659       M.coeff[i][j] = 0.;
0660     }
0661   }
0662   //printf(" creation de matrice --->  nlignes %d ncolonnes %d  \n",
0663   // M.nb_lignes,M.nb_colonnes) ;
0664   return (M);
0665 }
0666 
0667 void fill_mat(matrice A, matrice M) {
0668   int i, j;
0669   /* on remplit la matrice M avec la matrice A */
0670 
0671   M.nb_lignes = A.nb_lignes;
0672   M.nb_colonnes = A.nb_colonnes;
0673   for (i = 0; i < M.nb_lignes; i++) {
0674     for (j = 0; j < M.nb_colonnes; j++) {
0675       M.coeff[i][j] = A.coeff[i][j];
0676       printf("matrice remplie %e \n", M.coeff[i][j]);
0677     }
0678   }
0679   return;
0680 }
0681 void TFParams::print_mat(matrice M) {
0682   int i, j;
0683   if (M.coeff == nullptr) {
0684     printf(" erreur : affichage d'une matrice vide \n");
0685     return;
0686   }
0687   printf(" m_nli %d M_ncol %d \n", M.nb_lignes, M.nb_colonnes);
0688   for (i = 0; i < M.nb_lignes; i++) {
0689     for (j = 0; j < M.nb_colonnes; j++)
0690       printf(" MATRICE i= %d j= %d ---> %e \n", i, j, M.coeff[i][j]);
0691   }
0692   //printf(" apres passage d'impression \n") ;
0693   return;
0694 }
0695 void TFParams::zero_mat(matrice M) {
0696   int i, j;
0697   for (i = 0; i < M.nb_lignes; i++) {
0698     for (j = 0; j < M.nb_colonnes; j++)
0699       M.coeff[i][j] = 0.;
0700   }
0701   return;
0702 }
0703 void TFParams::zero_mat_nk(matrice M, int nk) {
0704   int j;
0705   for (j = 0; j < M.nb_colonnes; j++)
0706     M.coeff[nk][j] = 0.;
0707   return;
0708 }
0709 void TFParams::print_mat_nk(matrice M, int nk) {
0710   int j;
0711   if (M.coeff == nullptr) {
0712     printf(" erreur : affichage d'une matrice vide \n");
0713   } else {
0714     printf(" nk = %d m_nli %d M_ncol %d \n", nk, M.nb_lignes, M.nb_colonnes);
0715     for (j = 0; j < M.nb_colonnes; j++)
0716       printf(" MATRICE nk= %d j= %d  ---> %e \n", nk, j, M.coeff[nk][j]);
0717     printf(" apres passage d'impression \n");
0718   }
0719   return;
0720 }
0721 void TFParams::inverse_mat(matrice A, matrice M) {
0722   /*   A[ligne][colonne]   B[ligne][colonne]   */
0723   int i, j;
0724   double deter = 0.;
0725   /*  M est la matrice inverse de A */
0726 
0727   if (A.nb_lignes != A.nb_colonnes) {
0728     printf(" attention matrice non inversible !!!! %d lignes %d colonnes \n", A.nb_lignes, A.nb_colonnes);
0729     return;
0730   }
0731   zero_mat(M);
0732   if (A.nb_lignes == 2) {
0733     deter = A.coeff[0][0] * A.coeff[1][1] - A.coeff[0][1] * A.coeff[1][0];
0734     M.coeff[0][0] = A.coeff[1][1] / deter;
0735     M.coeff[0][1] = -A.coeff[0][1] / deter;
0736     M.coeff[1][0] = -A.coeff[1][0] / deter;
0737     M.coeff[1][1] = A.coeff[0][0] / deter;
0738   } else if (A.nb_lignes == 3) {
0739     M.coeff[0][0] = A.coeff[1][1] * A.coeff[2][2] - A.coeff[2][1] * A.coeff[1][2];
0740     M.coeff[1][1] = A.coeff[0][0] * A.coeff[2][2] - A.coeff[2][0] * A.coeff[0][2];
0741 
0742     M.coeff[2][2] = A.coeff[0][0] * A.coeff[1][1] - A.coeff[0][1] * A.coeff[1][0];
0743     M.coeff[0][1] = A.coeff[2][1] * A.coeff[0][2] - A.coeff[0][1] * A.coeff[2][2];
0744     M.coeff[0][2] = A.coeff[0][1] * A.coeff[1][2] - A.coeff[1][1] * A.coeff[0][2];
0745     M.coeff[1][0] = A.coeff[1][2] * A.coeff[2][0] - A.coeff[1][0] * A.coeff[2][2];
0746     M.coeff[1][2] = A.coeff[1][0] * A.coeff[0][2] - A.coeff[0][0] * A.coeff[1][2];
0747     M.coeff[2][0] = A.coeff[1][0] * A.coeff[2][1] - A.coeff[1][1] * A.coeff[2][0];
0748     M.coeff[2][1] = A.coeff[0][1] * A.coeff[2][0] - A.coeff[0][0] * A.coeff[2][1];
0749     deter = A.coeff[0][0] * M.coeff[0][0] + A.coeff[1][0] * M.coeff[0][1] + A.coeff[2][0] * M.coeff[0][2];
0750     for (i = 0; i < 3; i++) {
0751       for (j = 0; j < 3; j++)
0752         M.coeff[i][j] = M.coeff[i][j] / deter;
0753     }
0754   } else {
0755     printf(" Attention , on ne peut inverser la MATRICE %d \n", A.nb_lignes);
0756     return;
0757   }
0758 
0759   return;
0760 }
0761 Double_t TFParams::polfit(Int_t ns, Int_t imax, Double_t par3d[dimout], Double_t errpj[dimmat][dimmat], double *adcpj) {
0762   double val, val2, val3, adfmx[dimmat], parfp3[dimout];
0763   double ius[dimmat] = {0.}, maskp3[dimmat];
0764   double deglib, fit3, tm, h, xki2;
0765   int i, nus, ilow, isup;
0766   val = adcpj[imax];
0767   val2 = val / 2.;
0768   val3 = val / 3.;
0769   ilow = 0;
0770   isup = ns;
0771   deglib = -4.;
0772   for (i = 0; i < ns; i++) {
0773     deglib = deglib + 1.;
0774     ius[i] = 1.;
0775     if ((adcpj[i] < val3) && (i < imax)) {
0776       ilow = i;
0777     }
0778     if (adcpj[i] > val2) {
0779       isup = i;
0780     }
0781   }
0782   ilow = ilow + 1;
0783   if (ilow == imax)
0784     ilow = ilow - 1;
0785   if (isup - ilow < 3)
0786     isup = ilow + 3;
0787   nus = 0;
0788   for (i = ilow; i <= isup; i++) {
0789     adfmx[nus] = adcpj[i];
0790     maskp3[nus] = 0.;
0791     if (ius[i] == 1) {
0792       maskp3[nus] = 1.;
0793       nus = nus + 1;
0794     }
0795   }
0796   if (nus < 4)
0797     return 10000.;
0798   xki2 = f3deg(nus, parfp3, maskp3, adfmx, errpj);
0799   tm = parfp3[4];
0800   h = parfp3[5];
0801   tm = tm + (double)ilow;
0802   par3d[0] = h;
0803   par3d[1] = tm;
0804   fit3 = xki2;
0805   return fit3;
0806 }
0807 double TFParams::f3deg(
0808     int nmxu, double parom[dimout], double mask[dimmat], double adcpj[dimmat], double errpj[dimmat][dimmat]) {
0809   /*                                                                   */
0810   /*  fit   3rd degree polynomial                                      */
0811   /*  nmxu = nb of samples in sample data array adcpj[]
0812     parom   values of parameters
0813     errpj  inverse of the error matrix
0814     fplo3dg uses only the diagonal terms of errpj[][]
0815 */
0816   int i, k, l /*,iworst*/;
0817   double h, t2, tm, delta, tmp;
0818   double xki2, dif, difmx, deglib;
0819   double t[dimmat], f[dimmat][4];
0820   double cov[dimmat][dimmat], bv[4], invcov[dimmat][dimmat], s /*, deter*/;
0821 
0822   deglib = (double)nmxu - 4.;
0823   for (i = 0; i < nmxu; i++) {
0824     t[i] = i;
0825     f[i][0] = 1.;
0826     f[i][1] = t[i];
0827     f[i][2] = t[i] * t[i];
0828     f[i][3] = f[i][2] * t[i];
0829   }
0830   /*   computation of covariance matrix     */
0831   for (k = 0; k < 4; k++) {
0832     for (l = 0; l < 4; l++) {
0833       s = 0.;
0834       for (i = 0; i < nmxu; i++) {
0835         s = s + f[i][k] * f[i][l] * errpj[i][i] * mask[i];
0836       }
0837       cov[k][l] = s;
0838     }
0839     s = 0.;
0840     for (i = 0; i < nmxu; i++) {
0841       s = s + f[i][k] * adcpj[i] * errpj[i][i] * mask[i];
0842     }
0843     bv[k] = s;
0844   }
0845   /*     parameters                          */
0846   /*deter =*/inverpj(4, cov, invcov);
0847   for (k = 0; k < 4; k++) {
0848     s = 0.;
0849     for (l = 0; l < 4; l++) {
0850       s = s + bv[l] * invcov[l][k];
0851     }
0852     parom[k] = s;
0853   }
0854 
0855   if (parom[3] == 0.) {
0856     parom[4] = -1000.;
0857     parom[5] = -1000.;
0858     parom[6] = -1000.;
0859     return 1000000.;
0860   }
0861   /*    worst hit and ki2                    */
0862   xki2 = 0.;
0863   difmx = 0.;
0864   for (i = 0; i < nmxu; i++) {
0865     t2 = t[i] * t[i];
0866     h = parom[0] + parom[1] * t[i] + parom[2] * t2 + parom[3] * t2 * t[i];
0867     dif = (adcpj[i] - h) * mask[i];
0868     if (dif > difmx) {
0869       // iworst=i  ;
0870       difmx = dif;
0871     }
0872   }
0873   if (deglib > 0.5)
0874     xki2 = xki2 / deglib;
0875   /*     amplitude and maximum position                    */
0876   delta = parom[2] * parom[2] - 3. * parom[3] * parom[1];
0877   if (delta > 0.) {
0878     delta = sqrt(delta);
0879     tm = -(delta + parom[2]) / (3. * parom[3]);
0880     tmp = (delta - parom[2]) / (3. * parom[3]);
0881   } else {
0882     parom[4] = -1000.;
0883     parom[5] = -1000.;
0884     parom[6] = -1000.;
0885     return xki2;
0886   }
0887   parom[4] = tm;
0888   parom[5] = parom[0] + parom[1] * tm + parom[2] * tm * tm + parom[3] * tm * tm * tm;
0889   parom[6] = tmp;
0890   // printf("par --------> %f %f %f %f \n",parom[3],parom[2],parom[1],parom[0]);
0891 
0892   return xki2;
0893 }
0894 /*------------------------------------------------------------------*/
0895 
0896 double TFParams::inverpj(int n, double g[dimmat][dimmat], double ginv[dimmat][dimmat]) {
0897   /*                                                                   */
0898   /*  inversion d une matrice symetrique definie positive de taille n  */
0899   /*  J.P. Pansart   Novembre 99                                       */
0900   /*                                                                   */
0901   int i, j, k, jj;
0902   double r, s;
0903   double deter = 0;
0904   double al[dimmat][dimmat], be[dimmat][dimmat];
0905   /*   initialisation                                                  */
0906   for (i = 0; i < n; i++) {
0907     for (j = 0; j < n; j++) {
0908       al[i][j] = 0.;
0909       be[i][j] = 0.;
0910     }
0911   }
0912   /*  decomposition en vecteurs sur une base orthonormee               */
0913   al[0][0] = sqrt(g[0][0]);
0914   for (i = 1; i < n; i++) {
0915     al[i][0] = g[0][i] / al[0][0];
0916     for (j = 1; j <= i; j++) {
0917       s = 0.;
0918       for (k = 0; k <= j - 1; k++) {
0919         s = s + al[i][k] * al[j][k];
0920       }
0921       r = g[i][j] - s;
0922       if (j < i)
0923         al[i][j] = r / al[j][j];
0924       if (j == i)
0925         al[i][j] = sqrt(r);
0926     }
0927   }
0928   /*  inversion de la matrice al                                       */
0929   be[0][0] = 1. / al[0][0];
0930   for (i = 1; i < n; i++) {
0931     be[i][i] = 1. / al[i][i];
0932     for (j = 0; j < i; j++) {
0933       jj = i - j - 1;
0934       s = 0.;
0935       for (k = jj + 1; k <= i; k++) {
0936         s = s + be[i][k] * al[k][jj];
0937       }
0938       be[i][jj] = -s / al[jj][jj];
0939     }
0940   }
0941   /*   calcul de la matrice ginv                                       */
0942   for (i = 0; i < n; i++) {
0943     for (j = 0; j < n; j++) {
0944       s = 0.;
0945       for (k = 0; k < n; k++) {
0946         s = s + be[k][i] * be[k][j];
0947       }
0948       ginv[i][j] = s;
0949       //    if (debug==1){
0950       //printf("valeur de la matrice %d %d %f \n",i,j,ginv[i][j]) ;
0951       //}
0952     }
0953   }
0954   return deter;
0955 }
0956 /*                                                                   */
0957 /*  inversion d une matrice 3x3                                      */
0958 /*                                                                   */
0959 double TFParams::inv3x3(double a[3][3], double b[3][3]) {
0960   /*   a[ligne][colonne]   b[ligne][colonne]   */
0961   int i, j;
0962   double deter = 0.;
0963   b[0][0] = a[1][1] * a[2][2] - a[2][1] * a[1][2];
0964   b[1][1] = a[0][0] * a[2][2] - a[2][0] * a[0][2];
0965   b[2][2] = a[0][0] * a[1][1] - a[0][1] * a[1][0];
0966   printf("a[x][x] %e %e %e %e %e %e %e \n",
0967          a[0][0],
0968          a[1][1],
0969          a[0][1],
0970          a[1][0],
0971          a[0][0] * a[1][1],
0972          a[0][1] * a[1][0],
0973          b[2][2]);
0974   b[0][1] = a[2][1] * a[0][2] - a[0][1] * a[2][2];
0975   b[0][2] = a[0][1] * a[1][2] - a[1][1] * a[0][2];
0976   b[1][0] = a[1][2] * a[2][0] - a[1][0] * a[2][2];
0977   b[1][2] = a[1][0] * a[0][2] - a[0][0] * a[1][2];
0978   b[2][0] = a[1][0] * a[2][1] - a[1][1] * a[2][0];
0979   b[2][1] = a[0][1] * a[2][0] - a[0][0] * a[2][1];
0980   deter = a[0][0] * b[0][0] + a[1][0] * b[0][1] + a[2][0] * b[0][2];
0981   printf(" deter = %e \n", deter);
0982   for (i = 0; i < 3; i++) {
0983     for (j = 0; j < 3; j++) {
0984       printf(" avant division a[3][3] %d %d  %e \n", i, j, a[i][j]);
0985       printf(" avant division b[3][3] %d %d  %e %e \n", i, j, b[i][j], deter);
0986       b[i][j] = b[i][j] / deter;
0987       printf(" valeur de b[3][3] apres division %d %d  %e %e \n", i, j, b[i][j], deter);
0988     }
0989   }
0990   return deter;
0991 }
0992 
0993 double TFParams::pulseShapepj(Double_t *x, Double_t *par) {
0994   Double_t fitfun;
0995   Double_t ped, h, tm, alpha, beta;
0996   Double_t dt, dtsbeta, albet, variab, puiss;
0997   Double_t b1, b2, a1, a2;
0998   b1 = par[0];
0999   b2 = par[1];
1000   a1 = par[2];
1001   a2 = par[3];
1002 
1003   ped = 0.;
1004   h = b1;
1005   tm = b2;
1006   alpha = a1;
1007   beta = a2;
1008   dt = x[0] - tm;
1009   //printf(" par %f %f %f %f dt = %f albet = %f",b1,b2,a1,a2,dt,albet) ;
1010   albet = alpha * beta;
1011   if (albet <= 0)
1012     return ((Double_t)0.);
1013 
1014   if (dt > -albet) {
1015     dtsbeta = dt / beta;
1016     variab = 1. + dt / albet;
1017     puiss = pow(variab, alpha);
1018     fitfun = h * puiss * exp(-dtsbeta) + ped;
1019     //printf(" dt = %f h = %f puiss = %f exp(-dtsbeta) %f \n",dt,h,puiss,
1020     // exp(-dtsbeta)) ;
1021   } else {
1022     fitfun = ped;
1023   }
1024 
1025   return fitfun;
1026 }
1027 
1028 double TFParams::lastShape(Double_t *x, Double_t *par) {
1029   Double_t fitfun;
1030   Double_t alpha, beta;
1031   Double_t dt, alphadt, exponent;
1032   Double_t b1, b2;
1033   b1 = par[0];
1034   b2 = par[1];
1035   alpha = par[2];
1036   beta = par[3];
1037   dt = x[0] - b2;
1038   alphadt = alpha * dt;
1039   exponent = -(alphadt + (exp(-alphadt) - 1.)) / beta;
1040   fitfun = b1 * exp(exponent);
1041   return fitfun;
1042 }
1043 double TFParams::lastShape2(Double_t *x, Double_t *par) {
1044   Double_t fitfun;
1045   Double_t alpha, beta;
1046   Double_t dt, expo1, dt2, exponent;
1047   Double_t b1, b2;
1048   b1 = par[0];
1049   b2 = par[1];
1050   alpha = par[2];
1051   beta = par[3];
1052   dt = x[0] - b2;
1053   expo1 = exp(-beta * dt);
1054   dt2 = dt * dt;
1055   exponent = -(alpha * dt2 + (expo1 - 1.));
1056   fitfun = b1 * exp(exponent);
1057   return fitfun;
1058 }
1059 
1060 Double_t TFParams::pulseShapepj2(Double_t *x, Double_t *par) {
1061   Double_t fitfun;
1062   Double_t ped, h, /*tm,*/ alpha, beta;
1063   Double_t dt, dtsbeta, albet, variab, puiss;
1064   Double_t b1, /*b2,*/ a1, a2;
1065   b1 = par[0];
1066   //b2 = par[1] ;
1067   a1 = par[2];
1068   a2 = par[3];
1069   ped = 0.;
1070   h = b1;
1071   //tm    =  b2 ;
1072   alpha = a1;
1073   beta = a2;
1074   dt = x[0];
1075   albet = alpha * beta;
1076   if (albet <= 0)
1077     return ((Double_t)0.);
1078 
1079   if (dt > -albet) {
1080     dtsbeta = dt / beta;
1081     variab = 1. + dt / albet;
1082     puiss = pow(variab, alpha);
1083     fitfun = h * puiss * exp(-dtsbeta) + ped;
1084   } else {
1085     fitfun = ped;
1086   }
1087 
1088   /*  printf( "fitfun %f %f %f %f, %f %f %f\n", ped, h, tm, alpha, beta, *x, fitfun );  */
1089 
1090   return fitfun;
1091 }
1092 
1093 double TFParams::parab(Double_t ampl[nsamp], Int_t nmin, Int_t nmax, Double_t parout[3]) {
1094   /* Now we calculate the parabolic adjustement in order to get        */
1095   /*    maximum and time max                                           */
1096 
1097   double denom, dt, amp1, amp2, amp3;
1098   double ampmax = 0.;
1099   int imax = 0;
1100   int k;
1101   /*
1102                                                                    */
1103   for (k = nmin; k < nmax; k++) {
1104     if (ampl[k] > ampmax) {
1105       ampmax = ampl[k];
1106       imax = k;
1107     }
1108   }
1109   amp1 = ampl[imax - 1];
1110   amp2 = ampl[imax];
1111   amp3 = ampl[imax + 1];
1112   denom = 2. * amp2 - amp1 - amp3;
1113   /*                                         */
1114   if (denom > 0.) {
1115     dt = 0.5 * (amp3 - amp1) / denom;
1116   } else {
1117     //printf("denom =%f\n",denom)  ;
1118     dt = 0.5;
1119   }
1120   /*                                             */
1121   /* ampmax correspond au maximum d'amplitude parabolique et dt        */
1122   /* decalage en temps par rapport au sample maximum soit k + dt       */
1123 
1124   parout[0] = amp2 + (amp3 - amp1) * dt * 0.25;
1125   parout[1] = (double)imax + dt;
1126   parout[2] = (double)imax;
1127   return denom;
1128 }
1129 
1130 double TFParams::mixShape(Double_t *x, Double_t *par) {
1131   Double_t fitval0, fitval;
1132   Double_t alpha, beta, fact, puiss;
1133   Double_t dt, alpha2dt, exponent;
1134   Double_t b1, b2, alpha2, t;
1135   b1 = par[0];
1136   b2 = par[1];
1137   alpha = par[2];
1138   alpha2 = par[3];
1139   beta = par[4];
1140   //
1141   t = x[0];
1142   dt = x[0] - b2;
1143   //
1144   if (t > 0.) {
1145     fact = t / b2;
1146     puiss = pow(fact, alpha);
1147     fitval0 = puiss * exp(-alpha * dt / b2);
1148   } else {
1149     fitval0 = 1.;
1150   }
1151   dt = x[0] - b2;
1152   alpha2dt = dt * alpha2;
1153   exponent = -(alpha2dt + (exp(-alpha2dt) - 1.)) / beta;
1154   fitval = b1 * fitval0 * exp(exponent);
1155   return fitval;
1156 }
1157 //=========================
1158 // Method computePulseWidth
1159 //=========================
1160 double TFParams::computePulseWidth(int methode, double alpha_here, double beta_here) {
1161   // level of amplitude where we calculate the width ( level = 0.5 if at 50 % )
1162   //   (level = 0.3 if at 30 % )
1163   double level = 0.30;
1164   // fixed parameters
1165   double amplitude = 1.00;
1166   double offset = 7.00;
1167   double amp_max = amplitude;
1168 
1169   // steps in time
1170   double t_min = offset - 4.50;
1171   double t_max = offset + 12.50;
1172 
1173   int t_step_max = 3000;
1174   double delta_t = (double)((t_max - t_min) / t_step_max);
1175 
1176   // Loop over time ( Loop 2  --> get width )
1177   int t_amp_half_flag = 0;
1178   double t_amp_half_min = 999.;
1179   double t_amp_half_max = -999.;
1180 
1181   for (int t_step = 0; t_step < t_step_max; t_step++) {
1182     double t_val = t_min + (double)t_step * delta_t;
1183     double albet = alpha_here * beta_here;
1184     double dt = t_val - offset;
1185     double amp = 0;
1186 
1187     if (methode == 2) {  // electronic function
1188       if ((t_val - offset) > -albet) {
1189         amp = amplitude * TMath::Power((1 + (dt / (alpha_here * beta_here))), alpha_here) *
1190               TMath::Exp(-1.0 * (dt / beta_here));
1191       } else {
1192         amp = 1.;
1193       }
1194     }
1195 
1196     if (amp > (amp_max * level) && t_amp_half_flag == 0) {
1197       t_amp_half_flag = 1;
1198       t_amp_half_min = t_val;
1199     }
1200 
1201     if (amp < (amp_max * level) && t_amp_half_flag == 1) {
1202       t_amp_half_flag = 2;
1203       t_amp_half_max = t_val;
1204     }
1205   }
1206 
1207   // Compute Width
1208   double width = (t_amp_half_max - t_amp_half_min);
1209 
1210   return width;
1211 }