CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
EcalUncalibRecHitRecAnalFitAlgo.h
Go to the documentation of this file.
1 #ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecAnalFitAlgo_HH
2 #define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecAnalFitAlgo_HH
3 
11 #include "CLHEP/Matrix/Matrix.h"
12 #include "CLHEP/Matrix/SymMatrix.h"
14 #include <vector>
15 #include <string>
16 
17 #include "TROOT.h"
18 #include "TGraph.h"
19 #include "TF1.h"
20 
21 #include "TMinuitMinimizer.h"
22 
23 template <class C>
25 private:
26  double pulseShapeFunction(double* var, double* par) {
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  };
36 
37  double pedestalFunction(double* var, double* par) {
38  double ped = par[0];
39  return ped;
40  };
41 
42 public:
44  //In order to make fitting ROOT histograms thread safe
45  // one must call this undocumented function
46  TMinuitMinimizer::UseStaticMinuit(false);
47  }
48 
49  // destructor
51 
54  const double* pedestals,
55  const double* gainRatios,
57  const EcalWeightSet::EcalChi2WeightMatrix** chi2Matrix) override {
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  }
156 };
157 #endif
double pedestalFunction(double *var, double *par)
math::Matrix< 10, 10 >::type EcalChi2WeightMatrix
Definition: EcalWeightSet.h:20
math::Matrix< 3, 10 >::type EcalWeightMatrix
Definition: EcalWeightSet.h:19
double f[11][100]
static constexpr int MAXSAMPLES
constexpr int gainId(sample_type sample)
get the gainId (2 bits)
EcalUncalibratedRecHit makeRecHit(const C &dataFrame, const double *pedestals, const double *gainRatios, const EcalWeightSet::EcalWeightMatrix **weights, const EcalWeightSet::EcalChi2WeightMatrix **chi2Matrix) override
Compute parameters.
bool isSaturated(const Digi &digi, const int &maxADCvalue, int ifirst, int n)
double pulseShapeFunction(double *var, double *par)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29