EcalRecAlgos/interface/EcalUncalibRecHitTimeWeightsAlgo.h

0001 #ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitTimeWeightsAlgo_HH
0002 #define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitTimeWeightsAlgo_HH
0003
0004 /** \class EcalUncalibRecHitTimeWeightsAlgo
0005   *  Template used to compute amplitude, pedestal, time jitter, chi2 of a pulse
0006   *  using a weights method
0007   *
0008   *  \author J. Bendavid, E. Di Marco
0009   *
0010   *  The chi2 computation with matrix is replaced by the chi2express which is  moved outside the weight algo
0011   *  (need to clean up the interface in next iteration so that we do not pass-by useless arrays)
0012   *
0013   */
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "CondFormats/EcalObjects/interface/EcalPedestals.h"
0017 #include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
0018 #include "CondFormats/EcalObjects/interface/EcalWeightSet.h"
0019
0020 #include "RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h"
0021 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
0022
0023 #include "TVectorD.h"
0024 #include <vector>
0025
0026 template <class C>
0027 class EcalUncalibRecHitTimeWeightsAlgo {
0028 public:
0029   EcalUncalibRecHitTimeWeightsAlgo() {}
0030   virtual ~EcalUncalibRecHitTimeWeightsAlgo() {}
0031
0032   /// Compute time
0033   double time(const C &dataFrame,
0034               const std::vector<double> &amplitudes,
0035               const EcalPedestals::Item *aped,
0036               const EcalMGPAGainRatio *aGain,
0037               const FullSampleVector &fullpulse,
0038               const EcalWeightSet::EcalWeightMatrix **weights) {
0039     const unsigned int nsample = EcalDataFrame::MAXSAMPLES;
0040
0041     double maxamplitude = -std::numeric_limits<double>::max();
0042
0043     double pulsenorm = 0.;
0044     int iGainSwitch = 0;
0045
0046     ROOT::Math::SVector<double, nsample> pedSubSamples;
0047     for (unsigned int iSample = 0; iSample < nsample; iSample++) {
0048       const EcalMGPASample &sample = dataFrame.sample(iSample);
0049
0050       double amplitude = 0.;
0051       int gainId = sample.gainId();
0052
0053       double pedestal = 0.;
0054       double gainratio = 1.;
0055
0056       if (gainId == 0 || gainId == 3) {
0057         pedestal = aped->mean_x1;
0058         gainratio = aGain->gain6Over1() * aGain->gain12Over6();
0059         iGainSwitch = 1;
0060       } else if (gainId == 1) {
0061         pedestal = aped->mean_x12;
0062         gainratio = 1.;
0063         iGainSwitch = 0;
0064       } else if (gainId == 2) {
0065         pedestal = aped->mean_x6;
0066         gainratio = aGain->gain12Over6();
0067         iGainSwitch = 1;
0068       }
0069
0070       amplitude = ((double)(sample.adc()) - pedestal) * gainratio;
0071
0072       if (gainId == 0) {
0073         //saturation
0074         amplitude = (4095. - pedestal) * gainratio;
0075       }
0076
0077       pedSubSamples(iSample) = amplitude;
0078
0079       if (amplitude > maxamplitude) {
0080         maxamplitude = amplitude;
0081       }
0082       pulsenorm += fullpulse(iSample);
0083     }
0084
0085     std::vector<double>::const_iterator amplit;
0086     for (amplit = amplitudes.begin(); amplit < amplitudes.end(); ++amplit) {
0087       int ipulse = std::distance(amplitudes.begin(), amplit);
0088       int bx = ipulse - 5;
0089       int firstsamplet = std::max(0, bx + 3);
0090       int offset = 7 - 3 - bx;
0091
0092       TVectorD pulse;
0093       pulse.ResizeTo(nsample);
0094       for (unsigned int isample = firstsamplet; isample < nsample; ++isample) {
0095         pulse(isample) = fullpulse(isample + offset);
0096         pedSubSamples(isample) = std::max(0., pedSubSamples(isample) - amplitudes[ipulse] * pulse(isample) / pulsenorm);
0097       }
0098     }
0099
0100     // Compute parameters
0101     double amplitude_(-1.), jitter_(-1.);
0102     ROOT::Math::SVector<double, 3> param = (*(weights[iGainSwitch])) * pedSubSamples;
0103     amplitude_ = param(EcalUncalibRecHitRecAbsAlgo<C>::iAmplitude);
0104     if (amplitude_)
0105       jitter_ = -param(EcalUncalibRecHitRecAbsAlgo<C>::iTime) / amplitude_;
0106     else
0107       jitter_ = 0.;
0108
0109     return jitter_;
0110   }
0111 };
0112 #endif