CMS 3D CMS Logo

ESRecHitAnalyticAlgo.cc
Go to the documentation of this file.
3 
4 #include <iostream>
5 #include <cmath>
6 
8 
10 
11 double* ESRecHitAnalyticAlgo::EvalAmplitude(const ESDataFrame& digi, double ped) const {
12  double* fitresults = new double[3]{};
13  double adc[3]{};
14 
15  for (int i = 0; i < digi.size(); i++)
16  adc[i] = digi.sample(i).adc() - ped;
17 
18  double status = 0;
19  if (adc[0] > 20)
20  status = 14;
21  if (adc[1] < 0 || adc[2] < 0)
22  status = 10;
23  if (adc[0] > adc[1] && adc[0] > adc[2])
24  status = 8;
25  if (adc[2] > adc[1] && adc[2] > adc[0])
26  status = 9;
27  double r12 = (adc[1] != 0) ? adc[0] / adc[1] : 99999;
28  double r23 = (adc[2] != 0) ? adc[1] / adc[2] : 99999;
29  if (r12 > ratioCuts_->getR12High())
30  status = 5;
31  if (r23 > ratioCuts_->getR23High())
32  status = 6;
33  if (r23 < ratioCuts_->getR23Low())
34  status = 7;
35 
36  if (int(status) == 0) {
37  double A1 = adc[1];
38  double A2 = adc[2];
39 
40  // t0 from analytical formula:
41  double n = 1.798;
42  double w = 0.07291;
43  double DeltaT = 25.;
44  double aaa = log(A2 / A1) / n;
45  double bbb = w / n * DeltaT;
46  double ccc = exp(aaa + bbb);
47 
48  //double t0 = (2.-ccc)/(ccc-1) * DeltaT + 5;
49  double t0 = (2. - ccc) / (1. - ccc) * DeltaT - 5;
50 
51  // A from analytical formula:
52  double t1 = 20.;
53  //double t2 = 45.;
54  double A_1 = pow(w / n * (t1), n) * exp(n - w * (t1));
55  //double A_2 = pow(w/n*(t2),n) * exp(n-w*(t2));
56  double AA1 = A1 / A_1;
57  //double AA2 = A2 / A_2;
58 
59  fitresults[0] = AA1;
60  fitresults[1] = t0;
61 
62  if (adc[1] > 2800 && adc[2] > 2800)
63  status = 11;
64  else if (adc[1] > 2800)
65  status = 12;
66  else if (adc[2] > 2800)
67  status = 13;
68 
69  } else {
70  fitresults[0] = 0;
71  fitresults[1] = -999;
72  }
73 
74  fitresults[2] = status;
75 
76  return fitresults;
77 }
78 
80  ESPedestals::const_iterator it_ped = peds_->find(digi.id());
81 
84 
86 
87  double* results;
88 
89  results = EvalAmplitude(digi, it_ped->getMean());
90 
91  double energy = results[0];
92  double t0 = results[1];
93  int status = (int)results[2];
94  delete[] results;
95 
96  double mipCalib = (fabs(cos(*it_ang)) != 0.) ? (*it_mip) / fabs(cos(*it_ang)) : 0.;
97  energy *= (mipCalib != 0.) ? MIPGeV_ / mipCalib : 0.;
98 
99  LogDebug("ESRecHitAnalyticAlgo") << "ESRecHitAnalyticAlgo : reconstructed energy " << energy << " T0 " << t0;
100 
101  EcalRecHit rechit(digi.id(), energy, t0);
102 
103  if (it_status->getStatusCode() == 1) {
105  } else {
106  if (status == 0)
107  rechit.setFlag(EcalRecHit::kESGood);
108  else if (status == 5)
109  rechit.setFlag(EcalRecHit::kESBadRatioFor12);
110  else if (status == 6)
111  rechit.setFlag(EcalRecHit::kESBadRatioFor23Upper);
112  else if (status == 7)
113  rechit.setFlag(EcalRecHit::kESBadRatioFor23Lower);
114  else if (status == 8)
115  rechit.setFlag(EcalRecHit::kESTS1Largest);
116  else if (status == 9)
117  rechit.setFlag(EcalRecHit::kESTS3Largest);
118  else if (status == 10)
119  rechit.setFlag(EcalRecHit::kESTS3Negative);
120  else if (status == 11)
121  rechit.setFlag(EcalRecHit::kESSaturated);
122  else if (status == 12)
123  rechit.setFlag(EcalRecHit::kESTS2Saturated);
124  else if (status == 13)
125  rechit.setFlag(EcalRecHit::kESTS3Saturated);
126  else if (status == 14)
127  rechit.setFlag(EcalRecHit::kESTS13Sigmas);
128  }
129 
130  return rechit;
131 }
const ESDetId & id() const
Definition: ESDataFrame.h:19
float getR12High() const
EcalRecHit reconstruct(const ESDataFrame &digi) const
int size() const
Definition: ESDataFrame.h:21
T w() const
const ESIntercalibConstants * mips_
void setFlag(int flag)
set the flags (from Flags or ESFlags)
Definition: EcalRecHit.h:183
const_iterator find(uint32_t rawId) const
const ESRecHitRatioCuts * ratioCuts_
float getR23High() const
int adc() const
get the ADC sample (singed 16 bits)
Definition: ESSample.h:16
const ESSample & sample(int i) const
Definition: ESDataFrame.h:24
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
const ESPedestals * peds_
const ESChannelStatus * channelStatus_
const ESAngleCorrectionFactors * ang_
std::vector< Item >::const_iterator const_iterator
double * EvalAmplitude(const ESDataFrame &digi, double ped) const
const self & getMap() const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
uint16_t *__restrict__ uint16_t const *__restrict__ adc
#define LogDebug(id)