Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:39 

0001 #ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRatioMethodAlgo_HH
0002 #define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRatioMethodAlgo_HH
0003 
0004 /** \class EcalUncalibRecHitRatioMethodAlgo
0005  *  Template used to compute amplitude, pedestal, time jitter, chi2 of a pulse
0006  *  using a ratio method
0007  *
0008  *  \author A. Ledovskoy (Design) - M. Balazs (Implementation)
0009  */
0010 
0011 #include "Math/SVector.h"
0012 #include "Math/SMatrix.h"
0013 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
0014 #include "CondFormats/EcalObjects/interface/EcalSampleMask.h"
0015 #include <vector>
0016 #include <array>
0017 
0018 //#include "vdt/vdtMath.h"
0019 #include "DataFormats/Math/interface/approx_exp.h"
0020 #include "DataFormats/Math/interface/approx_log.h"
0021 
0022 #define RANDOM_MAGIC
0023 
0024 #include <random>
0025 
0026 namespace myMath {
0027   inline float fast_expf(float x) { return unsafe_expf<6>(x); }
0028   inline float fast_logf(float x) { return unsafe_logf<7>(x); }
0029 }  // namespace myMath
0030 
0031 template <class C>
0032 class EcalUncalibRecHitRatioMethodAlgo {
0033 public:
0034   struct Ratio {
0035     unsigned int index;
0036     unsigned int step;
0037     double value;
0038     double error;
0039   };
0040   struct Tmax {
0041     unsigned int index;
0042     unsigned int step;
0043     double value;
0044     double error;
0045     double amplitude;
0046     double chi2;
0047   };
0048   struct CalculatedRecHit {
0049     double amplitudeMax;
0050     double timeMax;
0051     double timeError;
0052     double chi2;
0053   };
0054 
0055   EcalUncalibratedRecHit makeRecHit(const C &dataFrame,
0056                                     const EcalSampleMask &sampleMask,
0057                                     const double *pedestals,
0058                                     const double *pedestalRMSes,
0059                                     const double *gainRatios,
0060                                     std::vector<double> &timeFitParameters,
0061                                     std::vector<double> &amplitudeFitParameters,
0062                                     std::pair<double, double> &timeFitLimits);
0063 
0064   // more function to be able to compute
0065   // amplitude and time separately
0066   void init(const C &dataFrame,
0067             const EcalSampleMask &sampleMask,
0068             const double *pedestals,
0069             const double *pedestalRMSes,
0070             const double *gainRatios);
0071   void computeTime(std::vector<double> &timeFitParameters,
0072                    std::pair<double, double> &timeFitLimits,
0073                    std::vector<double> &amplitudeFitParameters);
0074   void computeAmplitude(std::vector<double> &amplitudeFitParameters);
0075   CalculatedRecHit getCalculatedRecHit() { return calculatedRechit_; }
0076   bool fixMGPAslew(const C &dataFrame);
0077 
0078   double computeAmplitudeImpl(std::vector<double> &amplitudeFitParameters, double, double);
0079 
0080 protected:
0081   EcalSampleMask sampleMask_;
0082   DetId theDetId_;
0083   std::array<double, C::MAXSAMPLES> amplitudes_;
0084   std::array<double, C::MAXSAMPLES> amplitudeErrors_;
0085   std::array<double, C::MAXSAMPLES> amplitudeIE2_;
0086   std::array<bool, C::MAXSAMPLES> useless_;
0087 
0088   double pedestal_;
0089   int num_;
0090   double ampMaxError_;
0091 
0092   CalculatedRecHit calculatedRechit_;
0093 };
0094 
0095 template <class C>
0096 void EcalUncalibRecHitRatioMethodAlgo<C>::init(const C &dataFrame,
0097                                                const EcalSampleMask &sampleMask,
0098                                                const double *pedestals,
0099                                                const double *pedestalRMSes,
0100                                                const double *gainRatios) {
0101   sampleMask_ = sampleMask;
0102   theDetId_ = DetId(dataFrame.id().rawId());
0103 
0104   calculatedRechit_.timeMax = 5;
0105   calculatedRechit_.amplitudeMax = 0;
0106   calculatedRechit_.timeError = -999;
0107 
0108   // to obtain gain 12 pedestal:
0109   // -> if it's in gain 12, use first sample
0110   // --> average it with second sample if in gain 12 and 3-sigma-noise
0111   // compatible (better LF noise cancellation)
0112   // -> else use pedestal from database
0113   pedestal_ = 0;
0114   num_ = 0;
0115   if (dataFrame.sample(0).gainId() == 1 && sampleMask_.useSample(0, theDetId_)) {
0116     pedestal_ += double(dataFrame.sample(0).adc());
0117     num_++;
0118   }
0119   if (num_ != 0 && dataFrame.sample(1).gainId() == 1 && sampleMask_.useSample(1, theDetId_) &&
0120       std::abs(dataFrame.sample(1).adc() - dataFrame.sample(0).adc()) < 3 * pedestalRMSes[0]) {
0121     pedestal_ += double(dataFrame.sample(1).adc());
0122     num_++;
0123   }
0124   if (num_ != 0)
0125     pedestal_ /= num_;
0126   else
0127     pedestal_ = pedestals[0];
0128 
0129   // fill vector of amplitudes, pedestal subtracted and vector
0130   // of amplitude uncertainties Also, find the uncertainty of a
0131   // sample with max amplitude. We will use it later.
0132 
0133   ampMaxError_ = 0;
0134   double ampMaxValue = -1000;
0135 
0136   // ped-subtracted and gain-renormalized samples. It is VERY
0137   // IMPORTANT to have samples one clock apart which means to
0138   // have vector size equal to MAXSAMPLES
0139   double sample;
0140   double sampleError;
0141   int GainId;
0142   for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
0143     GainId = dataFrame.sample(iSample).gainId();
0144 
0145     bool bad = false;
0146     // only use normally samples which are desired; if sample not to be used
0147     // inflate error so won't generate ratio considered for the measurement
0148     if (!sampleMask_.useSample(iSample, theDetId_)) {
0149       sample = 0;
0150       sampleError = 0;
0151       bad = true;
0152     } else if (GainId == 1) {
0153       sample = double(dataFrame.sample(iSample).adc() - pedestal_);
0154       sampleError = pedestalRMSes[0];
0155     } else if (GainId == 2 || GainId == 3) {
0156       sample = (double(dataFrame.sample(iSample).adc() - pedestals[GainId - 1])) * gainRatios[GainId - 1];
0157       sampleError = pedestalRMSes[GainId - 1] * gainRatios[GainId - 1];
0158     } else {
0159       sample = 0;       // GainId=0 case falls here, from saturation
0160       sampleError = 0;  // inflate error so won't generate ratio considered for
0161                         // the measurement
0162       bad = true;
0163     }
0164 
0165     useless_[iSample] = (sampleError <= 0) | bad;
0166     amplitudes_[iSample] = sample;
0167     // inflate error for useless samples
0168     amplitudeErrors_[iSample] = useless_[iSample] ? 1e+9 : sampleError;
0169     amplitudeIE2_[iSample] = useless_[iSample] ? 0 : 1 / (amplitudeErrors_[iSample] * amplitudeErrors_[iSample]);
0170     if (sampleError > 0) {
0171       if (ampMaxValue < sample) {
0172         ampMaxValue = sample;
0173         ampMaxError_ = sampleError;
0174       }
0175     }
0176   }
0177 }
0178 template <class C>
0179 bool EcalUncalibRecHitRatioMethodAlgo<C>::fixMGPAslew(const C &dataFrame) {
0180   // This fuction finds sample(s) preceeding gain switching and
0181   // inflates errors on this sample, therefore, making this sample
0182   // invisible for Ratio Method. Only gain switching DOWN is
0183   // considered Only gainID=1,2,3 are considered. In case of the
0184   // saturation (gainID=0), we keep "distorted" sample because it is
0185   // the only chance to make time measurement; the qualilty of it will
0186   // be bad anyway.
0187 
0188   bool result = false;
0189 
0190   int GainIdPrev;
0191   int GainIdNext;
0192   for (int iSample = 1; iSample < C::MAXSAMPLES; iSample++) {
0193     // only use samples which are desired
0194     if (!sampleMask_.useSample(iSample, theDetId_))
0195       continue;
0196 
0197     GainIdPrev = dataFrame.sample(iSample - 1).gainId();
0198     GainIdNext = dataFrame.sample(iSample).gainId();
0199     if (GainIdPrev >= 1 && GainIdPrev <= 3 && GainIdNext >= 1 && GainIdNext <= 3 && GainIdPrev < GainIdNext) {
0200       amplitudes_[iSample - 1] = 0;
0201       amplitudeErrors_[iSample - 1] = 1e+9;
0202       amplitudeIE2_[iSample - 1] = 0;
0203       useless_[iSample - 1] = true;
0204       result = true;
0205     }
0206   }
0207   return result;
0208 }
0209 
0210 template <class C>
0211 void EcalUncalibRecHitRatioMethodAlgo<C>::computeTime(std::vector<double> &timeFitParameters,
0212                                                       std::pair<double, double> &timeFitLimits,
0213                                                       std::vector<double> &amplitudeFitParameters) {
0214   //////////////////////////////////////////////////////////////
0215   //                                                          //
0216   //              RATIO METHOD FOR TIME STARTS HERE           //
0217   //                                                          //
0218   //////////////////////////////////////////////////////////////
0219   double ampMaxAlphaBeta = 0;
0220   double tMaxAlphaBeta = 5;
0221   double tMaxErrorAlphaBeta = 999;
0222   double tMaxRatio = 5;
0223   double tMaxErrorRatio = 999;
0224 
0225   double sumAA = 0;
0226   double sumA = 0;
0227   double sum1 = 0;
0228   double sum0 = 0;
0229   double sumAf = 0;
0230   double sumff = 0;
0231   double NullChi2 = 0;
0232 
0233   // null hypothesis = no pulse, pedestal only
0234   for (unsigned int i = 0; i < amplitudes_.size(); i++) {
0235     if (useless_[i])
0236       continue;
0237     double inverr2 = amplitudeIE2_[i];
0238     sum0 += 1;
0239     sum1 += inverr2;
0240     sumA += amplitudes_[i] * inverr2;
0241     sumAA += amplitudes_[i] * (amplitudes_[i] * inverr2);
0242   }
0243   if (sum0 > 0) {
0244     NullChi2 = (sumAA - sumA * sumA / sum1) / sum0;
0245   } else {
0246     // not enough samples to reconstruct the pulse
0247     return;
0248   }
0249 
0250   // Make all possible Ratio's based on any pair of samples i and j
0251   // (j>i) with positive amplitudes_
0252   //
0253   //       Ratio[k] = Amp[i]/Amp[j]
0254   //       where Amp[i] is pedestal subtracted ADC value in a time sample [i]
0255   //
0256   double alphabeta = amplitudeFitParameters[0] * amplitudeFitParameters[1];
0257   double invalphabeta = 1. / alphabeta;
0258   double alpha = amplitudeFitParameters[0];
0259   double beta = amplitudeFitParameters[1];
0260 
0261   Ratio ratios_[C::MAXSAMPLES * (C::MAXSAMPLES - 1) / 2];
0262   unsigned int ratios_size = 0;
0263 
0264   double Rlim[amplitudes_.size()];
0265   for (unsigned int k = 1; k != amplitudes_.size(); ++k)
0266     Rlim[k] = myMath::fast_expf(double(k) / beta) - 0.001;
0267 
0268   double relErr2[amplitudes_.size()];
0269   double invampl[amplitudes_.size()];
0270   for (unsigned int i = 0; i < amplitudes_.size(); i++) {
0271     invampl[i] = (useless_[i]) ? 0 : 1. / amplitudes_[i];
0272     relErr2[i] = (useless_[i]) ? 0 : (amplitudeErrors_[i] * invampl[i]) * (amplitudeErrors_[i] * invampl[i]);
0273   }
0274 
0275   for (unsigned int i = 0; i < amplitudes_.size() - 1; i++) {
0276     if (useless_[i])
0277       continue;
0278     for (unsigned int j = i + 1; j < amplitudes_.size(); j++) {
0279       if (useless_[j])
0280         continue;
0281 
0282       if (amplitudes_[i] > 1 && amplitudes_[j] > 1) {
0283         // ratio
0284         double Rtmp = amplitudes_[i] / amplitudes_[j];
0285 
0286         // error^2 due to stat fluctuations of time samples
0287         // (uncorrelated for both samples)
0288 
0289         double err1 = Rtmp * Rtmp * (relErr2[i] + relErr2[j]);
0290 
0291         // error due to fluctuations of pedestal (common to both samples)
0292         double stat;
0293         if (num_ > 0)
0294           stat = num_;  // num presampeles used to compute pedestal
0295         else
0296           stat = 1;  // pedestal from db
0297         double err2 =
0298             amplitudeErrors_[j] * (amplitudes_[i] - amplitudes_[j]) * (invampl[j] * invampl[j]);  // /sqrt(stat);
0299         err2 *= err2 / stat;
0300 
0301         //error due to integer round-down. It is relevant to low
0302         //amplitudes_ in gainID=1 and negligible otherwise.
0303         double err3 = (0.289 * 0.289) * (invampl[j] * invampl[j]);
0304 
0305         double totalError = std::sqrt(err1 + err2 + err3);
0306 
0307         // don't include useless ratios
0308         if (totalError < 1.0 && Rtmp > 0.001 && Rtmp < Rlim[j - i]) {
0309           Ratio currentRatio = {i, (j - i), Rtmp, totalError};
0310           ratios_[ratios_size++] = currentRatio;
0311         }
0312       }
0313     }
0314   }
0315 
0316   // No useful ratios, return zero amplitude and no time measurement
0317   if (0 == ratios_size)
0318     return;
0319 
0320   //  std::array < Tmax, C::MAXSAMPLES*(C::MAXSAMPLES-1)/2 > times_;
0321   Tmax timesAB_[C::MAXSAMPLES * (C::MAXSAMPLES - 1) / 2];
0322   unsigned int timesAB_size = 0;
0323 
0324   // make a vector of Tmax measurements that correspond to each ratio
0325   // and based on Alpha-Beta parameterization of the pulse shape
0326 
0327   for (unsigned int i = 0; i < ratios_size; i++) {
0328     double stepOverBeta = double(ratios_[i].step) / beta;
0329     double offset = double(ratios_[i].index) + alphabeta;
0330 
0331     double Rmin = ratios_[i].value - ratios_[i].error;
0332     if (Rmin < 0.001)
0333       Rmin = 0.001;
0334 
0335     double Rmax = ratios_[i].value + ratios_[i].error;
0336     double RLimit = Rlim[ratios_[i].step];
0337     if (Rmax > RLimit)
0338       Rmax = RLimit;
0339 
0340     double time1 =
0341         offset - ratios_[i].step / (myMath::fast_expf((stepOverBeta - myMath::fast_logf(Rmin)) / alpha) - 1.0);
0342     double time2 =
0343         offset - ratios_[i].step / (myMath::fast_expf((stepOverBeta - myMath::fast_logf(Rmax)) / alpha) - 1.0);
0344 
0345     // this is the time measurement based on the ratios[i]
0346     double tmax = 0.5 * (time1 + time2);
0347     double tmaxerr = 0.5 * std::sqrt((time1 - time2) * (time1 - time2));
0348 
0349     // calculate chi2
0350     sumAf = 0;
0351     sumff = 0;
0352     int itmin = std::max(-1, int(std::floor(tmax - alphabeta)));
0353     double loffset = (double(itmin) - tmax) * invalphabeta;
0354     for (unsigned int it = itmin + 1; it < amplitudes_.size(); it++) {
0355       loffset += invalphabeta;
0356       if (useless_[it])
0357         continue;
0358       double inverr2 = amplitudeIE2_[it];
0359       //      double offset = (double(it) - tmax)/alphabeta;
0360       double term1 = 1.0 + loffset;
0361       // assert(term1>1e-6);
0362       double f = (term1 > 1e-6) ? myMath::fast_expf(alpha * (myMath::fast_logf(term1) - loffset)) : 0;
0363       sumAf += amplitudes_[it] * (f * inverr2);
0364       sumff += f * (f * inverr2);
0365     }
0366 
0367     double chi2 = sumAA;
0368     double amp = 0;
0369     if (sumff > 0) {
0370       chi2 = sumAA - sumAf * (sumAf / sumff);
0371       amp = (sumAf / sumff);
0372     }
0373     chi2 /= sum0;
0374 
0375     // choose reasonable measurements. One might argue what is
0376     // reasonable and what is not.
0377     if (chi2 > 0 && tmaxerr > 0 && tmax > 0) {
0378       Tmax currentTmax = {ratios_[i].index, ratios_[i].step, tmax, tmaxerr, amp, chi2};
0379       timesAB_[timesAB_size++] = currentTmax;
0380     }
0381   }
0382 
0383   // no reasonable time measurements!
0384   if (0 == timesAB_size)
0385     return;
0386 
0387   // find minimum chi2
0388   double chi2min = 1.0e+9;
0389   //double timeMinimum = 5;
0390   //double errorMinimum = 999;
0391   for (unsigned int i = 0; i < timesAB_size; i++) {
0392     if (timesAB_[i].chi2 < chi2min) {
0393       chi2min = timesAB_[i].chi2;
0394       //timeMinimum = timesAB_[i].value;
0395       //errorMinimum = timesAB_[i].error;
0396     }
0397   }
0398 
0399   // make a weighted average of tmax measurements with "small" chi2
0400   // (within 1 sigma of statistical uncertainty :-)
0401   double chi2Limit = chi2min + 1.0;
0402   double time_max = 0;
0403   double time_wgt = 0;
0404   for (unsigned int i = 0; i < timesAB_size; i++) {
0405     if (timesAB_[i].chi2 < chi2Limit) {
0406       double inverseSigmaSquared = 1.0 / (timesAB_[i].error * timesAB_[i].error);
0407       time_wgt += inverseSigmaSquared;
0408       time_max += timesAB_[i].value * inverseSigmaSquared;
0409     }
0410   }
0411 
0412   tMaxAlphaBeta = time_max / time_wgt;
0413   tMaxErrorAlphaBeta = 1.0 / sqrt(time_wgt);
0414 
0415   // find amplitude and chi2
0416   sumAf = 0;
0417   sumff = 0;
0418   for (unsigned int i = 0; i < amplitudes_.size(); i++) {
0419     if (useless_[i])
0420       continue;
0421     double inverr2 = amplitudeIE2_[i];
0422     double offset = (double(i) - tMaxAlphaBeta) * invalphabeta;
0423     double term1 = 1.0 + offset;
0424     if (term1 > 1e-6) {
0425       double f = myMath::fast_expf(alpha * (myMath::fast_logf(term1) - offset));
0426       sumAf += amplitudes_[i] * (f * inverr2);
0427       sumff += f * (f * inverr2);
0428     }
0429   }
0430 
0431   if (sumff > 0) {
0432     ampMaxAlphaBeta = sumAf / sumff;
0433     double chi2AlphaBeta = (sumAA - sumAf * sumAf / sumff) / sum0;
0434     if (chi2AlphaBeta > NullChi2) {
0435       // null hypothesis is better
0436       return;
0437     }
0438 
0439   } else {
0440     // no visible pulse here
0441     return;
0442   }
0443 
0444   // if we got to this point, we have a reconstructied Tmax
0445   // using RatioAlphaBeta Method. To summarize:
0446   //
0447   //     tMaxAlphaBeta      - Tmax value
0448   //     tMaxErrorAlphaBeta - error on Tmax, but I would not trust it
0449   //     ampMaxAlphaBeta    - amplitude of the pulse
0450   //     ampMaxError_        - uncertainty of the time sample with max amplitude
0451   //
0452 
0453   // Do Ratio's Method with "large" pulses
0454   if (ampMaxAlphaBeta / ampMaxError_ > 5.0) {
0455     // make a vector of Tmax measurements that correspond to each
0456     // ratio. Each measurement have it's value and the error
0457 
0458     double time_max = 0;
0459     double time_wgt = 0;
0460 
0461     for (unsigned int i = 0; i < ratios_size; i++) {
0462       if (ratios_[i].step == 1 && ratios_[i].value >= timeFitLimits.first && ratios_[i].value <= timeFitLimits.second) {
0463         double time_max_i = ratios_[i].index;
0464 
0465         // calculate polynomial for Tmax
0466 
0467         double u = timeFitParameters[timeFitParameters.size() - 1];
0468         for (int k = timeFitParameters.size() - 2; k >= 0; k--) {
0469           u = u * ratios_[i].value + timeFitParameters[k];
0470         }
0471 
0472         // calculate derivative for Tmax error
0473         double du = (timeFitParameters.size() - 1) * timeFitParameters[timeFitParameters.size() - 1];
0474         for (int k = timeFitParameters.size() - 2; k >= 1; k--) {
0475           du = du * ratios_[i].value + k * timeFitParameters[k];
0476         }
0477 
0478         // running sums for weighted average
0479         double errorsquared = ratios_[i].error * ratios_[i].error * du * du;
0480         if (errorsquared > 0) {
0481           time_max += (time_max_i - u) / errorsquared;
0482           time_wgt += 1.0 / errorsquared;
0483           //            Tmax currentTmax =
0484           //  { ratios_[i].index, 1, (time_max_i - u),
0485           //sqrt(errorsquared),0,1 };
0486           // times_.push_back(currentTmax);
0487         }
0488       }
0489     }
0490 
0491     // calculate weighted average of all Tmax measurements
0492     if (time_wgt > 0) {
0493       tMaxRatio = time_max / time_wgt;
0494       tMaxErrorRatio = 1.0 / sqrt(time_wgt);
0495 
0496       // combine RatioAlphaBeta and Ratio Methods
0497 
0498       if (ampMaxAlphaBeta / ampMaxError_ > 10.0) {
0499         // use pure Ratio Method
0500         calculatedRechit_.timeMax = tMaxRatio;
0501         calculatedRechit_.timeError = tMaxErrorRatio;
0502 
0503       } else {
0504         // combine two methods
0505         calculatedRechit_.timeMax = (tMaxAlphaBeta * (10.0 - ampMaxAlphaBeta / ampMaxError_) +
0506                                      tMaxRatio * (ampMaxAlphaBeta / ampMaxError_ - 5.0)) /
0507                                     5.0;
0508         calculatedRechit_.timeError = (tMaxErrorAlphaBeta * (10.0 - ampMaxAlphaBeta / ampMaxError_) +
0509                                        tMaxErrorRatio * (ampMaxAlphaBeta / ampMaxError_ - 5.0)) /
0510                                       5.0;
0511       }
0512 
0513     } else {
0514       // use RatioAlphaBeta Method
0515       calculatedRechit_.timeMax = tMaxAlphaBeta;
0516       calculatedRechit_.timeError = tMaxErrorAlphaBeta;
0517     }
0518 
0519   } else {
0520     // use RatioAlphaBeta Method
0521     calculatedRechit_.timeMax = tMaxAlphaBeta;
0522     calculatedRechit_.timeError = tMaxErrorAlphaBeta;
0523   }
0524 }
0525 
0526 template <class C>
0527 void EcalUncalibRecHitRatioMethodAlgo<C>::computeAmplitude(std::vector<double> &amplitudeFitParameters) {
0528   calculatedRechit_.amplitudeMax = computeAmplitudeImpl(amplitudeFitParameters, 1., 1.);
0529 }
0530 
0531 template <class C>
0532 double EcalUncalibRecHitRatioMethodAlgo<C>::computeAmplitudeImpl(std::vector<double> &amplitudeFitParameters,
0533                                                                  double corr4,
0534                                                                  double corr6) {
0535   ////////////////////////////////////////////////////////////////
0536   //                                                            //
0537   //             CALCULATE AMPLITUDE                            //
0538   //                                                            //
0539   ////////////////////////////////////////////////////////////////
0540 
0541   double alpha = amplitudeFitParameters[0];
0542   double beta = amplitudeFitParameters[1];
0543 
0544   // calculate pedestal, again
0545 
0546   double pedestalLimit = calculatedRechit_.timeMax - (alpha * beta) - 1.0;
0547   double sumA = 0;
0548   double sumF = 0;
0549   double sumAF = 0;
0550   double sumFF = 0;
0551   double sum1 = 0;
0552   for (unsigned int i = 0; i < amplitudes_.size(); i++) {
0553     if (useless_[i])
0554       continue;
0555     double inverr2 = amplitudeIE2_[i];
0556     double f = 0;
0557     double termOne = 1 + (i - calculatedRechit_.timeMax) / (alpha * beta);
0558     if (termOne > 1.e-5)
0559       f = myMath::fast_expf(alpha * myMath::fast_logf(termOne) - (i - calculatedRechit_.timeMax) / beta);
0560 
0561     // apply range of interesting samples
0562 
0563     if ((i < pedestalLimit) || (f > 0.6 * corr6 && i <= calculatedRechit_.timeMax) ||
0564         (f > 0.4 * corr4 && i >= calculatedRechit_.timeMax)) {
0565       sum1 += inverr2;
0566       sumA += amplitudes_[i] * inverr2;
0567       sumF += f * inverr2;
0568       sumAF += (f * inverr2) * amplitudes_[i];
0569       sumFF += f * (f * inverr2);
0570     }
0571   }
0572 
0573   double amplitudeMax = 0;
0574   if (sum1 > 0) {
0575     double denom = sumFF * sum1 - sumF * sumF;
0576     if (std::abs(denom) > 1.0e-20) {
0577       amplitudeMax = (sumAF * sum1 - sumA * sumF) / denom;
0578     }
0579   }
0580   return amplitudeMax;
0581 }
0582 
0583 template <class C>
0584 EcalUncalibratedRecHit EcalUncalibRecHitRatioMethodAlgo<C>::makeRecHit(const C &dataFrame,
0585                                                                        const EcalSampleMask &sampleMask,
0586                                                                        const double *pedestals,
0587                                                                        const double *pedestalRMSes,
0588                                                                        const double *gainRatios,
0589                                                                        std::vector<double> &timeFitParameters,
0590                                                                        std::vector<double> &amplitudeFitParameters,
0591                                                                        std::pair<double, double> &timeFitLimits) {
0592   init(dataFrame, sampleMask, pedestals, pedestalRMSes, gainRatios);
0593   computeTime(timeFitParameters, timeFitLimits, amplitudeFitParameters);
0594   computeAmplitude(amplitudeFitParameters);
0595 
0596   // 1st parameters is id
0597   //
0598   // 2nd parameters is amplitude. It is calculated by this method.
0599   //
0600   // 3rd parameter is pedestal. It is not calculated. This method
0601   // relies on input parameters for pedestals and gain ratio. Return
0602   // zero.
0603   //
0604   // 4th parameter is jitter which is a bad choice to call Tmax. It is
0605   // calculated by this method (in 25 nsec clock units)
0606   //
0607   // GF subtract 5 so that jitter==0 corresponds to synchronous hit
0608   //
0609   //
0610   // 5th parameter is chi2. It is possible to calculate chi2 for
0611   // Tmax. It is possible to calculate chi2 for Amax. However, these
0612   // values are not very useful without TmaxErr and AmaxErr. This
0613   // method can return one value for chi2 but there are 4 different
0614   // parameters that have useful information about the quality of Amax
0615   // ans Tmax. For now we can return TmaxErr. The quality of Tmax and
0616   // Amax can be judged from the magnitude of TmaxErr
0617 
0618   return EcalUncalibratedRecHit(dataFrame.id(),
0619                                 calculatedRechit_.amplitudeMax,
0620                                 pedestal_,
0621                                 calculatedRechit_.timeMax - 5,
0622                                 calculatedRechit_.timeError);
0623 }
0624 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!