EcalUncalibRecHitRecAnalFitAlgo< C > Class Template Reference

#include <EcalUncalibRecHitRecAnalFitAlgo.h>

Inheritance diagram for EcalUncalibRecHitRecAnalFitAlgo< C >:

Public Member Functions
virtual EcalUncalibratedRecHit	makeRecHit (const C &dataFrame, const double *pedestals, const double *gainRatios, const EcalWeightSet::EcalWeightMatrix **weights, const EcalWeightSet::EcalChi2WeightMatrix **chi2Matrix)
	Compute parameters. More...
virtual	~EcalUncalibRecHitRecAnalFitAlgo ()

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

Id:

EcalUncalibRecHitRecAnalFitAlgo.h,v 1.11 2009/03/27 18:07:38 ferriff Exp

Date:

2009/03/27 18:07:38

Revision:

1.11

Author

A. Palma, Sh. Rahatlou Roma1

Definition at line 25 of file EcalUncalibRecHitRecAnalFitAlgo.h.

Constructor & Destructor Documentation

template<class C>

virtual EcalUncalibRecHitRecAnalFitAlgo< C >::~EcalUncalibRecHitRecAnalFitAlgo ( )

inlinevirtual

Definition at line 50 of file EcalUncalibRecHitRecAnalFitAlgo.h.

{ };

Member Function Documentation

template<class C>

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

inlinevirtual

Compute parameters.

Implements EcalUncalibRecHitRecAbsAlgo< C >.

Definition at line 54 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 = 0;
    uint32_t flag = 0;
    for(int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
      int gainId = dataFrame.sample(iSample).gainId(); 
      if ( dataFrame.isSaturated() != -1 ) 
    {
      gainId = 3;
      isSaturated = 1;
    }

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

template<class C>

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

inlineprivate

Definition at line 42 of file EcalUncalibRecHitRecAnalFitAlgo.h.

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

template<class C>

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

inlineprivate

Definition at line 31 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;
  };