Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:24:20

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<C>(){};
0028 
0029   EcalUncalibRecHitRecChi2Algo<C>(){};
0030   EcalUncalibRecHitRecChi2Algo<C>(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