169{
170
173 ATH_MSG_WARNING(
"Model not loaded or Egamma pointer is null, returning cluster energy");
175 } else {
176 ATH_MSG_FATAL(
"Model not loaded or Egamma pointer is null, and useClusterIf0 is false, cannot proceed");
177 return 0.0f;
178 }
179 }
180
181
182 IegammaCellRecoveryTool::Info recoveryInfo;
183 bool recoverySucceeded = true;
185 egammaCellUtils::MaxECell maxECell(&clus);
186 if (maxECell.sc == StatusCode::FAILURE) {
188 recoverySucceeded = false;
189 } else {
190 recoveryInfo.
etamax = maxECell.etaCell;
191 recoveryInfo.
phimax = maxECell.phiCell;
193 ATH_MSG_WARNING(
"Cell Recovery Tool failed. Proceeding without recovered cells.");
194 recoverySucceeded = false;
195 }
196 }
197 }
198
199
201 auto array_layer_scales = std::array<double, 4>{1.0, 1.0, 1.0, 1.0};
202
204 ATH_MSG_DEBUG(
"Applying layer recalibration for GNN on data.");
205
206
209 }
210
212 ATH_MSG_DEBUG(
"Applying extra layer scales for systematic studies, normally this is for MC events.");
214 ATH_MSG_WARNING(
"You are applying extra layer scales but the input is not MC! Are you sure this is intended?");
215 }
216
217 for (std::size_t i = 0;
i < 4; ++
i)
218 array_layer_scales[i] *= gei.
scaleEs[i];
219 }
220
221
222
223
224
228
229
230 std::vector<float> cells_E, cells_eta, cells_phi, cells_x, cells_y, cells_z;
231 std::vector<int> cells_layer;
232 std::vector<Identifier> included_cells;
233
234
235 double sum_cell_E_L0 = 0.0, sum_cell_E_L1 = 0.0, sum_cell_E_L2 = 0.0, sum_cell_E_L3 = 0.0, sum_cell_E_Gap = 0.0;
236
237
238 const CaloClusterCellLink* cellLinks = clus.
getCellLinks();
239 if (cellLinks) {
240 for (const CaloCell* cell : *cellLinks) {
241 if (!cell) continue;
242
243 int sampling =
cell->caloDDE()->getSampling();
244 double scale_factor = 1.0;
246
247 switch (sampling) {
248 case CaloCell_ID::PreSamplerB: case CaloCell_ID::PreSamplerE:
249 scale_factor = array_layer_scales[0];
layer_idx = 0;
break;
250 case CaloCell_ID::EMB1: case CaloCell_ID::EME1:
251 scale_factor = array_layer_scales[1];
layer_idx = 1;
break;
252 case CaloCell_ID::EMB2: case CaloCell_ID::EME2:
253 scale_factor = array_layer_scales[2];
layer_idx = 2;
254
256 included_cells.push_back(
cell->ID());
257 }
258 break;
259 case CaloCell_ID::EMB3: case CaloCell_ID::EME3:
260 scale_factor = array_layer_scales[3];
layer_idx = 3;
262 included_cells.push_back(
cell->ID());
263 }
264 break;
265 case CaloCell_ID::TileGap3:
266 scale_factor = 1.0;
layer_idx = 4;
break;
267 default: continue;
268 }
269
270 double final_E =
cell->e() * scale_factor;
271
272 cells_E.push_back(final_E);
273 cells_eta.push_back(
cell->eta());
274 cells_phi.push_back(
cell->phi());
275 cells_x.push_back(
cell->x());
276 cells_y.push_back(
cell->y());
277 cells_z.push_back(
cell->z());
278 cells_layer.push_back(layer_idx);
279
280
281 switch(layer_idx) {
282 case 0: sum_cell_E_L0 += final_E; break;
283 case 1: sum_cell_E_L1 += final_E; break;
284 case 2: sum_cell_E_L2 += final_E; break;
285 case 3: sum_cell_E_L3 += final_E; break;
286 case 4: sum_cell_E_Gap += final_E; break;
287 }
288 }
289 }
290
291
292
293 for (
const CaloCell* cell : recoveryInfo.
addedCells) {
294 if (!cell || !
cell->caloDDE())
continue;
295
296
297 if (std::find(included_cells.begin(), included_cells.end(),
cell->ID()) != included_cells.end()) {
298 ATH_MSG_WARNING(
"Recovered cell " <<
cell->ID() <<
" already included in cluster. Skipping to avoid double counting.");
299 continue;
300 }
301 else {
303 }
304
305 int sampling =
cell->caloDDE()->getSampling();
306 double scale_factor = 1.0;
308
309 if (sampling == CaloCell_ID::EMB2 || sampling == CaloCell_ID::EME2) {
310 scale_factor = array_layer_scales[2];
layer_idx = 2;
311 } else if (sampling == CaloCell_ID::EMB3 || sampling == CaloCell_ID::EME3) {
312 scale_factor = array_layer_scales[3];
layer_idx = 3;
313 } else {
314
315 if (sampling == CaloCell_ID::PreSamplerB || sampling == CaloCell_ID::PreSamplerE) {
316 scale_factor = array_layer_scales[0];
layer_idx = 0;
317 } else if (sampling == CaloCell_ID::EMB1 || sampling == CaloCell_ID::EME1) {
318 scale_factor = array_layer_scales[1];
layer_idx = 1;
319 } else {
320 continue;
321 }
322 }
323
324 double final_E =
cell->e() * scale_factor;
325
326 cells_E.push_back(final_E);
327 cells_eta.push_back(
cell->eta());
328 cells_phi.push_back(
cell->phi());
329 cells_x.push_back(
cell->x());
330 cells_y.push_back(
cell->y());
331 cells_z.push_back(
cell->z());
332 cells_layer.push_back(layer_idx);
333
334 switch(layer_idx) {
335 case 0: sum_cell_E_L0 += final_E; break;
336 case 1: sum_cell_E_L1 += final_E; break;
337 case 2: sum_cell_E_L2 += final_E; break;
338 case 3: sum_cell_E_L3 += final_E; break;
339 case 4: sum_cell_E_Gap += final_E; break;
340 }
341 }
342
343
344 const size_t nCells = cells_E.size();
345 if (nCells == 0) return 0.0f;
346
347 double sum_cell_E_total = sum_cell_E_L0 + sum_cell_E_L1 + sum_cell_E_L2 + sum_cell_E_L3;
348 const double cluster_eta = clus.
eta();
349 const double cluster_phi = clus.
phi();
350
351 std::vector<float> cells_deta, cells_dphi, cells_eFrac;
352 cells_deta.reserve(nCells);
353 cells_dphi.reserve(nCells);
354 cells_eFrac.reserve(nCells);
355
356 for (
size_t i = 0;
i <
nCells; ++
i) {
357 float deta = cells_eta[
i] - cluster_eta;
358 float dphi = cells_phi[
i] - cluster_phi;
359 dphi = std::fmod(dphi + 3.0f *
M_PI, 2.0f *
M_PI) -
M_PI;
360
361 cells_deta.push_back(deta);
362 cells_dphi.push_back(dphi);
363
364 float eFrac_layer = 0.0f;
365 switch (cells_layer[i]) {
366 case 0: eFrac_layer = (sum_cell_E_L0 != 0) ? (cells_E[i] / sum_cell_E_L0) : 0.0f; break;
367 case 1: eFrac_layer = (sum_cell_E_L1 != 0) ? (cells_E[i] / sum_cell_E_L1) : 0.0f; break;
368 case 2: eFrac_layer = (sum_cell_E_L2 != 0) ? (cells_E[i] / sum_cell_E_L2) : 0.0f; break;
369 case 3: eFrac_layer = (sum_cell_E_L3 != 0) ? (cells_E[i] / sum_cell_E_L3) : 0.0f; break;
370 case 4: eFrac_layer = (sum_cell_E_Gap != 0) ? (cells_E[i] / sum_cell_E_Gap) : 0.0f; break;
371 }
372 cells_eFrac.push_back(eFrac_layer);
373 }
374
375
376 double ratio_L1_L2 = (sum_cell_E_L2 != 0) ? (sum_cell_E_L1 / sum_cell_E_L2) : 0.0;
377 double main_layers_sum = sum_cell_E_L1 + sum_cell_E_L2 + sum_cell_E_L3;
378 double ratio_L0_total = (main_layers_sum != 0) ? (sum_cell_E_L0 / main_layers_sum) : 0.0;
379 double ratio_Tile_total = (main_layers_sum != 0) ? (sum_cell_E_Gap / main_layers_sum) : 0.0;
380
382
383
384 std::vector<float> cluster_feats = {
385 static_cast<float>(sum_cell_E_total),
386 static_cast<float>(sum_cell_E_L0),
387 static_cast<float>(sum_cell_E_L1),
388 static_cast<float>(sum_cell_E_L2),
389 static_cast<float>(sum_cell_E_L3),
390 static_cast<float>(sum_cell_E_Gap),
391 static_cast<float>(cluster_eta),
392 static_cast<float>(cluster_phi),
393 static_cast<float>(ratio_L1_L2),
394 static_cast<float>(ratio_L0_total),
395 static_cast<float>(ratio_Tile_total)
396 };
397
398
401 if (photon) {
402
403 float convR = 799.0f;
406 }
407
408
409 float convEtOverPt = 0.0f;
413 (raw_Es1 * array_layer_scales[1] + raw_Es2 * array_layer_scales[2] + raw_Es3 * array_layer_scales[3]) :
414 (raw_Es1 + raw_Es2 + raw_Es3));
416 convEtOverPt = std::max(0.0f, eacc / (std::cosh(cl_eta) * ptconv));
417 }
418 convEtOverPt = std::min(convEtOverPt, 2.0f);
419
420
421 float convPtRatio = 1.0f;
425 if ((pt1 + pt2) > 0.0f) {
426 convPtRatio = std::max(pt1, pt2) / (pt1 + pt2);
427 }
428 }
429
430
432
433 cluster_feats.push_back(convR);
434 cluster_feats.push_back(convEtOverPt);
435 cluster_feats.push_back(convPtRatio);
436 cluster_feats.push_back(conversionType);
437 } else {
438 cluster_feats.push_back(0.0f);
439 cluster_feats.push_back(0.0f);
440 cluster_feats.push_back(0.0f);
441 cluster_feats.push_back(0.0f);
442 }
443 }
444
445
446 gnn_input["cluster_features"] = FlavorTagInference::Inputs(cluster_feats, {1, (int64_t)cluster_feats.size()});
447
448
449 std::vector<float> cell_feats_flat;
451 for (
size_t i = 0;
i <
nCells; ++
i) {
452 cell_feats_flat.push_back(cells_eFrac[i]);
453 cell_feats_flat.push_back(cells_deta[i]);
454 cell_feats_flat.push_back(cells_dphi[i]);
455 cell_feats_flat.push_back(cells_x[i]);
456 cell_feats_flat.push_back(cells_y[i]);
457 cell_feats_flat.push_back(cells_z[i]);
458 cell_feats_flat.push_back(static_cast<float>(cells_layer[i]));
459 }
460 gnn_input[
"cell_features"] = FlavorTagInference::Inputs(cell_feats_flat, {(int64_t)nCells,
m_num_cell_features});
461
462
463 auto [out_f, out_vc, out_vf] =
m_saltModel->runInference(gnn_input);
464
465 float el_gnn_score = 0.0f;
466 if (out_vf.empty() || out_vf.begin()->second.empty()) {
468 } else {
469 el_gnn_score = out_vf.begin()->second.front();
470 }
471
472
473 if (el_gnn_score == 0.0f) {
475 }
476
477 return el_gnn_score * static_cast<float>(sum_cell_E_total);
478}
#define ATH_MSG_WARNING(x)
const CaloClusterCellLink * getCellLinks() const
Get a pointer to the CaloClusterCellLink object (const version).
virtual double eta() const
The pseudorapidity ( ) of the particle.
virtual double e() const
The total energy of the particle.
float energyBE(const unsigned layer) const
Get the energy in one layer of the EM Calo.
virtual double phi() const
The azimuthal angle ( ) of the particle.
const xAOD::CaloCluster * caloCluster(size_t index=0) const
Pointer to the xAOD::CaloCluster/s that define the electron candidate.
std::map< std::string, Inputs, std::less<> > InputMap
float compute_ptconv(const xAOD::Photon *ph)
This ptconv is the old one used by MVACalib.
float compute_pt2conv(const xAOD::Photon *ph)
float compute_pt1conv(const xAOD::Photon *ph)
ConversionType conversionType(const bool hasTrk1, const bool hasTrk2, const std::uint8_t nSiHits1, const std::uint8_t nSiHits2)
return the photon conversion type (see EgammaEnums)
std::size_t numberOfSiTracks(const xAOD::Photon *eg)
return the number of Si tracks in the conversion
float conversionRadius(const xAOD::Vertex *vx)
return the conversion radius or 9999.
EventInfo_v1 EventInfo
Definition of the latest event info version.
Photon_v1 Photon
Definition of the current "egamma version".
setRawEt setRawPhi nCells
std::array< float, 4 > scaleEs
const xAOD::EventInfo * eventInfo