CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions
EcalUncalibRecHitRecAnalFitAlgo< C > Class Template Reference

#include <EcalUncalibRecHitRecAnalFitAlgo.h>

Inheritance diagram for EcalUncalibRecHitRecAnalFitAlgo< C >:
EcalUncalibRecHitRecAbsAlgo< C >

Public Member Functions

 EcalUncalibRecHitRecAnalFitAlgo ()
 
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 24 of file EcalUncalibRecHitRecAnalFitAlgo.h.

Constructor & Destructor Documentation

◆ EcalUncalibRecHitRecAnalFitAlgo()

Definition at line 43 of file EcalUncalibRecHitRecAnalFitAlgo.h.

43  {
44  //In order to make fitting ROOT histograms thread safe
45  // one must call this undocumented function
46  TMinuitMinimizer::UseStaticMinuit(false);
47  }

◆ ~EcalUncalibRecHitRecAnalFitAlgo()

template<class C>
EcalUncalibRecHitRecAnalFitAlgo< C >::~EcalUncalibRecHitRecAnalFitAlgo ( )
inlineoverride

Definition at line 50 of file EcalUncalibRecHitRecAnalFitAlgo.h.

50 {};

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 53 of file EcalUncalibRecHitRecAnalFitAlgo.h.

Referenced by EcalUncalibRecHitWorkerAnalFit::run().

57  {
58  double amplitude_(-1.), pedestal_(-1.), jitter_(-1.), chi2_(-1.);
59 
60  // Get time samples
61  //HepMatrix frame(C::MAXSAMPLES, 1);
62  double frame[C::MAXSAMPLES];
63  // int gainId0 = dataFrame.sample(0).gainId();
64  int gainId0 = 1;
65  int iGainSwitch = 0;
66  double maxsample(-1);
67  int imax(-1);
68  bool isSaturated = false;
69  uint32_t flag = 0;
70  for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
71  int gainId = dataFrame.sample(iSample).gainId();
72  if (dataFrame.isSaturated()) {
73  gainId = 3;
74  isSaturated = true;
75  }
76 
77  if (gainId != gainId0)
78  ++iGainSwitch;
79  if (!iGainSwitch)
80  frame[iSample] = double(dataFrame.sample(iSample).adc());
81  else
82  frame[iSample] =
83  double(((double)(dataFrame.sample(iSample).adc()) - pedestals[gainId - 1]) * gainRatios[gainId - 1]);
84 
85  if (frame[iSample] > maxsample) {
86  maxsample = frame[iSample];
87  imax = iSample;
88  }
89  }
90 
91  // Compute parameters
92  //std::cout << "EcalUncalibRecHitRecAnalFitAlgo::makeRecHit() not yey implemented. returning dummy rechit" << std::endl;
93 
94  // prepare TGraph for analytic fit
95  double xarray[10] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9.};
96  TGraph graph(10, xarray, frame);
97 
98  // fit functions
99  TF1 pulseShape = TF1("pulseShape",
100  "[0]*pow((x - [3])/[1],[2])*exp(-[2]*(x - [1] - [3])/[1])",
101  imax - 1.,
102  imax + 3.,
103  TF1::EAddToList::kNo);
104  TF1 pedestal = TF1("pedestal", "[0]", 0., 2., TF1::EAddToList::kNo);
105 
106  //pulseShape parameters
107  // Amplitude
108  double FIT_A = (double)maxsample; //Amplitude
109  pulseShape.SetParameter(0, FIT_A);
110  pulseShape.SetParName(0, "Amplitude");
111  // T peak
112  double FIT_Tp = (double)imax; //T peak
113  pulseShape.SetParameter(1, FIT_Tp);
114  pulseShape.SetParName(1, "t_{P}");
115  // Alpha
116  double FIT_ALFA = 1.5; //Alpha
117  pulseShape.SetParameter(2, FIT_ALFA);
118  pulseShape.SetParName(2, "\\alpha");
119  // T off
120  double FIT_To = 3.; //T off
121  pulseShape.SetParameter(3, FIT_To);
122  pulseShape.SetParName(3, "t_{0}");
123 
124  // pedestal
125  pedestal.SetParameter(0, frame[0]);
126  pedestal.SetParName(0, "Pedestal");
127 
128  int result = graph.Fit(&pulseShape, "QRMN SERIAL");
129 
130  if (0 == result) {
131  double amplitude_value = pulseShape.GetParameter(0);
132 
133  graph.Fit(&pedestal, "QRLN SERIAL");
134  double pedestal_value = pedestal.GetParameter(0);
135 
136  if (!iGainSwitch)
137  amplitude_ = amplitude_value - pedestal_value;
138  else
139  amplitude_ = amplitude_value;
140 
141  pedestal_ = pedestal_value;
142  jitter_ = pulseShape.GetParameter(3);
143  chi2_ = 1.; // successful fit
144  if (isSaturated)
146  /*
147  std::cout << "separate fits\nA: " << amplitude_value << ", Ped: " << pedestal_value
148  << ", t0: " << jitter_ << ", tp: " << pulseShape.GetParameter(1)
149  << ", alpha: " << pulseShape.GetParameter(2)
150  << std::endl;
151  */
152  }
153 
154  return EcalUncalibratedRecHit(dataFrame.id(), amplitude_, pedestal_, jitter_ - 6, chi2_, flag);
155  }
constexpr int gainId(sample_type sample)
get the gainId (2 bits)
bool isSaturated(const Digi &digi, const int &maxADCvalue, int ifirst, int n)

◆ pedestalFunction()

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

Definition at line 37 of file EcalUncalibRecHitRecAnalFitAlgo.h.

37  {
38  double ped = par[0];
39  return ped;
40  };

◆ pulseShapeFunction()

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

Definition at line 26 of file EcalUncalibRecHitRecAnalFitAlgo.h.

26  {
27  double x = var[0];
28  double ampl = par[0];
29  double tp = par[1];
30  double alpha = par[2];
31  double t0 = par[3];
32 
33  double f = pow((x - t0) / tp, alpha) * exp(-alpha * (x - tp - t0) / tp);
34  return ampl * f;
35  };
float alpha
Definition: AMPTWrapper.h:105
double f[11][100]
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29