EcalRecAlgos/interface/EcalUncalibRecHitRecChi2Algo.h

0001 #ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecChi2Algo_HH
0002 #define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecChi2Algo_HH
0003
0004 /** \class EcalUncalibRecHitRecChi2Algo
0005   *
0006   *  Template used to compute the chi2 of an MGPA pulse for in-time and out-of-time signals, algorithm based on the chi2express.
0007   *  The in-time chi2 is calculated against the time intercalibrations from the DB while the out-of-time chi2 is calculated
0008   *  against the Tmax measurement on event by event basis.
0009   *
0010   *  \author Konstantinos Theofilatos 02 Feb 2010
0011   */
0012
0013 #include "Math/SVector.h"
0014 #include "Math/SMatrix.h"
0015 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
0016 #include "CondFormats/EcalObjects/interface/EcalWeightSet.h"
0017 #include "SimCalorimetry/EcalSimAlgos/interface/EcalShapeBase.h"
0018 #include "CondFormats/EcalObjects/interface/EcalTimeCalibConstants.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020
0021 #include <vector>
0022
0023 template <class C>
0024 class EcalUncalibRecHitRecChi2Algo {
0025 public:
0026   // destructor
0027   virtual ~EcalUncalibRecHitRecChi2Algo() {}
0028
0029   EcalUncalibRecHitRecChi2Algo() {}
0030   EcalUncalibRecHitRecChi2Algo(const C& dataFrame,
0031                                const double amplitude,
0032                                const EcalTimeCalibConstant& timeIC,
0033                                const double amplitudeOutOfTime,
0034                                const double jitter,
0035                                const double* pedestals,
0036                                const double* pedestalsRMS,
0037                                const double* gainRatios,
0038                                const EcalShapeBase& testbeamPulseShape,
0039                                const std::vector<double>& chi2Parameters);
0040
0041   virtual double chi2() { return chi2_; }
0042   virtual double chi2OutOfTime() { return chi2OutOfTime_; }
0043
0044 private:
0045   double chi2_;
0046   double chi2OutOfTime_;
0047 };
0048
0049 template <class C>
0050 EcalUncalibRecHitRecChi2Algo<C>::EcalUncalibRecHitRecChi2Algo(const C& dataFrame,
0051                                                               const double amplitude,
0052                                                               const EcalTimeCalibConstant& timeIC,
0053                                                               const double amplitudeOutOfTime,
0054                                                               const double jitter,
0055                                                               const double* pedestals,
0056                                                               const double* pedestalsRMS,
0057                                                               const double* gainRatios,
0058                                                               const EcalShapeBase& testbeamPulseShape,
0059                                                               const std::vector<double>& chi2Parameters) {
0060   double noise_A = chi2Parameters[0];  // noise term for in-time chi2
0061   double const_A = chi2Parameters[1];  // constant term for in-time chi2
0062   double noise_B = chi2Parameters[2];  // noise term for out-of-time chi2
0063   double const_B = chi2Parameters[3];  // constant term for out-of-time chi2
0064
0065   chi2_ = 0;
0066   chi2OutOfTime_ = 0;
0067   double S_0 = 0;      // will store the first mgpa sample
0068   double ped_ave = 0;  // will store the average pedestal
0069
0070   int gainId0 = 1;
0071   int iGainSwitch = 0;
0072   bool isSaturated = false;
0073   for (int iSample = 0; iSample < C::MAXSAMPLES;
0074        iSample++)  // if gain switch use later the pedestal RMS, otherwise we use the pedestal from the DB
0075   {
0076     int gainId = dataFrame.sample(iSample).gainId();
0077     if (gainId == 0) {
0078       gainId = 3;  // if saturated, treat it as G1
0079       isSaturated = true;
0080     }
0081     if (gainId != gainId0)
0082       iGainSwitch = 1;
0083
0084     if (gainId == 1 && iSample == 0)
0085       S_0 = dataFrame.sample(iSample).adc();  // take only first presample to estimate the pedestal
0086     if (gainId == 1 && iSample < 3)
0087       ped_ave += (1 / 3.0) * dataFrame.sample(iSample).adc();  // take first 3 presamples to estimate the pedestal
0088   }
0089
0090   // compute testbeamPulseShape shape parameters
0091   double ADC_clock = 25;  // 25 ns
0092   double risingTime = testbeamPulseShape.timeToRise();
0093   double tzero = risingTime - 5 * ADC_clock;  // 5 samples before the peak
0094
0095   double shiftTime = +timeIC;                       // we put positive here
0096   double shiftTimeOutOfTime = -jitter * ADC_clock;  // we put negative here
0097
0098   bool readoutError = false;
0099
0100   for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
0101     int gainId = dataFrame.sample(iSample).gainId();
0102     if (dataFrame.sample(iSample).adc() == 0)
0103       readoutError = true;
0104     if (gainId == 0)
0105       continue;  // skip saturated samples
0106
0107     double ped =
0108         !iGainSwitch ? ped_ave : pedestals[gainId - 1];  // use dynamic pedestal for G12 and average pedestal for G6,G1
0109                                                          //double pedRMS = pedestalsRMS[gainId-1];
0110     double S_i = double(dataFrame.sample(iSample).adc());
0111
0112     // --- calculate in-time chi2
0113
0114     double f_i = (testbeamPulseShape)(tzero + shiftTime + iSample * ADC_clock);
0115     double R_i = (S_i - ped) * gainRatios[gainId - 1] - f_i * amplitude;
0116     double R_iErrorSquare = noise_A * noise_A + const_A * const_A * amplitude * amplitude;
0117
0118     chi2_ += R_i * R_i / R_iErrorSquare;
0119
0120     // --- calculate out-of-time chi2
0121
0122     double g_i = (testbeamPulseShape)(tzero + shiftTimeOutOfTime + iSample * ADC_clock);  // calculate out of time chi2
0123
0124     double R_iOutOfTime = (S_i - S_0) * gainRatios[gainId - 1] - g_i * amplitudeOutOfTime;
0125     double R_iOutOfTimeErrorSquare = noise_B * noise_B + const_B * const_B * amplitudeOutOfTime * amplitudeOutOfTime;
0126
0127     chi2OutOfTime_ += R_iOutOfTime * R_iOutOfTime / R_iOutOfTimeErrorSquare;
0128   }
0129
0130   if (!isSaturated && !iGainSwitch && chi2_ > 0 && chi2OutOfTime_ > 0) {
0131     chi2_ = 7 * (3 + log(chi2_));
0132     chi2_ = chi2_ < 0 ? 0 : chi2_;  // this is just a convinient mapping for storing in the calibRecHit bit map
0133     chi2OutOfTime_ = 7 * (3 + log(chi2OutOfTime_));
0134     chi2OutOfTime_ = chi2OutOfTime_ < 0 ? 0 : chi2OutOfTime_;
0135   } else {
0136     chi2_ = 0;
0137     chi2OutOfTime_ = 0;
0138   }
0139
0140   if (readoutError)  // rare situation
0141   {
0142     chi2_ = 99.0;  // chi2 is very large in these cases, put a code value to discriminate against normal noise
0143     chi2OutOfTime_ = 99.0;
0144   }
0145 }
0146
0147 #endif