EcalRecAlgos/interface/EcalUncalibRecHitRecWeightsAlgo.h

0001 #ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecWeightsAlgo_HH
0002 #define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecWeightsAlgo_HH
0003
0004 /** \class EcalUncalibRecHitRecWeightsAlgo
0005   *  Template used to compute amplitude, pedestal, time jitter, chi2 of a pulse
0006   *  using a weights method
0007   *
0008   *  \author R. Bruneliere - A. Zabi
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 "Math/SVector.h"
0016 #include "Math/SMatrix.h"
0017 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
0018 #include "CondFormats/EcalObjects/interface/EcalWeightSet.h"
0019 #include "SimCalorimetry/EcalSimAlgos/interface/EcalShapeBase.h"
0020
0021 #include <vector>
0022
0023 template <class C>
0024 class EcalUncalibRecHitRecWeightsAlgo {
0025 public:
0026   // destructor
0027   virtual ~EcalUncalibRecHitRecWeightsAlgo(){};
0028
0029   /// Compute parameters
0030   virtual EcalUncalibratedRecHit makeRecHit(const C& dataFrame,
0031                                             const double* pedestals,
0032                                             const double* pedestalsRMS,
0033                                             const double* gainRatios,
0034                                             const EcalWeightSet::EcalWeightMatrix** weights,
0035                                             const EcalShapeBase& testbeamPulseShape) {
0036     double amplitude_(-1.), pedestal_(-1.), jitter_(-1.), chi2_(-1.);
0037     uint32_t flag = 0;
0038
0039     // Get time samples
0040     ROOT::Math::SVector<double, C::MAXSAMPLES> frame;
0041     int gainId0 = 1;
0042     int iGainSwitch = 0;
0043     bool isSaturated = false;
0044     for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
0045       int gainId = dataFrame.sample(iSample).gainId();
0046       //Handling saturation (treating saturated gainId as maximum gain)
0047       if (gainId == 0) {
0048         gainId = 3;
0049         //isSaturated = 1;
0050         // in pileup run May2012 samples 7,8,9,10 have gainid ==0
0051         // fix it like this: it won't hurt for the future SA20120512
0052         if (iSample == 4 || iSample == 5 || iSample == 6)
0053           isSaturated = true;
0054       }
0055
0056       //      if (gainId != gainId0) iGainSwitch = 1;
0057       // same problem as above: mark saturation only when physically
0058       // expected to occur SA20120513
0059       if ((gainId != gainId0) && (iSample == 4 || iSample == 5 || iSample == 6))
0060         iGainSwitch = 1;
0061       if (!iGainSwitch)
0062         frame(iSample) = double(dataFrame.sample(iSample).adc());
0063       else
0064         frame(iSample) =
0065             double(((double)(dataFrame.sample(iSample).adc()) - pedestals[gainId - 1]) * gainRatios[gainId - 1]);
0066     }
0067
0068     // Compute parameters
0069     ROOT::Math::SVector<double, 3> param = (*(weights[iGainSwitch])) * frame;
0070     amplitude_ = param(EcalUncalibRecHitRecAbsAlgo<C>::iAmplitude);
0071     pedestal_ = param(EcalUncalibRecHitRecAbsAlgo<C>::iPedestal);
0072     if (amplitude_)
0073       jitter_ = -param(EcalUncalibRecHitRecAbsAlgo<C>::iTime) / amplitude_;
0074     else
0075       jitter_ = 0.;
0076
0077     //When saturated gain flag i
0078     if (isSaturated) {
0079       flag = EcalUncalibratedRecHit::kSaturated;
0080       amplitude_ = double((4095. - pedestals[2]) * gainRatios[2]);
0081     }
0082     return EcalUncalibratedRecHit(dataFrame.id(), amplitude_, pedestal_, jitter_, chi2_, flag);
0083   }
0084
0085   ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
0086 };
0087 #endif