ATLAS Offline Software
eFEXtauBDT.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 //*************************************************************************
6 // eFEXtauBDT - description
7 // --------------------
8 // begin : 15 09 2022
9 // email : david.reikher@gmail.com
10 //*************************************************************************
11 
13 #include <string>
14 
15 #define ENERGY_WIDTH 16
16 #define PARAM_WIDTH 8
17 
18 // default constructor for persistency
19 LVL1::eFEXtauBDT::eFEXtauBDT(AthAlgTool *log, std::string config_path)
20  : m_bdt(config_path), m_log(log) {
21  m_log->msg(MSG::DEBUG) << "Configured BDT with this file: " << config_path
22  << endmsg;
23  m_bdtVars.resize(m_bdt.getVariables().size());
24  m_towers.resize(m_bdt.getNTowers());
25  int n_multipliers = sizeof(m_fracMultipliers) / sizeof(m_fracMultipliers[0]);
26  m_emEtXMultiplier.resize(n_multipliers);
27  m_emEtXMultiplierOverflow.resize(n_multipliers);
28 }
29 
32 
33 void LVL1::eFEXtauBDT::setPointerToSCell(int eta, int phi, int layer,
34  unsigned int *sCellPtr) {
35  switch (layer) {
36  case 0:
37  m_em0cells[eta][phi] = sCellPtr;
38  break;
39  case 1:
40  m_em1cells[eta][phi] = sCellPtr;
41  break;
42  case 2:
43  m_em2cells[eta][phi] = sCellPtr;
44  break;
45  case 3:
46  m_em3cells[eta][phi] = sCellPtr;
47  break;
48  case 4:
49  m_hadcells[eta][phi] = sCellPtr;
50  break;
51  }
52 }
53 
55  int index, unsigned int *fracMultiplierPtr) {
56  m_fracMultipliers[index] = fracMultiplierPtr;
57 }
58 
60  int index, unsigned int *bdtThresholdPtr) {
61  m_bdtThresholds[index] = bdtThresholdPtr;
62 }
63 
64 void LVL1::eFEXtauBDT::setPointerToETThresholdParam(unsigned int *etThreshold) {
65  m_etThreshold = etThreshold;
66 }
67 
69  unsigned int *maxEtThreshold) {
70  m_maxEtThreshold = maxEtThreshold;
71 }
72 
74  unsigned int *maxEtThreshold) {
75  m_maxEtThresholdFrac = maxEtThreshold;
76 }
77 
79  unsigned int *bdtMinEtThreshold) {
80  m_bdtMinEtThreshold = bdtMinEtThreshold;
81 }
82 
83 unsigned int *LVL1::eFEXtauBDT::superCellToPtr(int eta, int phi, int layer) {
84  unsigned int *ptr = 0;
85  switch (layer) {
86  case 0:
87  ptr = m_em0cells[eta][phi];
88  break;
89  case 1:
90  ptr = m_em1cells[eta][phi];
91  break;
92  case 2:
93  ptr = m_em2cells[eta][phi];
94  break;
95  case 3:
96  ptr = m_em3cells[eta][phi];
97  break;
98  case 4:
99  ptr = m_hadcells[eta][phi];
100  break;
101  }
102  return ptr;
103 }
104 
106  initPointers(m_bdt.getETSCells(), m_eTComputeSCellPointers);
107 }
108 
110  initPointers(m_bdt.getEMETSCells(), m_EM_eTComputeSCellPointers);
111 }
112 
114  initPointers(m_bdt.getHADETSCells(), m_HAD_eTComputeSCellPointers);
115 }
116 
117 void LVL1::eFEXtauBDT::initPointers(const std::vector<std::vector<int>> &scells,
118  std::vector<unsigned int *> &ptr_list) {
119  m_log->msg(MSG::DEBUG) << "Will use sum of supercells: " << endmsg;
120  for (const auto& scell : scells) {
121  int eta = scell[0];
122  int phi = scell[1];
123  int layer = scell[2];
124  m_log->msg(MSG::DEBUG) << "\teta=" << eta << "\tphi=" << phi
125  << "\tlayer=" << layer << endmsg;
126  unsigned int *ptr = superCellToPtr(eta, phi, layer);
127  if (ptr == 0) {
128  m_log->msg(MSG::DEBUG)
129  << "Could not convert eta=" << eta << " phi=" << phi
130  << " layer=" << layer
131  << " to a pointer to supercell. Are they within range?" << endmsg;
132  throw std::domain_error(
133  std::string("Could not convert eta=") + std::to_string(eta) +
134  " phi=" + std::to_string(phi) + " layer=" + std::to_string(layer) +
135  " to a pointer to supercell. Are they within range?");
136  }
137  ptr_list.push_back(ptr);
138  }
139 }
140 
142  m_towersComputeSCellPointers.resize(m_bdt.getNTowers());
143  for (size_t i = 0; i < m_towers.size(); i++) {
144  m_log->msg(MSG::DEBUG) << "Tower " << i << endmsg;
145  initPointers(m_bdt.getTowerSCells(i), m_towersComputeSCellPointers[i]);
146  }
147 }
148 
150  for (size_t i = 0; i < m_bdt.getVariables().size(); i++) {
151  BDTVariable var = m_bdt.getVariables()[i];
152 
153  m_log->msg(MSG::DEBUG) << i << " is " << var.m_name << ", sum of supercells"
154  << endmsg;
155  std::vector<unsigned int *> pointersToSCells;
156  for (size_t j = 0; j < var.m_scells.size(); j++) {
157  int eta = var.m_scells[j][0];
158  int phi = var.m_scells[j][1];
159  int layer = var.m_scells[j][2];
160  m_log->msg(MSG::DEBUG) << "\teta=" << eta << "\tphi=" << phi
161  << "\tlayer=" << layer << endmsg;
162  unsigned int *ptr = superCellToPtr(eta, phi, layer);
163  if (ptr == 0) {
164  m_log->msg(MSG::DEBUG)
165  << "Could not convert eta=" << eta << " phi=" << phi
166  << " layer=" << layer
167  << " to a pointer to supercell. Are they within range?" << endmsg;
168  throw std::domain_error(
169  std::string("Could not convert eta=") + std::to_string(eta) +
170  " phi=" + std::to_string(phi) + " layer=" + std::to_string(layer) +
171  " to a pointer to supercell. Are they within range?");
172  }
173  pointersToSCells.push_back(ptr);
174  }
175 
176  m_bdtVarComputeSCellPointers.push_back(pointersToSCells);
177  }
178 }
179 
181  buildBDTVariables();
182  computeTowers();
183  computeETEstimate();
184  computeHADETEstimate();
185  computeEMETEstimate();
186  computeBDTCondition();
187  computeFracCondition();
188  computeIsCentralTowerSeed();
189 }
190 
192  std::string bdtVariables = "";
193  for (size_t i = 0; i < m_bdtVars.size(); i++) {
194  bdtVariables += std::to_string(m_bdtVars[i]) + " ";
195  }
196 
197  m_log->msg(MSG::DEBUG) << "BDT Variables: " << bdtVariables << endmsg;
198 }
199 
200 // Build BDT Variables
202  for (size_t i = 0; i < m_bdtVarComputeSCellPointers.size(); i++) {
203  bool overflow;
204  m_bdtVars[i] = computeEstimate(m_bdtVarComputeSCellPointers[i], overflow,
205  ENERGY_WIDTH);
206  if (overflow) {
207  m_bdtVars[i] = (1 << ENERGY_WIDTH) - 1;
208  }
209  }
210  debugPrintBDTVariables();
211  m_bdtVarsComputed = true;
212 }
213 
215  // Need BDT variables to be computed
216  if (m_bdtVarsComputed == false) {
217  m_log->msg(MSG::DEBUG)
218  << "BDT Variables not computed. BDT score will be garbage." << endmsg;
219  }
220 
221  m_bdtScore = m_bdt.getBDT().decision_function(m_bdtVars)[0];
222  m_log->msg(MSG::DEBUG) << "BDT Score: " << m_bdtScore << endmsg;
223 }
224 
226  m_eTEstimate = computeEstimate(m_eTComputeSCellPointers, m_eTEstimateOverflow,
227  ENERGY_WIDTH);
228  m_log->msg(MSG::DEBUG) << "ET Estimate: " << m_eTEstimate << endmsg;
229 }
230 
232  m_EM_eTEstimate = computeEstimate(m_EM_eTComputeSCellPointers,
233  m_EM_eTEstimateOverflow, ENERGY_WIDTH);
234  m_log->msg(MSG::DEBUG) << "EM ET Estimate: " << m_EM_eTEstimate << endmsg;
235 }
236 
238  m_HAD_eTEstimate = computeEstimate(m_HAD_eTComputeSCellPointers,
239  m_HAD_eTEstimateOverflow, ENERGY_WIDTH);
240  m_log->msg(MSG::DEBUG) << "HAD ET Estimate: " << m_HAD_eTEstimate << endmsg;
241 }
242 
244  m_log->msg(MSG::DEBUG) << "Towers Estimate: " << endmsg;
245  for (int eta = 0; eta < 3; eta++) {
246  for (int phi = 0; phi < 3; phi++) {
247  int flatIndex = flatTowerIndex(eta, phi);
248  m_log->msg(MSG::DEBUG)
249  << "Tower " << flatIndex << " ET (eta=" << eta << ", phi=" << phi
250  << "): " << m_towers[flatIndex] << endmsg;
251  }
252  }
253 }
254 
255 int LVL1::eFEXtauBDT::flatTowerIndex(int eta, int phi) { return 3 * phi + eta; }
256 
258  for (size_t i = 0; i < m_towers.size(); i++) {
259  bool overflow;
260  m_towers[i] = computeEstimate(m_towersComputeSCellPointers[i], overflow,
261  ENERGY_WIDTH);
262  if (overflow) {
263  m_towers[i] = (1 << ENERGY_WIDTH) - 1;
264  }
265  }
266 
267  debugPrintTowers();
268 }
269 
270 bool LVL1::eFEXtauBDT::isOverflow(unsigned int number, int nBits) {
271  if ((number >> nBits) != 0) {
272  return true;
273  }
274  return false;
275 }
276 
277 unsigned int
278 LVL1::eFEXtauBDT::computeEstimate(std::vector<unsigned int *> &ptr_list,
279  bool &overflow, int resultNBits) {
280  unsigned int estimate = 0;
281  overflow = false;
282  for (unsigned int *it : ptr_list) {
283  estimate += *it;
284  if (isOverflow(estimate, resultNBits)) {
285  overflow = true;
286  }
287  }
288  return estimate;
289 }
290 
291 unsigned int LVL1::eFEXtauBDT::multWithOverflow(unsigned int a, unsigned int b,
292  bool &overflow,
293  int resultNBits) {
294  overflow = false;
295  if (b == 0) {
296  return 0;
297  }
298  unsigned int result = a * b;
299  if (a != result / b) {
300  // This shouldn't happen (a and b are in reality 16 bit numbers), but just
301  // in case.
302  overflow = true;
303  }
304 
305  overflow = isOverflow(result, resultNBits);
306 
307  return result;
308 }
309 
310 unsigned int LVL1::eFEXtauBDT::BitLeftShift(unsigned int number, int by,
311  int totalNBits) {
312  if ((number >> (totalNBits - by)) != 0) {
313  return (1 << totalNBits) - 1;
314  }
315  return number << by;
316 }
317 
319  int n_multipliers = sizeof(m_fracMultipliers) / sizeof(m_fracMultipliers[0]);
320 
321  if ((m_eTEstimate >= *m_maxEtThresholdFrac) or m_eTEstimateOverflow or
322  m_HAD_eTEstimateOverflow) {
323 
324  m_fracCondition = (1 << (n_multipliers - 1)) - 1;
325  return;
326  }
327 
328  if (m_EM_eTEstimateOverflow) {
329  m_fracCondition = 0;
330  return;
331  }
332 
333  m_hadEstimateShifted = BitLeftShift(m_HAD_eTEstimate, 3, ENERGY_WIDTH);
334  int i = 0;
335  for (; i < n_multipliers; i++) {
336 
337  bool overflow;
338  m_emEtXMultiplier[i] = multWithOverflow(
339  *(m_fracMultipliers[i]), m_EM_eTEstimate, overflow, ENERGY_WIDTH);
340  m_emEtXMultiplierOverflow[i] = (int)overflow;
341 
342  if ((m_hadEstimateShifted < m_emEtXMultiplier[i]) or
343  m_emEtXMultiplierOverflow[i]) {
344  break;
345  }
346  }
347  m_fracCondition = i;
348 }
349 
351  computeBDTScore();
352  int n_thresholds = sizeof(m_bdtThresholds) / sizeof(m_bdtThresholds[0]);
353 
354  int toShiftRight = m_bdt.getScorePrecision() - PARAM_WIDTH;
355  // Only compare the MSB bits of the BDT score to the thresholds provided in
356  // the parameters
357  m_bdtScoreShifted = (m_bdtScore >> toShiftRight);
358 
359  if ((m_eTEstimate >= *m_maxEtThreshold) or m_eTEstimateOverflow or
360  m_eTEstimate < *m_bdtMinEtThreshold) {
361 
362  m_bdtCondition = (1 << (n_thresholds - 1)) - 1;
363  return;
364  }
365 
366  int i = 0;
367  for (; i < n_thresholds; i++) {
368  if (m_bdtScoreShifted < *(m_bdtThresholds[i])) {
369  break;
370  }
371  }
372  m_bdtCondition = i;
373 }
374 
375 // Check if central tower qualifies as a seed tower for the tau algorithm
376 // Bitwise computation (as opposed to implementation in eFEXtauBDTAlgo.cxx)
378  m_isSeeded = true;
379 
380  // Get central tower ET
381  unsigned int centralET = m_towers[4];
382 
383  // Loop over all cells and check that the central tower is a local maximum
384  for (unsigned int beta = 0; beta < 3; beta++) {
385  for (unsigned int bphi = 0; bphi < 3; bphi++) {
386  int flatIndex = flatTowerIndex(beta, bphi);
387  // Don't need to compare central cell with itself
388  if ((beta == 1) && (bphi == 1)) {
389  continue;
390  }
391 
392  // Cells to the up and right must have strictly lesser ET
393  if (beta == 2 || (beta == 1 && bphi == 2)) {
394  if (centralET <= m_towers[flatIndex]) {
395  m_isSeeded = false;
396  }
397  }
398  // Cells down and to the left must have lesser or equal ET. If strictly
399  // lesser would create zero TOB if two adjacent cells had equal energy
400  else if (beta == 0 || (beta == 1 && bphi == 0)) {
401  if (centralET < m_towers[flatIndex]) {
402  m_isSeeded = false;
403  }
404  }
405  }
406  }
407 
408  if (m_eTEstimate < *m_etThreshold) {
409  m_isSeeded = false;
410  }
411 
412  m_log->msg(MSG::DEBUG) << "Seeded: " << (int)m_isSeeded << endmsg;
413 }
414 
415 unsigned int LVL1::eFEXtauBDT::getET() const {
416  if (m_eTEstimateOverflow) {
417  return (1 << ENERGY_WIDTH) - 1;
418  }
419  if (m_eTEstimate < *m_etThreshold) {
420  return 0;
421  }
422  return m_eTEstimate;
423 }
LVL1::eFEXtauBDT::computeBDTCondition
void computeBDTCondition()
Definition: eFEXtauBDT.cxx:350
LVL1::eFEXtauBDT::m_towers
std::vector< unsigned int > m_towers
Definition: eFEXtauBDT.h:131
beamspotnt.var
var
Definition: bin/beamspotnt.py:1393
LVL1::eFEXtauBDT::m_bdt
eFEXBDT m_bdt
Definition: eFEXtauBDT.h:149
LVL1::eFEXtauBDT::flatTowerIndex
int flatTowerIndex(int eta, int phi)
Definition: eFEXtauBDT.cxx:255
get_generator_info.result
result
Definition: get_generator_info.py:21
LVL1::eFEXtauBDT::setPointerToMaxETParam
void setPointerToMaxETParam(unsigned int *maxEtThreshold)
Definition: eFEXtauBDT.cxx:68
LVL1::eFEXtauBDT::initPointers
void initPointers(const std::vector< std::vector< int >> &scells, std::vector< unsigned int * > &ptr_list)
Definition: eFEXtauBDT.cxx:117
LVL1::eFEXtauBDT::initHADETPointers
void initHADETPointers()
Definition: eFEXtauBDT.cxx:113
LVL1::eFEXtauBDT::setPointerToBDTMinETParam
void setPointerToBDTMinETParam(unsigned int *bdtMinEtThreshold)
Definition: eFEXtauBDT.cxx:78
index
Definition: index.py:1
eFEXBDT::getVariables
const std::vector< BDTVariable > & getVariables() const
Definition: eFEXBDT.h:58
LVL1::eFEXtauBDT::setPointerToMaxETParamFrac
void setPointerToMaxETParamFrac(unsigned int *maxEtThreshold)
Definition: eFEXtauBDT.cxx:73
skel.it
it
Definition: skel.GENtoEVGEN.py:407
LVL1::eFEXtauBDT::getET
unsigned int getET() const
Definition: eFEXtauBDT.cxx:415
LVL1::eFEXtauBDT::debugPrintBDTVariables
void debugPrintBDTVariables()
Definition: eFEXtauBDT.cxx:191
LVL1::eFEXtauBDT::debugPrintTowers
void debugPrintTowers()
Definition: eFEXtauBDT.cxx:243
LVL1::eFEXtauBDT::isOverflow
bool isOverflow(unsigned int number, int nBits)
Definition: eFEXtauBDT.cxx:270
dbg::ptr
void * ptr(T *p)
Definition: SGImplSvc.cxx:74
LVL1::eFEXtauBDT::setPointerToBDTThresholdsParam
void setPointerToBDTThresholdsParam(int index, unsigned int *bdtThresholds)
Definition: eFEXtauBDT.cxx:59
LVL1::eFEXtauBDT::computeETEstimate
void computeETEstimate()
Definition: eFEXtauBDT.cxx:225
LVL1::eFEXtauBDT::computeEMETEstimate
void computeEMETEstimate()
Definition: eFEXtauBDT.cxx:231
ENERGY_WIDTH
#define ENERGY_WIDTH
Definition: eFEXtauBDT.cxx:15
LVL1::eFEXtauBDT::computeTowers
void computeTowers()
Definition: eFEXtauBDT.cxx:257
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
LVL1::eFEXtauBDT::m_log
AthAlgTool * m_log
Definition: eFEXtauBDT.h:150
LVL1::eFEXtauBDT::eFEXtauBDT
eFEXtauBDT(AthAlgTool *log, std::string config_path)
Constructors.
Definition: eFEXtauBDT.cxx:19
lumiFormat.i
int i
Definition: lumiFormat.py:85
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
LVL1::eFEXtauBDT::setPointerToETThresholdParam
void setPointerToETThresholdParam(unsigned int *etThreshold)
Definition: eFEXtauBDT.cxx:64
LVL1::eFEXtauBDT::setPointerToSCell
void setPointerToSCell(int eta, int phi, int layer, unsigned int *sCellPtr)
Definition: eFEXtauBDT.cxx:33
LVL1::eFEXtauBDT::initETPointers
void initETPointers()
Definition: eFEXtauBDT.cxx:105
fitman.by
by
Definition: fitman.py:411
LVL1::eFEXtauBDT::initTowersPointers
void initTowersPointers()
Definition: eFEXtauBDT.cxx:141
internal_poltrig::estimate
REAL estimate(const int &elen, REAL *e)
Definition: PolygonTriangulator.cxx:451
LVL1::eFEXtauBDT::initBDTVars
void initBDTVars()
Definition: eFEXtauBDT.cxx:149
python.selection.number
number
Definition: selection.py:20
LVL1::eFEXtauBDT::m_emEtXMultiplier
std::vector< unsigned int > m_emEtXMultiplier
Definition: eFEXtauBDT.h:129
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
LVL1::eFEXtauBDT::m_fracMultipliers
unsigned int * m_fracMultipliers[3]
Definition: eFEXtauBDT.h:111
eFEXBDT::getNTowers
int getNTowers() const
Definition: eFEXBDT.h:80
LVL1::eFEXtauBDT::buildBDTVariables
void buildBDTVariables()
Definition: eFEXtauBDT.cxx:201
LVL1::eFEXtauBDT::setPointerToFracMultipliersParam
void setPointerToFracMultipliersParam(int index, unsigned int *fracMultipliers)
Definition: eFEXtauBDT.cxx:54
DeMoScan.index
string index
Definition: DeMoScan.py:362
eFEXtauBDT.h
a
TList * a
Definition: liststreamerinfos.cxx:10
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
LVL1::eFEXtauBDT::computeBDTScore
void computeBDTScore()
Definition: eFEXtauBDT.cxx:214
DEBUG
#define DEBUG
Definition: page_access.h:11
AthCommonMsg::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
LVL1::eFEXtauBDT::computeFracCondition
void computeFracCondition()
Definition: eFEXtauBDT.cxx:318
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
PARAM_WIDTH
#define PARAM_WIDTH
Definition: eFEXtauBDT.cxx:16
LVL1::eFEXtauBDT::superCellToPtr
unsigned int * superCellToPtr(int eta, int phi, int layer)
Definition: eFEXtauBDT.cxx:83
LVL1::eFEXtauBDT::m_emEtXMultiplierOverflow
std::vector< unsigned int > m_emEtXMultiplierOverflow
Definition: eFEXtauBDT.h:130
LVL1::eFEXtauBDT::~eFEXtauBDT
virtual ~eFEXtauBDT()
Destructor.
Definition: eFEXtauBDT.cxx:31
AthAlgTool
Definition: AthAlgTool.h:26
LVL1::eFEXtauBDT::BitLeftShift
unsigned int BitLeftShift(unsigned int number, int by, int totalNBits)
Definition: eFEXtauBDT.cxx:310
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
LVL1::eFEXtauBDT::m_bdtVars
std::vector< unsigned int > m_bdtVars
Definition: eFEXtauBDT.h:134
LVL1::eFEXtauBDT::initEMETPointers
void initEMETPointers()
Definition: eFEXtauBDT.cxx:109
LVL1::eFEXtauBDT::computeIsCentralTowerSeed
void computeIsCentralTowerSeed()
Definition: eFEXtauBDT.cxx:377
LVL1::eFEXtauBDT::computeEstimate
unsigned int computeEstimate(std::vector< unsigned int * > &ptr_list, bool &overflow, int resultNBits)
Definition: eFEXtauBDT.cxx:278
LVL1::eFEXtauBDT::computeHADETEstimate
void computeHADETEstimate()
Definition: eFEXtauBDT.cxx:237
LVL1::eFEXtauBDT::multWithOverflow
unsigned int multWithOverflow(unsigned int a, unsigned int b, bool &overflow, int resultNBits)
Definition: eFEXtauBDT.cxx:291
BDTVariable
Definition: eFEXBDT.h:16
LVL1::eFEXtauBDT::next
void next()
Definition: eFEXtauBDT.cxx:180