CMS 3D CMS Logo
ESRecHitSimAlgo.cc
Go to the documentation of this file.
1 #include "RecoLocalCalo/EcalRecAlgos/interface/ESRecHitSimAlgo.h"
2 #include "FWCore/MessageLogger/interface/MessageLogger.h"
3 
4 #include <cmath>
5 #include <vdt/vdtMath.h>
6 
7 #include <cstdlib>
8 #include <cstring>
9 #include <cassert>
10 #include <iostream>
11 
12 EcalRecHit::ESFlags ESRecHitSimAlgo::evalAmplitude(float* results, const ESDataFrame& digi, float ped) const {
13  float energy = 0;
14  float adc[3]{};
15  float pw[3]{};
16  pw[0] = w0_;
17  pw[1] = w1_;
18  pw[2] = w2_;
19 
20  for (int i = 0; i < digi.size(); i++) {
21  energy += pw[i] * (digi.sample(i).adc() - ped);
22  LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : Digi " << i << " ADC counts " << digi.sample(i).adc() << " Ped "
23  << ped;
24  //std::cout<<i<<" "<<digi.sample(i).adc()<<" "<<ped<<" "<<pw[i]<<std::endl;
25  adc[i] = digi.sample(i).adc() - ped;
26  }
27 
28  EcalRecHit::ESFlags status = EcalRecHit::kESGood;
29  if (adc[0] > 20.f)
30  status = EcalRecHit::kESTS13Sigmas; // 14;
31  if (adc[1] <= 0 || adc[2] <= 0)
32  status = EcalRecHit::kESTS3Negative; // 10;
33  if (adc[0] > adc[1] && adc[0] > adc[2])
34  status = EcalRecHit::kESTS1Largest; // 8;
35  if (adc[2] > adc[1] && adc[2] > adc[0])
36  status = EcalRecHit::kESTS3Largest; // 9;
37  auto r12 = (adc[1] != 0) ? adc[0] / adc[1] : 99.f;
38  auto r23 = (adc[2] != 0) ? adc[1] / adc[2] : 99.f;
39  if (r12 > ratioCuts_->getR12High())
40  status = EcalRecHit::kESBadRatioFor12; // 5;
41  if (r23 > ratioCuts_->getR23High())
42  status = EcalRecHit::kESBadRatioFor23Upper; // 6;
43  if (r23 < ratioCuts_->getR23Low())
44  status = EcalRecHit::kESBadRatioFor23Lower; // 7
45 
46  auto A1 = adc[1];
47  auto A2 = adc[2];
48 
49  // t0 from analytical formula:
50  constexpr float n = 1.798;
51  constexpr float w = 0.07291;
52  constexpr float DeltaT = 25.;
53  auto aaa = (A2 > 0 && A1 > 0) ? std::log(A2 / A1) / n : 20.f; // if A1=0, t0=20
54  constexpr float bbb = w / n * DeltaT;
55  auto ccc = std::exp(aaa + bbb);
56 
57  auto t0 = (2.f - ccc) / (1.f - ccc) * DeltaT - 5.f;
58 
59  // A from analytical formula:
60  constexpr float t1 = 20.;
61 #if defined(__clang__) || defined(__INTEL_COMPILER)
62  const float A_1 = 1. / (std::pow(w / n * (t1), n) * std::exp(n - w * (t1)));
63 #else
64  constexpr float A_1 = 1. / (std::pow(w / n * (t1), n) * std::exp(n - w * (t1)));
65 #endif
66  auto AA1 = A1 * A_1;
67 
68  if (adc[1] > 2800.f && adc[2] > 2800.f)
69  status = EcalRecHit::kESSaturated;
70  else if (adc[1] > 2800.f)
71  status = EcalRecHit::kESTS2Saturated;
72  else if (adc[2] > 2800.f)
73  status = EcalRecHit::kESTS3Saturated;
74 
75  results[0] = energy; // energy with weight method
76  results[1] = t0; // timing
77  results[2] = AA1; // energy with analytic method
78 
79  return status;
80 }
81 
82 EcalRecHit ESRecHitSimAlgo::reconstruct(const ESDataFrame& digi) const {
83  auto ind = digi.id().hashedIndex();
84 
85  auto const& ped = peds_->preshower(ind);
86  auto const& mip = mips_->getMap().preshower(ind);
87  auto const& ang = ang_->getMap().preshower(ind);
88  auto const& statusCh = channelStatus_->getMap().preshower(ind);
89 
90  float results[3]{};
91 
92  auto status = evalAmplitude(results, digi, ped.getMean());
93 
94  auto energy = results[0];
95  auto t0 = results[1];
96  auto otenergy = results[2] * 1000000.f; // set out-of-time energy to keV
97 
98  auto mipCalib = (mip != 0.f) ? MIPGeV_ * std::abs(vdt::fast_cosf(ang)) / (mip) : 0.f;
99  energy *= mipCalib;
100  otenergy *= mipCalib;
101 
102  LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : reconstructed energy " << energy;
103 
104  EcalRecHit rechit(digi.id(), energy, t0);
105  // edm: this is just a placeholder for alternative energy reconstruction,
106  // so put it in the same float, with different name
107  // rechit.setOutOfTimeEnergy(otenergy);
108  rechit.setEnergyError(otenergy);
109 
110  rechit.setFlag(statusCh.getStatusCode() == 1 ? EcalRecHit::kESDead : status);
111 
112  return rechit;
113 }
114 
115 /*
116 
117  auto oldHit = oldreconstruct(digi);
118 
119  assert(rechit.recoFlag()==oldHit.recoFlag());
120 
121  if (oldHit.energy()>0)
122  std::cout << "ESd " << digi.id() <<" : "
123  << rechit.energy()<<"," << oldHit.energy() << " "
124  << rechit.time()<<"," << oldHit.time() << " "
125  << rechit.outOfTimeEnergy()<<"," << oldHit.outOfTimeEnergy() << std::endl;
126 
127 
128  auto bitdiff = [](float a, float b)->int {
129  int ia, ib;
130  memcpy(&ia,&a,4);
131  memcpy(&ib,&b,4);
132  return std::abs(ia-ib);
133  };
134 
135  auto d0 = bitdiff( rechit.energy(), rechit.energy() );
136  auto d1 = bitdiff( rechit.time(), rechit.time() );
137  auto d2 = bitdiff(rechit.outOfTimeEnergy(), oldHit.outOfTimeEnergy());
138 
139  struct aMax {
140  ~aMax() { std::cout << "\n\nmax deviation " << m << " " << mp << std::endl; }
141  int m=0;
142  int mp=0;
143  };
144 
145  static aMax am;
146  am.m = std::max(am.m,std::max(d0,std::max(d1,d2)));
147  if (oldHit.energy()>0) am.mp = std::max(am.m,std::max(d0,std::max(d1,d2)));
148 
149 
150  return rechit;
151 }
152 */
153 
154 //
155 
156 double* ESRecHitSimAlgo::oldEvalAmplitude(
157  const ESDataFrame& digi, const double& ped, const double& w0, const double& w1, const double& w2) const {
158  double* results = new double[4]{};
159  float energy = 0;
160  double adc[3]{};
161  float pw[3]{};
162  pw[0] = w0;
163  pw[1] = w1;
164  pw[2] = w2;
165 
166  for (int i = 0; i < digi.size(); i++) {
167  energy += pw[i] * (digi.sample(i).adc() - ped);
168  LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : Digi " << i << " ADC counts " << digi.sample(i).adc() << " Ped "
169  << ped;
170  //std::cout<<i<<" "<<digi.sample(i).adc()<<" "<<ped<<" "<<pw[i]<<std::endl;
171  adc[i] = digi.sample(i).adc() - ped;
172  }
173 
174  double status = 0;
175  if (adc[0] > 20)
176  status = 14;
177  if (adc[1] <= 0 || adc[2] <= 0)
178  status = 10;
179  if (adc[0] > adc[1] && adc[0] > adc[2])
180  status = 8;
181  if (adc[2] > adc[1] && adc[2] > adc[0])
182  status = 9;
183  double r12 = (adc[1] != 0) ? adc[0] / adc[1] : 99;
184  double r23 = (adc[2] != 0) ? adc[1] / adc[2] : 99;
185  if (r12 > ratioCuts_->getR12High())
186  status = 5;
187  if (r23 > ratioCuts_->getR23High())
188  status = 6;
189  if (r23 < ratioCuts_->getR23Low())
190  status = 7;
191 
192  double A1 = adc[1];
193  double A2 = adc[2];
194 
195  // t0 from analytical formula:
196  double n = 1.798;
197  double w = 0.07291;
198  double DeltaT = 25.;
199  double aaa = (A2 > 0 && A1 > 0) ? log(A2 / A1) / n : 20.; // if A1=0, t0=20
200  double bbb = w / n * DeltaT;
201  double ccc = exp(aaa + bbb);
202 
203  double t0 = (2. - ccc) / (1. - ccc) * DeltaT - 5;
204 
205  // A from analytical formula:
206  double t1 = 20.;
207  double A_1 = pow(w / n * (t1), n) * exp(n - w * (t1));
208  double AA1 = (A_1 != 0.) ? A1 / A_1 : 0.;
209 
210  if (adc[1] > 2800 && adc[2] > 2800)
211  status = 11;
212  else if (adc[1] > 2800)
213  status = 12;
214  else if (adc[2] > 2800)
215  status = 13;
216 
217  results[0] = energy; // energy with weight method
218  results[1] = t0; // timing
219  results[2] = status; // hit status
220  results[3] = AA1; // energy with analytic method
221 
222  return results;
223 }
224 
225 EcalRecHit ESRecHitSimAlgo::oldreconstruct(const ESDataFrame& digi) const {
226  ESPedestals::const_iterator it_ped = peds_->find(digi.id());
227 
228  ESIntercalibConstantMap::const_iterator it_mip = mips_->getMap().find(digi.id());
229  ESAngleCorrectionFactors::const_iterator it_ang = ang_->getMap().find(digi.id());
230 
231  ESChannelStatusMap::const_iterator it_status = channelStatus_->getMap().find(digi.id());
232 
233  double* results;
234 
235  results = oldEvalAmplitude(digi, it_ped->getMean(), w0_, w1_, w2_);
236 
237  double energy = results[0];
238  double t0 = results[1];
239  int status = (int)results[2];
240  double otenergy = results[3] * 1000000.; // set out-of-time energy to keV
241  delete[] results;
242 
243  double mipCalib = (fabs(cos(*it_ang)) != 0.) ? (*it_mip) / fabs(cos(*it_ang)) : 0.;
244  energy *= (mipCalib != 0.) ? MIPGeV_ / mipCalib : 0.;
245  otenergy *= (mipCalib != 0.) ? MIPGeV_ / mipCalib : 0.;
246 
247  LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : reconstructed energy " << energy;
248 
249  EcalRecHit rechit(digi.id(), energy, t0);
250  // edm: this is just a placeholder for alternative energy reconstruction,
251  // so put it in the same float, with different name
252  // rechit.setOutOfTimeEnergy(otenergy);
253  rechit.setEnergyError(otenergy);
254 
255  if (it_status->getStatusCode() == 1) {
256  rechit.setFlag(EcalRecHit::kESDead);
257  } else {
258  if (status == 0)
259  rechit.setFlag(EcalRecHit::kESGood);
260  else if (status == 5)
261  rechit.setFlag(EcalRecHit::kESBadRatioFor12);
262  else if (status == 6)
263  rechit.setFlag(EcalRecHit::kESBadRatioFor23Upper);
264  else if (status == 7)
265  rechit.setFlag(EcalRecHit::kESBadRatioFor23Lower);
266  else if (status == 8)
267  rechit.setFlag(EcalRecHit::kESTS1Largest);
268  else if (status == 9)
269  rechit.setFlag(EcalRecHit::kESTS3Largest);
270  else if (status == 10)
271  rechit.setFlag(EcalRecHit::kESTS3Negative);
272  else if (status == 11)
273  rechit.setFlag(EcalRecHit::kESSaturated);
274  else if (status == 12)
275  rechit.setFlag(EcalRecHit::kESTS2Saturated);
276  else if (status == 13)
277  rechit.setFlag(EcalRecHit::kESTS3Saturated);
278  else if (status == 14)
279  rechit.setFlag(EcalRecHit::kESTS13Sigmas);
280  }
281 
282  return rechit;
283 }
EcalRecHit
Definition: EcalRecHit.h:15
w2
common ppss p3p6s2 common epss epspn46 common const1 w2
Definition: inclppp.h:1
mps_fire.i
i
Definition: mps_fire.py:428
EcalRecHit::kESBadRatioFor23Lower
Definition: EcalRecHit.h:52
MessageLogger.h
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
ESSample::adc
int adc() const
get the ADC sample (singed 16 bits)
Definition: ESSample.h:16
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
mps_update.status
status
Definition: mps_update.py:69
gpuClustering::adc
uint16_t *__restrict__ uint16_t const *__restrict__ adc
Definition: gpuClusterChargeCut.h:20
ESCondObjectContainer::find
const_iterator find(uint32_t rawId) const
Definition: ESCondObjectContainer.h:33
ESDataFrame::size
int size() const
Definition: ESDataFrame.h:21
ESRecHitSimAlgo::w2_
float w2_
Definition: ESRecHitSimAlgo.h:42
ESDetId::hashedIndex
int hashedIndex() const
get a compact index for arrays [TODO: NEEDS WORK]
Definition: ESDetId.cc:21
EcalRecHit::kESTS3Saturated
Definition: EcalRecHit.h:58
bookConverter.results
results
Definition: bookConverter.py:144
ESRecHitSimAlgo::w1_
float w1_
Definition: ESRecHitSimAlgo.h:41
ESRecHitRatioCuts::getR12High
float getR12High() const
Definition: ESRecHitRatioCuts.h:19
ESRecHitSimAlgo::ang_
const ESAngleCorrectionFactors * ang_
Definition: ESRecHitSimAlgo.h:39
ESRecHitSimAlgo::mips_
const ESIntercalibConstants * mips_
Definition: ESRecHitSimAlgo.h:36
ESDataFrame
Definition: ESDataFrame.h:10
ESDataFrame::sample
const ESSample & sample(int i) const
Definition: ESDataFrame.h:24
RandomServiceHelper.t1
t1
Definition: RandomServiceHelper.py:256
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
ESRecHitSimAlgo::MIPGeV_
float MIPGeV_
Definition: ESRecHitSimAlgo.h:43
FrontierCondition_GT_autoExpress_cfi.t0
t0
Definition: FrontierCondition_GT_autoExpress_cfi.py:149
w
const double w
Definition: UKUtility.cc:23
ESCondObjectContainer< ESPedestal >::const_iterator
std::vector< Item >::const_iterator const_iterator
Definition: ESCondObjectContainer.h:17
HCALHighEnergyHPDFilter_cfi.energy
energy
Definition: HCALHighEnergyHPDFilter_cfi.py:5
ESRecHitSimAlgo::evalAmplitude
EcalRecHit::ESFlags evalAmplitude(float *result, const ESDataFrame &digi, float ped) const
Definition: ESRecHitSimAlgo.cc:12
ESRecHitSimAlgo.h
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
EcalRecHit::kESSaturated
Definition: EcalRecHit.h:56
EcalRecHit::kESTS1Largest
Definition: EcalRecHit.h:53
EcalRecHit::kESTS3Negative
Definition: EcalRecHit.h:55
ESCondObjectContainer::preshower
const Item & preshower(size_t hashedIndex) const
Definition: ESCondObjectContainer.h:25
EcalRecHit::kESTS13Sigmas
Definition: EcalRecHit.h:59
ESRecHitSimAlgo::peds_
const ESPedestals * peds_
Definition: ESRecHitSimAlgo.h:35
createfilelist.int
int
Definition: createfilelist.py:10
ESRecHitSimAlgo::oldreconstruct
EcalRecHit oldreconstruct(const ESDataFrame &digi) const
Definition: ESRecHitSimAlgo.cc:225
ESRecHitSimAlgo::ratioCuts_
const ESRecHitRatioCuts * ratioCuts_
Definition: ESRecHitSimAlgo.h:38
EcalRecHit::ESFlags
ESFlags
Definition: EcalRecHit.h:44
EcalRecHit::kESBadRatioFor12
Definition: EcalRecHit.h:50
ESCondObjectContainer::getMap
const self & getMap() const
Definition: ESCondObjectContainer.h:41
ESRecHitRatioCuts::getR23High
float getR23High() const
Definition: ESRecHitRatioCuts.h:21
ESRecHitSimAlgo::oldEvalAmplitude
double * oldEvalAmplitude(const ESDataFrame &digi, const double &ped, const double &w0, const double &w1, const double &w2) const
Definition: ESRecHitSimAlgo.cc:156
EcalRecHit::setEnergyError
void setEnergyError(float energy)
Definition: EcalRecHit.h:144
ESRecHitSimAlgo::channelStatus_
const ESChannelStatus * channelStatus_
Definition: ESRecHitSimAlgo.h:37
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
ESRecHitSimAlgo::reconstruct
EcalRecHit reconstruct(const ESDataFrame &digi) const
Definition: ESRecHitSimAlgo.cc:82
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
ESRecHitSimAlgo::w0_
float w0_
Definition: ESRecHitSimAlgo.h:40
EcalRecHit::kESTS2Saturated
Definition: EcalRecHit.h:57
EcalRecHit::kESBadRatioFor23Upper
Definition: EcalRecHit.h:51
EcalRecHit::kESDead
Definition: EcalRecHit.h:46
EcalRecHit::kESGood
Definition: EcalRecHit.h:45
ESDataFrame::id
const ESDetId & id() const
Definition: ESDataFrame.h:19
EcalRecHit::kESTS3Largest
Definition: EcalRecHit.h:54
Generated for CMSSW Reference Manual by doxygen 1.8.16