File indexing completed on 2024-04-06 11:57:52
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
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
0020
0021 using namespace std;
0022
0023 TFParams::TFParams(int size, int size_sh) {
0024
0025
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
0041
0042
0043
0044
0045
0046 double a1, a2, a3, b1, b2;
0047 int iter, nevt;
0048
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 ;
0065
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
0083
0084
0085 a1 = a1ini;
0086 a2 = a2ini;
0087 a3 = a3ini;
0088 if (METHODE == 2) {
0089 a2 = a3ini;
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
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
0136
0137
0138 for (iter = 0; iter < 6; iter++) {
0139 chi2tot = 0;
0140
0141
0142
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
0167
0168
0169 for (nevt = 0; nevt < nevtmax; nevt++) {
0170
0171
0172
0173
0174
0175
0176
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
0197
0198
0199
0200
0201
0202 parab(amplu, 2, 9, par3degre);
0203
0204
0205
0206
0207
0208
0209
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
0223
0224
0225
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
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
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
0258
0259
0260
0261
0262 unsurs1 = 1. / noise_val;
0263
0264
0265
0266
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
0288
0289
0290
0291
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
0322
0323
0324 double wk1wk2;
0325
0326
0327
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
0344
0345 produit_mat_int(DA, DAT, BK);
0346
0347
0348
0349 produit_mat_int(DA, DBT, C);
0350
0351
0352
0353 produit_mat_int(DB, DBT, D);
0354
0355
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
0367
0368
0369
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
0389
0390 transpose_mat(C, CT);
0391
0392
0393
0394 inverse_mat(D, DM1);
0395
0396
0397
0398
0399
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
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
0422
0423 produit_mat(C, DM1, CDM1);
0424
0425
0426
0427 produit_mat_int(CDM1, CT, CDM1CT);
0428
0429
0430
0431 produit_mat_int(CDM1, Z, CDM1Z);
0432 }
0433
0434
0435
0436
0437
0438
0439 diff_mat(B, CDM1CT, X);
0440 inverse_mat(X, XINV);
0441 diff_mat(Y, CDM1Z, RES2);
0442
0443
0444
0445 produit_mat(XINV, RES2, A_CROISS);
0446
0447
0448
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
0468
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
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
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
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
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
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
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
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;
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
0645
0646 return (M);
0647 }
0648 matrice cree_mat(int n_lignes, int n_colonnes) {
0649 int i, j;
0650 matrice M;
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
0663
0664 return (M);
0665 }
0666
0667 void fill_mat(matrice A, matrice M) {
0668 int i, j;
0669
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
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
0723 int i, j;
0724 double deter = 0.;
0725
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
0811
0812
0813
0814
0815
0816 int i, k, l ;
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 ;
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
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
0846 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
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
0870 difmx = dif;
0871 }
0872 }
0873 if (deglib > 0.5)
0874 xki2 = xki2 / deglib;
0875
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
0891
0892 return xki2;
0893 }
0894
0895
0896 double TFParams::inverpj(int n, double g[dimmat][dimmat], double ginv[dimmat][dimmat]) {
0897
0898
0899
0900
0901 int i, j, k, jj;
0902 double r, s;
0903 double deter = 0;
0904 double al[dimmat][dimmat], be[dimmat][dimmat];
0905
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
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
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
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
0950
0951
0952 }
0953 }
0954 return deter;
0955 }
0956
0957
0958
0959 double TFParams::inv3x3(double a[3][3], double b[3][3]) {
0960
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
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
1020
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, alpha, beta;
1063 Double_t dt, dtsbeta, albet, variab, puiss;
1064 Double_t b1, a1, a2;
1065 b1 = par[0];
1066
1067 a1 = par[2];
1068 a2 = par[3];
1069 ped = 0.;
1070 h = b1;
1071
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
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
1095
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
1118 dt = 0.5;
1119 }
1120
1121
1122
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
1159
1160 double TFParams::computePulseWidth(int methode, double alpha_here, double beta_here) {
1161
1162
1163 double level = 0.30;
1164
1165 double amplitude = 1.00;
1166 double offset = 7.00;
1167 double amp_max = amplitude;
1168
1169
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
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) {
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
1208 double width = (t_amp_half_max - t_amp_half_min);
1209
1210 return width;
1211 }