EcalUncalibRecHitRecAnalFitAlgo< C > Class Template Reference

#include <EcalUncalibRecHitRecAnalFitAlgo.h>

Inheritance diagram for EcalUncalibRecHitRecAnalFitAlgo< C >:

Public Member Functions
EcalUncalibratedRecHit	makeRecHit (const C &dataFrame, const double *pedestals, const double *gainRatios, const EcalWeightSet::EcalWeightMatrix **weights, const EcalWeightSet::EcalChi2WeightMatrix **chi2Matrix) override
	Compute parameters. More...
	~EcalUncalibRecHitRecAnalFitAlgo () override
Public Member Functions inherited from EcalUncalibRecHitRecAbsAlgo< C >
virtual	~EcalUncalibRecHitRecAbsAlgo ()=default
	Constructor. More...

Private Member Functions
double	pedestalFunction (double *var, double *par)
double	pulseShapeFunction (double *var, double *par)

Additional Inherited Members
Public Types inherited from EcalUncalibRecHitRecAbsAlgo< C >
enum	{ nWeightsRows = 3, iAmplitude = 0, iPedestal = 1, iTime = 2 }

Detailed Description

template<class C>
class EcalUncalibRecHitRecAnalFitAlgo< C >

Template used to compute amplitude, pedestal, time jitter, chi2 of a pulse using an analytical fit

Author

A. Palma, Sh. Rahatlou Roma1

Definition at line 23 of file EcalUncalibRecHitRecAnalFitAlgo.h.

Constructor & Destructor Documentation

~EcalUncalibRecHitRecAnalFitAlgo()

template<class C>

EcalUncalibRecHitRecAnalFitAlgo< C >::~EcalUncalibRecHitRecAnalFitAlgo ( )

inlineoverride

Definition at line 43 of file EcalUncalibRecHitRecAnalFitAlgo.h.

{};

Member Function Documentation

makeRecHit()

template<class C>

EcalUncalibratedRecHit EcalUncalibRecHitRecAnalFitAlgo< C >::makeRecHit ( const C &  dataFrame, const double *  pedestals, const double *  gainRatios, const EcalWeightSet::EcalWeightMatrix **  weights, const EcalWeightSet::EcalChi2WeightMatrix **  chi2Matrix  )

inlineoverridevirtual

Compute parameters.

Implements EcalUncalibRecHitRecAbsAlgo< C >.

Definition at line 46 of file EcalUncalibRecHitRecAnalFitAlgo.h.

Referenced by EcalUncalibRecHitWorkerAnalFit::run().

                                                                                                   {
    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);
  }

pedestalFunction()

template<class C>

double EcalUncalibRecHitRecAnalFitAlgo< C >::pedestalFunction ( double *  var, double *  par  )

inlineprivate

Definition at line 36 of file EcalUncalibRecHitRecAnalFitAlgo.h.

                                                    {
    double ped = par[0];
    return ped;
  };

pulseShapeFunction()

template<class C>

double EcalUncalibRecHitRecAnalFitAlgo< C >::pulseShapeFunction ( double *  var, double *  par  )

inlineprivate

Definition at line 25 of file EcalUncalibRecHitRecAnalFitAlgo.h.

                                                      {
    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;
  };