Virtual function for the correction-specific code.
65{
66
67
68
69 float the_aeta;
70 if (m_use_raw_eta (myctx))
71 the_aeta = std::abs (adj_eta);
72 else
73 the_aeta = std::abs (
eta);
74
75 const float etamax = m_etamax (myctx);
76 if (the_aeta >= etamax) return;
77
78 const float eta_start_crack = m_eta_start_crack (myctx);
79 const float eta_end_crack = m_eta_end_crack (myctx);
80 int si;
81 if (the_aeta < eta_start_crack)
82 si = 0;
83 else if (the_aeta > eta_end_crack)
84 si = 1;
85 else {
86
87 return;
88 }
89
91
92
93
94
95 ATH_MSG_DEBUG(
"************************************************************************************************" <<
endmsg);
98 ATH_MSG_DEBUG(
"************************************************************************************************" <<
endmsg);
99
100 unsigned int shape[] = {2};
101 CaloRec::WritableArrayData<1> interp_barriers (shape);
102 interp_barriers[0] = eta_start_crack;
103 interp_barriers[1] = eta_end_crack;
104
105 int ibin = (static_cast<int> (the_aeta / etamax * 100)) ;
106
107 int ibin_frontCorr = ibin;
108 if (m_fix_v6_pathologies (myctx)) {
109 if (the_aeta>2.35) ibin_frontCorr = (static_cast<int> (2.35 / etamax * 100)) ;
110 }
111
112
113
114
115
118 CaloRec::Array<1> lleak =
correction[2][ibin];
119 CaloRec::Array<1> froffset =
correction[3][ibin_frontCorr];
120 CaloRec::Array<1> frslope =
correction[4][ibin_frontCorr];
121 CaloRec::Array<1> sec =
correction[5][ibin_frontCorr];
122
123 ATH_MSG_DEBUG(
"Check etas -------------------------------------------------------------------"<<
endmsg);
124 ATH_MSG_DEBUG(
"Eta --> " << the_aeta <<
" Bin --> " << ibin <<
" Cluster eta = " << cluster->
eta() <<
endmsg);
126
127 ATH_MSG_DEBUG(
"Check calibration coefficients -----------------------------------------------"<<
endmsg);
129 <<
" " << acc[2] <<
" " << acc[3] <<
endmsg);
131 <<
" " << ooc[2] <<
" " << ooc[3] <<
endmsg);
133 <<
" " << lleak[2] <<
" " << lleak[3] <<
endmsg);
135 <<froffset[1] <<
" " << froffset[2] <<
" " << froffset[3] <<
endmsg);
136 ATH_MSG_DEBUG(
"Front Slope : " << frslope[0] <<
" " << frslope[1]
137 <<
" " << frslope[2] <<
" " << frslope[3] <<
endmsg);
139 <<
" " << sec[2] <<
" " << sec[3] <<
endmsg);
140
142 { CaloSampling::PreSamplerB,
143 CaloSampling::EMB1,
144 CaloSampling::EMB2,
145 CaloSampling::EMB3 },
146 { CaloSampling::PreSamplerE,
147 CaloSampling::EME1,
148 CaloSampling::EME2,
149 CaloSampling::EME3 }
150 };
151
152
153
155 eta_start_crack,
156 eta_end_crack,
158 etamax,
160 if (shower_lbary == 0) return;
161
162 if (shower_lbary < 5. || shower_lbary > 25.) {
163 shower_lbary =15.;
165 }
166
167
168
169
170 double eacc_base = 0;
171 for (int sampling=1; sampling<4; sampling++) {
172 eacc_base += cluster->
eSample(samps[si][sampling]);
173 ATH_MSG_DEBUG(
"Barrel/endcap = " << si <<
" Sampling = " << sampling <<
" Energy -->> " << cluster->
eSample(samps[si][sampling]) <<
endmsg);
174 }
175 double eps_base = cluster->
eSample (samps[si][0]);
176
177 ATH_MSG_DEBUG(
"E accordion base --->>>> " << eacc_base <<
" Eps base " << eps_base <<
endmsg);
178
179
180
181
182
183 double depth_max = 20. + (eacc_base+eps_base)*(3./TeV) ;
185 ATH_MSG_DEBUG(
"Bary max for this event ---->> " << depth_max
187
188 if ( shower_lbary > depth_max ) {
189 shower_lbary = 15.;
190
193 }
194
195
196
197
198
199
200
201 double acc_corr =
202 acc[1] +
acc[2] * shower_lbary +
acc[3] * shower_lbary * shower_lbary ;
203
204 double e_out_perc = 0;
205 if (the_aeta < eta_start_crack) {
206 e_out_perc = ooc[1] + ooc[2] * shower_lbary + ooc[3] * shower_lbary * shower_lbary ;
207 }
208 else if (the_aeta > eta_end_crack) {
209 e_out_perc = ooc[1] + ooc[2] * shower_lbary + ooc[3] / shower_lbary ;
210 }
211
212 double e_acc_reco=acc_corr*(eacc_base )*(1+(e_out_perc)*0.01);
213
214
215
216
217
218 double e_leak_perc = 0;
219
220
221
222
223
224
225
226 e_leak_perc = lleak[3] + lleak[1] * shower_lbary + lleak[2] *
exp(shower_lbary);
227
228 if (e_leak_perc < 0 ) e_leak_perc = 0.;
229 if (e_leak_perc > 100.) e_leak_perc = 100.;
230 double e_leak_reco = e_leak_perc * (e_acc_reco)*0.01;
231
232
233
234
235
236
237
238 double raw_energy=e_acc_reco*1
e-3;
239 double e_front_reco = eps_base;
240
241 if (raw_energy <= 1) {
242
243
244 }
245 else if (the_aeta < 1.8) {
246
247 if (the_aeta < eta_start_crack) {
248 double WpsOff = froffset[1] + froffset[2] * raw_energy +
249 froffset[3] * raw_energy * raw_energy;
250 double WpsSlo = frslope[1] *
pow(
log(raw_energy),
251 static_cast<double>(frslope[2])) +
252 frslope[3] *sqrt( raw_energy );
253 e_front_reco=WpsOff + WpsSlo*(eps_base );
254
256 ATH_MSG_DEBUG(
" froffset coeff " << froffset[1] <<
" " << froffset[2] <<
" " << froffset[3] <<
endmsg);
257 ATH_MSG_DEBUG(
" frslope coeff " << frslope[1] <<
" " << frslope[2] <<
" " << frslope[3] <<
endmsg);
260 }
261 else{
262
263 if (raw_energy<20.) raw_energy=20.;
264
265 double WpsOff = froffset[1] + froffset[2] * raw_energy +
266 froffset[3] * sqrt(raw_energy);
267 double WpsSlo = frslope[1] + frslope[2] *
log(raw_energy) +
268 frslope[3] * sqrt(raw_energy);
269 double WpsSlo2 = sec[1] + sec[2] * raw_energy -
270 sec[3] / (raw_energy * raw_energy) ;
271 e_front_reco=WpsOff + WpsSlo*(eps_base) + WpsSlo2*(eps_base)*(eps_base);
272 if (e_front_reco<0.) e_front_reco= eps_base;
273
276 << froffset[2] <<
" " << froffset[3] <<
" " << WpsOff <<
endmsg);
278 << frslope[2] <<
" " << frslope[3] <<
" " << WpsSlo <<
endmsg);
280 << sec[2] <<
" " << sec[3] <<
" " << WpsSlo2 <<
endmsg);
281 ATH_MSG_DEBUG(
" WpsOff, WpsSlo, WpsSlo2 " << WpsOff <<
" " << WpsSlo <<
" " << WpsSlo2 <<
endmsg);
283
284 }
285
286 }
287 else {
288
289 double p1 = froffset[1] + froffset[2] * raw_energy +
290 froffset[3]* raw_energy * raw_energy ;
291 double p2 = frslope[1] + frslope[2] * raw_energy +
292 frslope[3] * raw_energy * raw_energy ;
293 double p3 = sec[1] + sec[2] * raw_energy +
294 sec[3] * raw_energy * raw_energy ;
295
297 << froffset[2] <<
" " << froffset[3] <<
endmsg);
299 << frslope[2] <<
" " << frslope[3] <<
endmsg);
301 << sec[2] <<
" " << sec[3] <<
endmsg);
302
303 e_front_reco= (
p1 +
p2 * shower_lbary +
p3 * shower_lbary * shower_lbary);
304 if (e_front_reco<0.) e_front_reco=eps_base;
305 }
306
307
308
309
310
311 double e_calo_reco =e_front_reco + e_leak_reco + e_acc_reco;
312
313 ATH_MSG_DEBUG(
"CaloSwCalibrationHits::Final reco energy ---------------------- " << e_calo_reco <<
endmsg);
314 ATH_MSG_DEBUG(
"CaloSwCalibrationHits::Front ---------------------- " << e_front_reco <<
endmsg);
315 ATH_MSG_DEBUG(
"CaloSwCalibrationHits::Accordion ------------------ " << e_acc_reco <<
endmsg);
316 ATH_MSG_DEBUG(
"CaloSwCalibrationHits::out of cone ---------------- " << acc_corr*(eacc_base )*(e_out_perc)*0.01 <<
endmsg);
317 ATH_MSG_DEBUG(
"CaloSwCalibrationHits::Leakage -------------------- " << e_leak_reco <<
endmsg);
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334 if (m_updateSamplingEnergies (myctx))
335 {
336
337
338 cluster->
setEnergy (samps[si][0], e_front_reco );
339
340
341
342 for (int sampling=1; sampling<=2; sampling++){
343 cluster->
setEnergy (samps[si][sampling], cluster->
eSample(samps[si][sampling]) * acc_corr );
344 }
345
346
347
348 cluster->
setEnergy (samps[si][3], cluster->
eSample (samps[si][3]) * acc_corr + e_leak_reco);
349 }
350
351
352 double e_temp = 0;
353 for (
int nl = 0 ; nl< 4 ; nl++) e_temp += cluster->
eSample (samps [si][nl]);
354
355 ATH_MSG_DEBUG(
"---------- Sum of the sampling energy --- >> " << e_temp <<
" EcaloReco = " << e_calo_reco <<
endmsg);
356
357 cluster->
setE (e_calo_reco);
358
359 ATH_MSG_DEBUG(
"CaloSwCalibHitsCalibration::Energy after correction --> " << cluster->
e() <<
endmsg);
360
361}
constexpr int pow(int base, int exp) noexcept
static double depth(const float aeta, const float start_crack, const float end_crack, const CaloRec::Array< 2 > &sampling_depth, const float etamax, const xAOD::CaloCluster *cluster, MsgStream &log)
Calculate the depth of the cluster.
Constant< CxxUtils::Array< 2 > > m_sampling_depth
Read-only multidimensional array.
float eSample(const CaloSample sampling) const
bool setEnergy(const CaloSample sampling, const float e)
Set energy for a given sampling. Returns false if the sample isn't part of the cluster.