CMS 3D CMS Logo

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