EcalUncalibRecHitRecAnalFitAlgo.h

Go to the documentation of this file.

#ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecAnalFitAlgo_HH
#define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecAnalFitAlgo_HH

#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
#include <vector>
#include <string>

#include "TROOT.h"
#include "TMinuit.h"
#include "TGraph.h"
#include "TF1.h"

template <class C>
class EcalUncalibRecHitRecAnalFitAlgo : public EcalUncalibRecHitRecAbsAlgo<C> {
private:
  double pulseShapeFunction(double* var, double* par) {
    double x = var[0];
    double ampl = par[0];
    double tp = par[1];
    double alpha = par[2];
    double t0 = par[3];

    double f = pow((x - t0) / tp, alpha) * exp(-alpha * (x - tp - t0) / tp);
    return ampl * f;
  };

  double pedestalFunction(double* var, double* par) {
    double ped = par[0];
    return ped;
  };

public:
  // destructor
  ~EcalUncalibRecHitRecAnalFitAlgo<C>() override{};

  EcalUncalibratedRecHit makeRecHit(const C& dataFrame,
                                    const double* pedestals,
                                    const double* gainRatios,
                                    const EcalWeightSet::EcalWeightMatrix** weights,
                                    const EcalWeightSet::EcalChi2WeightMatrix** chi2Matrix) override {
    double amplitude_(-1.), pedestal_(-1.), jitter_(-1.), chi2_(-1.);

    // Get time samples
    //HepMatrix frame(C::MAXSAMPLES, 1);
    double frame[C::MAXSAMPLES];
    //    int gainId0 = dataFrame.sample(0).gainId();
    int gainId0 = 1;
    int iGainSwitch = 0;
    double maxsample(-1);
    int imax(-1);
    bool isSaturated = false;
    uint32_t flag = 0;
    for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
      int gainId = dataFrame.sample(iSample).gainId();
      if (dataFrame.isSaturated()) {
        gainId = 3;
        isSaturated = true;
      }

      if (gainId != gainId0)
        ++iGainSwitch;
      if (!iGainSwitch)
        frame[iSample] = double(dataFrame.sample(iSample).adc());
      else
        frame[iSample] =
            double(((double)(dataFrame.sample(iSample).adc()) - pedestals[gainId - 1]) * gainRatios[gainId - 1]);

      if (frame[iSample] > maxsample) {
        maxsample = frame[iSample];
        imax = iSample;
      }
    }

    // Compute parameters
    //std::cout << "EcalUncalibRecHitRecAnalFitAlgo::makeRecHit() not yey implemented. returning dummy rechit" << std::endl;

    // prepare TGraph for analytic fit
    double xarray[10] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9.};
    TGraph graph(10, xarray, frame);

    // fit functions
    TF1 pulseShape =
        TF1("pulseShape", "[0]*pow((x - [3])/[1],[2])*exp(-[2]*(x - [1] - [3])/[1])", imax - 1., imax + 3.);
    TF1 pedestal = TF1("pedestal", "[0]", 0., 2.);

    //TF1 pulseShape = TF1("pulseShape",pulseShapeFunction,imax-1.,imax+3.);
    //TF1 pedestal = TF1("pedestal",pedestalFunction,0.,2.);
    TF1 pluseAndPed = TF1("pulseAndPed", "pedestal+pulseShape");

    //pulseShape parameters
    // Amplitude
    double FIT_A = (double)maxsample;  //Amplitude
    pulseShape.SetParameter(0, FIT_A);
    pulseShape.SetParName(0, "Amplitude");
    // T peak
    double FIT_Tp = (double)imax;  //T peak
    pulseShape.SetParameter(1, FIT_Tp);
    pulseShape.SetParName(1, "t_{P}");
    // Alpha
    double FIT_ALFA = 1.5;  //Alpha
    pulseShape.SetParameter(2, FIT_ALFA);
    pulseShape.SetParName(2, "\\alpha");
    // T off
    double FIT_To = 3.;  //T off
    pulseShape.SetParameter(3, FIT_To);
    pulseShape.SetParName(3, "t_{0}");

    // pedestal
    pedestal.SetParameter(0, frame[0]);
    pedestal.SetParName(0, "Pedestal");

    graph.Fit(&pulseShape, "QRM");
    //TF1 *pulseShape2=graph.GetFunction("pulseShape");

    if (std::string(gMinuit->fCstatu.Data()) == std::string("CONVERGED ")) {
      double amplitude_value = pulseShape.GetParameter(0);

      graph.Fit(&pedestal, "QRL");
      //TF1 *pedestal2=graph.GetFunction("pedestal");
      double pedestal_value = pedestal.GetParameter(0);

      if (!iGainSwitch)
        amplitude_ = amplitude_value - pedestal_value;
      else
        amplitude_ = amplitude_value;

      pedestal_ = pedestal_value;
      jitter_ = pulseShape.GetParameter(3);
      chi2_ = 1.;  // successful fit
      if (isSaturated)
        flag = EcalUncalibratedRecHit::kSaturated;
      /*
      std::cout << "separate fits\nA: " <<  amplitude_value << ", Ped: " << pedestal_value
                << ", t0: " << jitter_ << ", tp: " << pulseShape.GetParameter(1)
                << ", alpha: " << pulseShape.GetParameter(2)
                << std::endl;
      */
    }

    return EcalUncalibratedRecHit(dataFrame.id(), amplitude_, pedestal_, jitter_ - 6, chi2_, flag);
  }
};
#endif