161 {
162 const double dR_CUT2 = dR_CUT * dR_CUT;
163
164 double etot_em{0}, etot_hd{0}, etot{0};
165
166 double s[24] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
167
170 double dphi = std::abs(theClus->phi() - phi_trkAtCalo);
171 if (dphi >
M_PI) dphi = 2 *
M_PI - dphi;
172 const double deta = std::abs(theClus->eta() - eta_trkAtCalo);
173 const double rij2 = deta * deta + dphi * dphi;
174 if (rij2 > dR_CUT2) continue;
175
176 double eemb0 = theClus->eSample(CaloSampling::PreSamplerB);
177 double eemb1 = theClus->eSample(CaloSampling::EMB1);
178 double eemb2 = theClus->eSample(CaloSampling::EMB2);
179 double eemb3 = theClus->eSample(CaloSampling::EMB3);
180 double eeme0 = theClus->eSample(CaloSampling::PreSamplerE);
181 double eeme1 = theClus->eSample(CaloSampling::EME1);
182 double eeme2 = theClus->eSample(CaloSampling::EME2);
183 double eeme3 = theClus->eSample(CaloSampling::EME3);
184 double etileg1 = theClus->eSample(CaloSampling::TileGap1);
185 double etileg2 = theClus->eSample(CaloSampling::TileGap2);
186 double etileg3 = theClus->eSample(CaloSampling::TileGap3);
187
199
200 etot_em += eemb0 + eemb1 + eemb2 + eemb3 + eeme0 + eeme1 + eeme2 + eeme3 + etileg1 + etileg2 + etileg3;
201
202 double etileb0 = theClus->eSample(CaloSampling::TileBar0);
203 double etileb1 = theClus->eSample(CaloSampling::TileBar1);
204 double etileb2 = theClus->eSample(CaloSampling::TileBar2);
205 double etilee0 = theClus->eSample(CaloSampling::TileExt0);
206 double etilee1 = theClus->eSample(CaloSampling::TileExt1);
207 double etilee2 = theClus->eSample(CaloSampling::TileExt2);
208 double ehec0 = theClus->eSample(CaloSampling::HEC0);
209 double ehec1 = theClus->eSample(CaloSampling::HEC1);
210 double ehec2 = theClus->eSample(CaloSampling::HEC2);
211 double ehec3 = theClus->eSample(CaloSampling::HEC3);
212 double efcal0 = theClus->eSample(CaloSampling::FCAL0);
213 double efcal1 = theClus->eSample(CaloSampling::FCAL1);
214 double efcal2 = theClus->eSample(CaloSampling::FCAL2);
215
229
230 etot_hd += etileb0 + etileb1 + etileb2 + etilee0 + etilee1 + etilee2 + ehec0 + ehec1 + ehec2 + ehec3 + efcal0 + efcal1 + efcal2;
231
232 etot += eemb0 + eemb1 + eemb2 + eemb3 + eeme0 + eeme1 + eeme2 + eeme3 + etileg1 + etileg2 + etileg3 + etileb0 + etileb1 + etileb2 +
233 etilee0 + etilee1 + etilee2 + ehec0 + ehec1 + ehec2 + ehec3 + efcal0 + efcal1 + efcal2;
234 }
235
236 ATH_MSG_VERBOSE(
"Values extracted from the CaloTopoClusters for a cone of: "
237 << dR_CUT << "\n"
238 << " - Energy in Calorimeter Total: " << etot << " EM: " << etot_em << " HAD: " << etot_hd
239 << "\n eemb0: " << s[0] << "\n eemb1: " << s[1] << "\n eemb2: " << s[2] << "\n eemb3: " << s[3]
240 << "\n eeme0: " << s[4] << "\n eeme1: " << s[5] << "\n eeme2: " << s[6] << "\n eeme3: " << s[7]
241 << "\n etileg1: " << s[8] << "\n etileg2: " << s[9] << "\n etileg3: " << s[10] << "\n etileb0: " << s[11]
242 << "\n etileb1: " << s[12] << "\n etileb2: " << s[13] << "\n etilee0: " << s[14] << "\n etilee1: " << s[15]
243 << "\n etilee2: " << s[16] << "\n ehec0: " << s[17] << "\n ehec1: " << s[18] << "\n ehec2: " << s[19]
244 << "\n ehec3: " << s[20] << "\n efcal0: " << s[21] << "\n efcal1: " << s[22]
245 << "\n efcal2: " << s[23]);
246
247 for (
int i = 0;
i < 24; ++
i) s[i] /= Gaudi::Units::GeV;
248 etot_em /= Gaudi::Units::GeV;
249 etot_hd /= Gaudi::Units::GeV;
250 etot /= Gaudi::Units::GeV;
251
252 double tmp_mx(-999), tmp_mxH(-999), tmp_mxE(-999);
253 for (
int i = 0;
i < 24; ++
i) {
254 if (s[i] > tmp_mx) { tmp_mx =
s[
i]; }
255 if (i < 11 && s[i] > tmp_mxE ) { tmp_mxE =
s[
i]; }
256 if (i >= 11 && s[i] > tmp_mxH ) { tmp_mxH =
s[
i]; }
257 }
258
259 double emFr{999}, mxFr{999}, mxEM{999}, mxHad{999}, EoverEtrk{999}, eemb1_wrtTotal{999}, eemb2_wrtTotal{999}, eemb3_wrtGroup{999},
260 etileb0_wrtGroup{999}, etileb1_wrtTotal{999}, etileb2_wrtGroup{999}, eeme1_wrtGroup{999}, eeme1_wrtTotal{999}, eeme2_wrtTotal{999},
261 eeme3_wrtGroup{999}, ehec0_wrtTotal{999};
262
263 if (etot > 0) {
264 mxFr = tmp_mx / etot;
265 mxEM = tmp_mxE / etot;
266 mxHad = tmp_mxH / etot;
267 emFr = etot_em / etot;
268 eemb1_wrtTotal =
s[1] / etot;
269 eemb2_wrtTotal =
s[2] / etot;
270 etileb1_wrtTotal =
s[12] / etot;
271 eeme1_wrtTotal =
s[5] / etot;
272 eeme2_wrtTotal =
s[6] / etot;
273 ehec0_wrtTotal =
s[17] / etot;
274 }
275
276
277 if (p_trk) EoverEtrk = etot / (p_trk / Gaudi::Units::GeV);
278
279 if (etot_em > 0) {
280 eemb3_wrtGroup =
s[3] / etot_em;
281 eeme1_wrtGroup =
s[5] / etot_em;
282 eeme3_wrtGroup =
s[7] / etot_em;
283 }
284
285 if (etot_hd > 0) {
286 etileb0_wrtGroup =
s[11] / etot_hd;
287 etileb2_wrtGroup =
s[13] / etot_hd;
288 }
289
290
291 std::map<std::string, double> vars;
292 vars["emFr"] = emFr;
293 vars["EoverEtrk"] = EoverEtrk;
294 vars["mxFr"] = mxFr;
295 vars["mxEM"] = mxEM;
296 vars["mxHad"] = mxHad;
297 vars["eemb1_wrtTotal"] = eemb1_wrtTotal;
298 vars["eemb2_wrtTotal"] = eemb2_wrtTotal;
299 vars["eemb3_wrtGroup"] = eemb3_wrtGroup;
300 vars["etileb0_wrtGroup"] = etileb0_wrtGroup;
301 vars["etileb1_wrtTotal"] = etileb1_wrtTotal;
302 vars["etileb2_wrtGroup"] = etileb2_wrtGroup;
303 vars["eeme1_wrtGroup"] = eeme1_wrtGroup;
304 vars["eeme1_wrtTotal"] = eeme1_wrtTotal;
305 vars["eeme2_wrtTotal"] = eeme2_wrtTotal;
306 vars["eeme3_wrtGroup"] = eeme3_wrtGroup;
307 vars["ehec0_wrtTotal"] = ehec0_wrtTotal;
308
309 ATH_MSG_DEBUG(
"likelihood discriminant variables values: " << dR_CUT);
310 if (msgLevel(MSG::DEBUG)) {
311 std::map<std::string, double>::const_iterator
iter;
312 for (iter = vars.begin(); iter != vars.end(); ++iter)
ATH_MSG_DEBUG(
" - " <<
iter->first <<
": " <<
iter->second);
313 }
314 double LR{0}, ProbS{1}, ProbB{1};
315
316 int iFile = -1;
317 if (std::abs(eta_trk) < 1.4) {
318 if (p_trk / Gaudi::Units::GeV < 11.)
319 iFile = 0;
320 else if (p_trk / Gaudi::Units::GeV < 51.)
321 iFile = 1;
322 else
323 iFile = 2;
324 } else if (std::abs(eta_trk) >= 1.4 && std::abs(eta_trk) <= 1.6) {
325 if (p_trk / Gaudi::Units::GeV < 11.)
326 iFile = 3;
327 else if (p_trk / Gaudi::Units::GeV < 51.)
328 iFile = 4;
329 else
330 iFile = 5;
331 } else if (std::abs(eta_trk) > 1.6 && std::abs(eta_trk) < 2.5) {
332 if (p_trk / Gaudi::Units::GeV < 11.)
333 iFile = 6;
334 else if (p_trk / Gaudi::Units::GeV < 51.)
335 iFile = 7;
336 else
337 iFile = 8;
338 }
339
340 if (iFile < 0) return LR;
341
343 ProbS = 0;
344 ProbB = 1;
345 } else {
347 std::map<std::string, double>::const_iterator
it = vars.find(
m_TH1F_key[iFile][i]);
348
350 <<
": m_TH1F_key[" << iFile <<
"][" << i <<
"(/" <<
m_numKeys[iFile] <<
")] = " <<
m_TH1F_key[iFile][i] <<
", "
351 << ((it != vars.end()) ? "found" : "not found"));
352
353 if (it != vars.end()) {
354 int bin_sig =
m_TH1F_sig[iFile][
i]->GetXaxis()->FindFixBin(
it->second);
355
357 <<
": it->first = " <<
it->first <<
", it->second = " <<
it->second <<
", bin_sig = " << bin_sig
358 <<
", NbinsX = " <<
m_TH1F_sig[iFile][i]->GetNbinsX());
359
360 if (bin_sig >= 1 && bin_sig <=
m_TH1F_sig[iFile][i]->GetNbinsX()) {
361 double SbinContent =
m_TH1F_sig[iFile][
i]->GetBinContent(bin_sig);
362 ATH_MSG_DEBUG(
"m_TH1F_sig Bin Content for " <<
it->first <<
": " << SbinContent);
363 if (SbinContent)
364 ProbS *= SbinContent;
365 else {
366 ProbS *= 0.0000001;
368 }
369 } else {
371 ProbS *= 1;
372 }
374 int bin_bkg =
m_TH1F_bkg[iFile][
i]->GetXaxis()->FindFixBin(
it->second);
375
377 <<
": it->first = " <<
it->first <<
", it->second = " <<
it->second <<
", bin_bkg = " << bin_bkg
378 <<
", NbinsX = " <<
m_TH1F_bkg[iFile][i]->GetNbinsX());
379
380 if (bin_bkg >= 1 && bin_bkg <=
m_TH1F_bkg[iFile][i]->GetNbinsX()) {
381 double BbinContent =
m_TH1F_bkg[iFile][
i]->GetBinContent(bin_bkg);
382 ATH_MSG_VERBOSE(
"m_TH1F_bkg Bin Content for " <<
it->first <<
": " << BbinContent);
383 if (BbinContent)
384 ProbB *= BbinContent;
385 else {
386 ProbB *= 0.0000001;
388 }
389 } else {
391 ProbB *= 1;
392 }
394 } else {
397 << "> is not found in the calculated variable list");
399 }
400 }
401 }
402 }
403
404 if (ProbS == 0)
405 LR = 0;
406 else
407 LR = ProbS / (ProbS + ProbB);
408
412 return LR;
413}
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
CaloCluster_v1 CaloCluster
Define the latest version of the calorimeter cluster class.