233 double ampMaxAlphaBeta = 0;
234 double tMaxAlphaBeta = 5;
235 double tMaxErrorAlphaBeta = 999;
236 double tMaxRatio = 5;
237 double tMaxErrorRatio = 999;
257 NullChi2 = (sumAA - sumA * sumA / sum1) / sum0;
269 double alphabeta = amplitudeFitParameters[0] * amplitudeFitParameters[1];
270 double invalphabeta = 1. / alphabeta;
271 double alpha = amplitudeFitParameters[0];
272 double beta = amplitudeFitParameters[1];
274 Ratio ratios_[C::MAXSAMPLES * (C::MAXSAMPLES - 1) / 2];
275 unsigned int ratios_size = 0;
291 for (
unsigned int j = i + 1; j <
amplitudes_.size(); j++) {
302 double err1 = Rtmp * Rtmp * (relErr2[
i] + relErr2[j]);
311 (invampl[j] * invampl[j]);
316 double err3 = (0.289 * 0.289) * (invampl[j] * invampl[j]);
318 double totalError =
std::sqrt(err1 + err2 + err3);
321 if (totalError < 1.0 && Rtmp > 0.001 && Rtmp < Rlim[j - i]) {
322 Ratio currentRatio = {
i, (j -
i), Rtmp, totalError };
323 ratios_[ratios_size++] = currentRatio;
330 if (0 == ratios_size)
return;
333 Tmax timesAB_[C::MAXSAMPLES * (C::MAXSAMPLES - 1) / 2];
334 unsigned int timesAB_size = 0;
339 for (
unsigned int i = 0; i < ratios_size; i++) {
341 double stepOverBeta = double(ratios_[i].
step) /
beta;
342 double offset = double(ratios_[i].
index) + alphabeta;
344 double Rmin = ratios_[
i].value - ratios_[
i].error;
345 if (Rmin < 0.001) Rmin = 0.001;
347 double Rmax = ratios_[
i].value + ratios_[
i].error;
348 double RLimit = Rlim[ratios_[
i].step];
349 if (Rmax > RLimit) Rmax = RLimit;
363 double tmax = 0.5 * (time1 + time2);
364 double tmaxerr = 0.5 *
std::sqrt((time1 - time2) * (time1 - time2));
369 int itmin =
std::max(-1,
int(std::floor(tmax - alphabeta)));
370 double loffset = (double(itmin) -
tmax) * invalphabeta;
371 for (
unsigned int it = itmin + 1; it <
amplitudes_.size(); it++) {
372 loffset += invalphabeta;
376 double term1 = 1.0 + loffset;
383 sumff += f * (f * inverr2);
389 chi2 = sumAA - sumAf * (sumAf / sumff);
390 amp = (sumAf / sumff);
396 if (chi2 > 0 && tmaxerr > 0 && tmax > 0) {
397 Tmax currentTmax = { ratios_[
i].index, ratios_[
i].step,
tmax, tmaxerr,
399 timesAB_[timesAB_size++] = currentTmax;
404 if (0 == timesAB_size)
return;
407 double chi2min = 1.0e+9;
410 for (
unsigned int i = 0; i < timesAB_size; i++) {
411 if (timesAB_[i].chi2 < chi2min) {
412 chi2min = timesAB_[
i].chi2;
420 double chi2Limit = chi2min + 1.0;
423 for (
unsigned int i = 0; i < timesAB_size; i++) {
424 if (timesAB_[i].chi2 < chi2Limit) {
425 double inverseSigmaSquared =
426 1.0 / (timesAB_[
i].error * timesAB_[
i].error);
427 time_wgt += inverseSigmaSquared;
428 time_max += timesAB_[
i].value * inverseSigmaSquared;
432 tMaxAlphaBeta = time_max / time_wgt;
433 tMaxErrorAlphaBeta = 1.0 /
sqrt(time_wgt);
438 for (
unsigned int i = 0; i <
amplitudes_.size(); i++) {
441 double offset = (double(i) - tMaxAlphaBeta) * invalphabeta;
442 double term1 = 1.0 +
offset;
446 sumff += f * (f * inverr2);
451 ampMaxAlphaBeta = sumAf / sumff;
452 double chi2AlphaBeta = (sumAA - sumAf * sumAf / sumff) / sum0;
453 if (chi2AlphaBeta > NullChi2) {
481 for (
unsigned int i = 0; i < ratios_size; i++) {
483 if (ratios_[i].
step == 1 && ratios_[i].
value >= timeFitLimits.first &&
484 ratios_[i].value <= timeFitLimits.second) {
486 double time_max_i = ratios_[
i].index;
490 double u = timeFitParameters[timeFitParameters.size() - 1];
491 for (
int k = timeFitParameters.size() - 2;
k >= 0;
k--) {
492 u = u * ratios_[
i].value + timeFitParameters[
k];
496 double du = (timeFitParameters.size() - 1) *
497 timeFitParameters[timeFitParameters.size() - 1];
498 for (
int k = timeFitParameters.size() - 2;
k >= 1;
k--) {
499 du = du * ratios_[
i].value +
k * timeFitParameters[
k];
503 double errorsquared = ratios_[
i].error * ratios_[
i].error * du * du;
504 if (errorsquared > 0) {
506 time_max += (time_max_i - u) / errorsquared;
507 time_wgt += 1.0 / errorsquared;
519 tMaxRatio = time_max / time_wgt;
520 tMaxErrorRatio = 1.0 /
sqrt(time_wgt);
534 (tMaxAlphaBeta * (10.0 - ampMaxAlphaBeta /
ampMaxError_) +
535 tMaxRatio * (ampMaxAlphaBeta /
ampMaxError_ - 5.0)) / 5.0;
537 (tMaxErrorAlphaBeta * (10.0 - ampMaxAlphaBeta /
ampMaxError_) +
538 tMaxErrorRatio * (ampMaxAlphaBeta /
ampMaxError_ - 5.0)) / 5.0;
std::array< double, C::MAXSAMPLES > amplitudeErrors_
CalculatedRecHit calculatedRechit_
std::array< double, C::MAXSAMPLES > amplitudeIE2_
std::array< double, C::MAXSAMPLES > amplitudes_
static const double tmax[3]
std::array< bool, C::MAXSAMPLES > useless_
const double Rmax[kNumberCalorimeter]
const double Rmin[kNumberCalorimeter]