File indexing completed on 2024-04-06 11:57:51
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithFunction.h"
0019
0020 #include <iostream>
0021 #include "TMath.h"
0022
0023
0024
0025
0026 PulseFitWithFunction::PulseFitWithFunction() {
0027 fNsamples = 0;
0028 fNum_samp_bef_max = 0;
0029 fNum_samp_after_max = 0;
0030 }
0031
0032
0033 PulseFitWithFunction::~PulseFitWithFunction() {}
0034
0035
0036
0037 void PulseFitWithFunction::init(int n_samples, int samplb, int sampla, int niter, double alfa, double beta) {
0038
0039
0040
0041
0042
0043
0044 fNsamples = n_samples;
0045 fAlpha_laser = alfa;
0046 fBeta_laser = beta;
0047 fAlpha_beam = 0.98;
0048 fBeta_beam = 2.04;
0049 fNb_iter = niter;
0050 fNum_samp_bef_max = samplb;
0051 fNum_samp_after_max = sampla;
0052
0053
0054
0055
0056 return;
0057 }
0058
0059
0060 double PulseFitWithFunction::doFit(double *adc) {
0061 double parout[4];
0062
0063 double chi2;
0064
0065
0066
0067 Fit_parab(&adc[0], 3, fNsamples, parout);
0068 amp_parab = parout[0];
0069 tim_parab = parout[1];
0070 imax = (int)parout[2];
0071 amp_max = parout[3];
0072
0073 fNumber_samp_max = imax;
0074
0075
0076 if (amp_parab < 1.) {
0077 tim_parab = (double)imax;
0078 amp_parab = amp_max;
0079 }
0080
0081 fValue_tim_max = tim_parab;
0082 fFunc_max = amp_parab;
0083 fTim_max = tim_parab;
0084
0085
0086
0087 chi2 = Fit_electronic(0, &adc[0], 8.);
0088
0089
0090
0091
0092
0093
0094
0095
0096 return chi2;
0097 }
0098
0099
0100
0101
0102
0103 double PulseFitWithFunction::Fit_electronic(int data, double *adc_to_fit, double sigmas_sample) {
0104
0105
0106
0107
0108
0109 double chi2 = 0;
0110 double d_alpha, d_beta;
0111
0112 fAlpha = fAlpha_laser;
0113 fBeta = fBeta_laser;
0114 if (data == 1) {
0115 fAlpha = fAlpha_beam;
0116 fBeta = fBeta_beam;
0117 }
0118
0119 fAmp_fitted_max = 0.;
0120 fTim_fitted_max = 0.;
0121 double un_sur_sigma = 1.;
0122 double variation_func_max = 0.;
0123 double variation_tim_max = 0.;
0124
0125 if (fValue_tim_max > 20. || fValue_tim_max < 3.) {
0126 fValue_tim_max = fNumber_samp_max;
0127 }
0128 int num_fit_min = (int)(fValue_tim_max - fNum_samp_bef_max);
0129 int num_fit_max = (int)(fValue_tim_max + fNum_samp_after_max);
0130
0131
0132 if (sigmas_sample > 0.)
0133 un_sur_sigma = 1. / sigmas_sample;
0134
0135 double func, delta;
0136
0137 for (int iter = 0; iter < fNb_iter; iter++) {
0138
0139 chi2 = 0.;
0140 double d11 = 0.;
0141 double d12 = 0.;
0142 double d22 = 0.;
0143 double z1 = 0.;
0144 double z2 = 0.;
0145 fFunc_max += variation_func_max;
0146 fTim_max += variation_tim_max;
0147 int nsamp_used = 0;
0148
0149
0150
0151
0152 for (int i = num_fit_min; i < num_fit_max + 1; i++) {
0153
0154
0155 func = Electronic_shape((double)i);
0156
0157
0158 double dt = (double)i - fTim_max;
0159 double alpha_beta = fAlpha * fBeta;
0160
0161 if (dt > -alpha_beta) {
0162 double dt_sur_beta = dt / fBeta;
0163
0164 double variable = (double)1. + dt / alpha_beta;
0165 double expo = TMath::Exp(-dt_sur_beta);
0166
0167 double puissance = TMath::Power(variable, fAlpha);
0168 d_alpha = un_sur_sigma * puissance * expo;
0169 d_beta = fFunc_max * d_alpha * dt_sur_beta / (alpha_beta * variable);
0170 } else {
0171 continue;
0172 }
0173
0174 nsamp_used++;
0175
0176 d11 += d_alpha * d_alpha;
0177 d12 += d_alpha * d_beta;
0178 d22 += d_beta * d_beta;
0179
0180 delta = (adc_to_fit[i] - func) * un_sur_sigma;
0181
0182 z1 += delta * d_alpha;
0183 z2 += delta * d_beta;
0184 chi2 += delta * delta;
0185 }
0186 double denom = d11 * d22 - d12 * d12;
0187 if (denom == 0.) {
0188
0189 return 101;
0190 }
0191 if (nsamp_used < 3) {
0192
0193 return 102;
0194 }
0195
0196 variation_func_max = (z1 * d22 - z2 * d12) / denom;
0197 variation_tim_max = (-z1 * d12 + z2 * d11) / denom;
0198 chi2 = chi2 / ((double)nsamp_used - 2.);
0199 }
0200
0201 fAmp_fitted_max = fFunc_max + variation_func_max;
0202 fTim_fitted_max = fTim_max + variation_tim_max;
0203
0204 return chi2;
0205 }
0206
0207
0208
0209
0210 double PulseFitWithFunction::Electronic_shape(double tim) {
0211
0212 double func_electronic, dtsbeta, variable, puiss;
0213 double albet = fAlpha * fBeta;
0214 if (albet <= 0)
0215 return ((Double_t)0.);
0216 double dt = tim - fTim_max;
0217 if (dt > -albet) {
0218 dtsbeta = dt / fBeta;
0219 variable = 1. + dt / albet;
0220 puiss = TMath::Power(variable, fAlpha);
0221 func_electronic = fFunc_max * puiss * TMath::Exp(-dtsbeta);
0222 } else
0223 func_electronic = 0.;
0224
0225 return func_electronic;
0226 }
0227 void PulseFitWithFunction::Fit_parab(Double_t *ampl, Int_t nmin, Int_t nmax, Double_t *parout) {
0228
0229
0230
0231 double denom, dt, amp1, amp2, amp3;
0232 double ampmax = 0.;
0233 int imax = 0;
0234 int k;
0235
0236
0237 for (k = nmin; k < nmax; k++) {
0238
0239 if (ampl[k] > ampmax) {
0240 ampmax = ampl[k];
0241 imax = k;
0242 }
0243 }
0244 amp1 = ampl[imax - 1];
0245 amp2 = ampl[imax];
0246 amp3 = ampl[imax + 1];
0247 denom = 2. * amp2 - amp1 - amp3;
0248
0249 if (denom > 0.) {
0250 dt = 0.5 * (amp3 - amp1) / denom;
0251 } else {
0252
0253 dt = 0.5;
0254 }
0255
0256
0257
0258
0259 parout[0] = amp2 + (amp3 - amp1) * dt * 0.25;
0260 parout[1] = (double)imax + dt;
0261 parout[2] = (double)imax;
0262 parout[3] = ampmax;
0263 return;
0264 }