EcalRecAlgos/interface/EcalUncalibRecHitRecAnalFitAlgo.h

0001 #ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecAnalFitAlgo_HH
0002 #define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecAnalFitAlgo_HH
0003
0004 /** \class EcalUncalibRecHitRecAnalFitAlgo
0005   *  Template used to compute amplitude, pedestal, time jitter, chi2 of a pulse
0006   *  using an analytical fit
0007   *
0008   *  \author A. Palma, Sh. Rahatlou Roma1
0009   */
0010
0011 #include "CLHEP/Matrix/Matrix.h"
0012 #include "CLHEP/Matrix/SymMatrix.h"
0013 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
0014 #include <vector>
0015 #include <string>
0016
0017 #include "TROOT.h"
0018 #include "TGraph.h"
0019 #include "TF1.h"
0020
0021 #include "TMinuitMinimizer.h"
0022
0023 template <class C>
0024 class EcalUncalibRecHitRecAnalFitAlgo : public EcalUncalibRecHitRecAbsAlgo<C> {
0025 private:
0026   double pulseShapeFunction(double* var, double* par) {
0027     double x = var[0];
0028     double ampl = par[0];
0029     double tp = par[1];
0030     double alpha = par[2];
0031     double t0 = par[3];
0032
0033     double f = pow((x - t0) / tp, alpha) * exp(-alpha * (x - tp - t0) / tp);
0034     return ampl * f;
0035   };
0036
0037   double pedestalFunction(double* var, double* par) {
0038     double ped = par[0];
0039     return ped;
0040   };
0041
0042 public:
0043   EcalUncalibRecHitRecAnalFitAlgo() {
0044     //In order to make fitting ROOT histograms thread safe
0045     // one must call this undocumented function
0046     TMinuitMinimizer::UseStaticMinuit(false);
0047   }
0048
0049   // destructor
0050   ~EcalUncalibRecHitRecAnalFitAlgo() override {}
0051
0052   /// Compute parameters
0053   EcalUncalibratedRecHit makeRecHit(const C& dataFrame,
0054                                     const double* pedestals,
0055                                     const double* gainRatios,
0056                                     const EcalWeightSet::EcalWeightMatrix** weights,
0057                                     const EcalWeightSet::EcalChi2WeightMatrix** chi2Matrix) override {
0058     double amplitude_(-1.), pedestal_(-1.), jitter_(-1.), chi2_(-1.);
0059
0060     // Get time samples
0061     //HepMatrix frame(C::MAXSAMPLES, 1);
0062     double frame[C::MAXSAMPLES];
0063     //    int gainId0 = dataFrame.sample(0).gainId();
0064     int gainId0 = 1;
0065     int iGainSwitch = 0;
0066     double maxsample(-1);
0067     int imax(-1);
0068     bool isSaturated = false;
0069     uint32_t flag = 0;
0070     for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
0071       int gainId = dataFrame.sample(iSample).gainId();
0072       if (dataFrame.isSaturated()) {
0073         gainId = 3;
0074         isSaturated = true;
0075       }
0076
0077       if (gainId != gainId0)
0078         ++iGainSwitch;
0079       if (!iGainSwitch)
0080         frame[iSample] = double(dataFrame.sample(iSample).adc());
0081       else
0082         frame[iSample] =
0083             double(((double)(dataFrame.sample(iSample).adc()) - pedestals[gainId - 1]) * gainRatios[gainId - 1]);
0084
0085       if (frame[iSample] > maxsample) {
0086         maxsample = frame[iSample];
0087         imax = iSample;
0088       }
0089     }
0090
0091     // Compute parameters
0092     //std::cout << "EcalUncalibRecHitRecAnalFitAlgo::makeRecHit() not yey implemented. returning dummy rechit" << std::endl;
0093
0094     // prepare TGraph for analytic fit
0095     double xarray[10] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9.};
0096     TGraph graph(10, xarray, frame);
0097
0098     // fit functions
0099     TF1 pulseShape = TF1("pulseShape",
0100                          "[0]*pow((x - [3])/[1],[2])*exp(-[2]*(x - [1] - [3])/[1])",
0101                          imax - 1.,
0102                          imax + 3.,
0103                          TF1::EAddToList::kNo);
0104     TF1 pedestal = TF1("pedestal", "[0]", 0., 2., TF1::EAddToList::kNo);
0105
0106     //pulseShape parameters
0107     // Amplitude
0108     double FIT_A = (double)maxsample;  //Amplitude
0109     pulseShape.SetParameter(0, FIT_A);
0110     pulseShape.SetParName(0, "Amplitude");
0111     // T peak
0112     double FIT_Tp = (double)imax;  //T peak
0113     pulseShape.SetParameter(1, FIT_Tp);
0114     pulseShape.SetParName(1, "t_{P}");
0115     // Alpha
0116     double FIT_ALFA = 1.5;  //Alpha
0117     pulseShape.SetParameter(2, FIT_ALFA);
0118     pulseShape.SetParName(2, "\\alpha");
0119     // T off
0120     double FIT_To = 3.;  //T off
0121     pulseShape.SetParameter(3, FIT_To);
0122     pulseShape.SetParName(3, "t_{0}");
0123
0124     // pedestal
0125     pedestal.SetParameter(0, frame[0]);
0126     pedestal.SetParName(0, "Pedestal");
0127
0128     int result = graph.Fit(&pulseShape, "QRMN SERIAL");
0129
0130     if (0 == result) {
0131       double amplitude_value = pulseShape.GetParameter(0);
0132
0133       graph.Fit(&pedestal, "QRLN SERIAL");
0134       double pedestal_value = pedestal.GetParameter(0);
0135
0136       if (!iGainSwitch)
0137         amplitude_ = amplitude_value - pedestal_value;
0138       else
0139         amplitude_ = amplitude_value;
0140
0141       pedestal_ = pedestal_value;
0142       jitter_ = pulseShape.GetParameter(3);
0143       chi2_ = 1.;  // successful fit
0144       if (isSaturated)
0145         flag = EcalUncalibratedRecHit::kSaturated;
0146       /*
0147       std::cout << "separate fits\nA: " <<  amplitude_value << ", Ped: " << pedestal_value
0148                 << ", t0: " << jitter_ << ", tp: " << pulseShape.GetParameter(1)
0149                 << ", alpha: " << pulseShape.GetParameter(2)
0150                 << std::endl;
0151       */
0152     }
0153
0154     return EcalUncalibratedRecHit(dataFrame.id(), amplitude_, pedestal_, jitter_ - 6, chi2_, flag);
0155   }
0156 };
0157 #endif