CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecChi2Algo.h

Go to the documentation of this file.
00001 #ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecChi2Algo_HH
00002 #define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecChi2Algo_HH
00003 
00013 #include "Math/SVector.h"
00014 #include "Math/SMatrix.h"
00015 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
00016 #include "CondFormats/EcalObjects/interface/EcalWeightSet.h"
00017 #include "SimCalorimetry/EcalSimAlgos/interface/EcalShapeBase.h"
00018 #include "CondFormats/EcalObjects/interface/EcalTimeCalibConstants.h"
00019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00020 
00021 
00022 #include <vector>
00023 
00024 template<class C> class EcalUncalibRecHitRecChi2Algo 
00025 {
00026   public:
00027 
00028   // destructor
00029   virtual ~EcalUncalibRecHitRecChi2Algo<C>() { };
00030 
00031   EcalUncalibRecHitRecChi2Algo<C>() { };
00032   EcalUncalibRecHitRecChi2Algo<C>(
00033                                   const C& dataFrame,
00034                                   const double amplitude,
00035                                   const EcalTimeCalibConstant& timeIC,
00036                                   const double amplitudeOutOfTime,
00037                                   const double jitter,
00038                                   const double* pedestals,
00039                                   const double* pedestalsRMS,
00040                                   const double* gainRatios,
00041                                   const EcalShapeBase & testbeamPulseShape,
00042                                   const std::vector<double>& chi2Parameters
00043   );
00044 
00045 
00046   virtual double chi2(){return chi2_;}
00047   virtual double chi2OutOfTime(){return chi2OutOfTime_;}
00048 
00049 
00050   private:
00051     double chi2_;
00052     double chi2OutOfTime_;
00053 };
00054 
00055 
00056 template <class C>
00057 EcalUncalibRecHitRecChi2Algo<C>::EcalUncalibRecHitRecChi2Algo(
00058                         const C& dataFrame, 
00059                         const double amplitude, 
00060                         const EcalTimeCalibConstant& timeIC, 
00061                         const double amplitudeOutOfTime, 
00062                         const double jitter, 
00063                         const double* pedestals, 
00064                         const double* pedestalsRMS,
00065                         const double* gainRatios,
00066                         const EcalShapeBase & testbeamPulseShape,
00067                         const std::vector<double>& chi2Parameters
00068 ) 
00069 {
00070 
00071     double noise_A = chi2Parameters[0]; // noise term for in-time chi2
00072     double const_A = chi2Parameters[1]; // constant term for in-time chi2
00073     double noise_B = chi2Parameters[2]; // noise term for out-of-time chi2
00074     double const_B = chi2Parameters[3]; // constant term for out-of-time chi2
00075 
00076 
00077     chi2_=0;
00078     chi2OutOfTime_=0;
00079     double S_0=0;  // will store the first mgpa sample
00080     double ped_ave=0;  // will store the average pedestal
00081 
00082     int gainId0 = 1;
00083     int iGainSwitch = 0;
00084     bool isSaturated = 0;
00085     for(int iSample = 0; iSample < C::MAXSAMPLES; iSample++) // if gain switch use later the pedestal RMS, otherwise we use the pedestal from the DB
00086     {
00087       int gainId = dataFrame.sample(iSample).gainId();
00088       if(gainId == 0)
00089       {
00090         gainId = 3; // if saturated, treat it as G1
00091         isSaturated = 1;
00092       }
00093       if(gainId != gainId0)iGainSwitch = 1;
00094 
00095       if(gainId==1 && iSample==0)S_0 = dataFrame.sample(iSample).adc(); // take only first presample to estimate the pedestal
00096       if(gainId==1 && iSample<3)ped_ave += (1/3.0)*dataFrame.sample(iSample).adc(); // take first 3 presamples to estimate the pedestal
00097     }
00098 
00099 
00100     // compute testbeamPulseShape shape parameters
00101     double ADC_clock = 25; // 25 ns
00102     double risingTime = testbeamPulseShape.timeToRise();
00103     double tzero = risingTime  - 5*ADC_clock;  // 5 samples before the peak
00104 
00105     double shiftTime = + timeIC; // we put positive here
00106     double shiftTimeOutOfTime = -jitter*ADC_clock; // we put negative here
00107     
00108 
00109     bool readoutError = false;
00110 
00111     for(int iSample = 0; iSample < C::MAXSAMPLES; iSample++)
00112     {
00113         int gainId = dataFrame.sample(iSample).gainId();
00114         if(dataFrame.sample(iSample).adc()==0)readoutError=true;
00115         if(gainId==0)continue; // skip saturated samples
00116 
00117 
00118         double ped = !iGainSwitch ? ped_ave:pedestals[gainId-1];  // use dynamic pedestal for G12 and average pedestal for G6,G1
00119         //double pedRMS = pedestalsRMS[gainId-1];
00120         double S_i = double(dataFrame.sample(iSample).adc());
00121 
00122         // --- calculate in-time chi2
00123 
00124 
00125         double f_i = (testbeamPulseShape)(tzero + shiftTime + iSample*ADC_clock);
00126         double R_i = (S_i- ped)*gainRatios[gainId-1] - f_i*amplitude;
00127         double R_iErrorSquare = noise_A*noise_A + const_A*const_A*amplitude*amplitude;
00128 
00129         chi2_ += R_i*R_i/R_iErrorSquare;
00130 
00131         // --- calculate out-of-time chi2
00132 
00133         double g_i = (testbeamPulseShape)(tzero + shiftTimeOutOfTime + iSample*ADC_clock); // calculate out of time chi2
00134 
00135         double R_iOutOfTime = (S_i- S_0)*gainRatios[gainId-1] - g_i*amplitudeOutOfTime;
00136         double R_iOutOfTimeErrorSquare = noise_B*noise_B + const_B*const_B*amplitudeOutOfTime*amplitudeOutOfTime;
00137 
00138  
00139         chi2OutOfTime_ += R_iOutOfTime*R_iOutOfTime/R_iOutOfTimeErrorSquare;
00140     }
00141 
00142 
00143     if(!isSaturated && !iGainSwitch && chi2_>0 && chi2OutOfTime_>0)     
00144     {   
00145         chi2_ = 7*(3+log(chi2_)); chi2_ = chi2_<0 ? 0:chi2_; // this is just a convinient mapping for storing in the calibRecHit bit map 
00146         chi2OutOfTime_ = 7*(3+log(chi2OutOfTime_)); chi2OutOfTime_ = chi2OutOfTime_<0 ? 0:chi2OutOfTime_; 
00147     }else
00148     {
00149         chi2_=0;
00150         chi2OutOfTime_=0;
00151     }
00152 
00153     if(readoutError) // rare situation
00154     {
00155         chi2_=99.0;  // chi2 is very large in these cases, put a code value to discriminate against normal noise
00156         chi2OutOfTime_=99.0;
00157     }
00158 }
00159 
00160 #endif