15 #define ENERGY_WIDTH 16
20 : m_bdt(config_path), m_log(
log) {
34 unsigned int *sCellPtr) {
37 m_em0cells[eta][phi] = sCellPtr;
40 m_em1cells[eta][phi] = sCellPtr;
43 m_em2cells[eta][phi] = sCellPtr;
46 m_em3cells[eta][phi] = sCellPtr;
49 m_hadcells[eta][phi] = sCellPtr;
55 int index,
unsigned int *fracMultiplierPtr) {
56 m_fracMultipliers[
index] = fracMultiplierPtr;
60 int index,
unsigned int *bdtThresholdPtr) {
61 m_bdtThresholds[
index] = bdtThresholdPtr;
65 m_etThreshold = etThreshold;
69 unsigned int *maxEtThreshold) {
70 m_maxEtThreshold = maxEtThreshold;
74 unsigned int *maxEtThreshold) {
75 m_maxEtThresholdFrac = maxEtThreshold;
79 unsigned int *bdtMinEtThreshold) {
80 m_bdtMinEtThreshold = bdtMinEtThreshold;
84 unsigned int *
ptr = 0;
87 ptr = m_em0cells[eta][phi];
90 ptr = m_em1cells[eta][phi];
93 ptr = m_em2cells[eta][phi];
96 ptr = m_em3cells[eta][phi];
99 ptr = m_hadcells[eta][phi];
106 initPointers(m_bdt.getETSCells(), m_eTComputeSCellPointers);
110 initPointers(m_bdt.getEMETSCells(), m_EM_eTComputeSCellPointers);
114 initPointers(m_bdt.getHADETSCells(), m_HAD_eTComputeSCellPointers);
118 std::vector<unsigned int *> &ptr_list) {
120 for (
const auto& scell : scells) {
123 int layer = scell[2];
124 m_log->msg(
MSG::DEBUG) <<
"\teta=" << eta <<
"\tphi=" << phi
126 unsigned int *
ptr = superCellToPtr(eta, phi,
layer);
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(
135 " to a pointer to supercell. Are they within range?");
137 ptr_list.push_back(
ptr);
142 m_towersComputeSCellPointers.resize(m_bdt.getNTowers());
143 for (
size_t i = 0;
i < m_towers.size();
i++) {
145 initPointers(m_bdt.getTowerSCells(
i), m_towersComputeSCellPointers[
i]);
150 for (
size_t i = 0;
i < m_bdt.getVariables().
size();
i++) {
153 m_log->msg(
MSG::DEBUG) <<
i <<
" is " <<
var.m_name <<
", sum of supercells"
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];
160 m_log->msg(
MSG::DEBUG) <<
"\teta=" << eta <<
"\tphi=" << phi
162 unsigned int *
ptr = superCellToPtr(eta, phi,
layer);
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(
171 " to a pointer to supercell. Are they within range?");
173 pointersToSCells.push_back(
ptr);
176 m_bdtVarComputeSCellPointers.push_back(pointersToSCells);
184 computeHADETEstimate();
185 computeEMETEstimate();
186 computeBDTCondition();
187 computeFracCondition();
188 computeIsCentralTowerSeed();
192 std::string bdtVariables =
"";
193 for (
size_t i = 0;
i < m_bdtVars.size();
i++) {
202 for (
size_t i = 0;
i < m_bdtVarComputeSCellPointers.size();
i++) {
204 m_bdtVars[
i] = computeEstimate(m_bdtVarComputeSCellPointers[
i], overflow,
210 debugPrintBDTVariables();
211 m_bdtVarsComputed =
true;
216 if (m_bdtVarsComputed ==
false) {
218 <<
"BDT Variables not computed. BDT score will be garbage." <<
endmsg;
221 m_bdtScore = m_bdt.getBDT().decision_function(m_bdtVars)[0];
226 m_eTEstimate = computeEstimate(m_eTComputeSCellPointers, m_eTEstimateOverflow,
232 m_EM_eTEstimate = computeEstimate(m_EM_eTComputeSCellPointers,
238 m_HAD_eTEstimate = computeEstimate(m_HAD_eTComputeSCellPointers,
240 m_log->msg(
MSG::DEBUG) <<
"HAD ET Estimate: " << m_HAD_eTEstimate <<
endmsg;
245 for (
int eta = 0; eta < 3; eta++) {
246 for (
int phi = 0; phi < 3; phi++) {
247 int flatIndex = flatTowerIndex(eta, phi);
249 <<
"Tower " << flatIndex <<
" ET (eta=" << eta <<
", phi=" << phi
250 <<
"): " << m_towers[flatIndex] <<
endmsg;
258 for (
size_t i = 0;
i < m_towers.size();
i++) {
260 m_towers[
i] = computeEstimate(m_towersComputeSCellPointers[
i], overflow,
271 if ((
number >> nBits) != 0) {
279 bool &overflow,
int resultNBits) {
282 for (
unsigned int *
it : ptr_list) {
284 if (isOverflow(
estimate, resultNBits)) {
305 overflow = isOverflow(
result, resultNBits);
312 if ((
number >> (totalNBits -
by)) != 0) {
313 return (1 << totalNBits) - 1;
319 int n_multipliers =
sizeof(m_fracMultipliers) /
sizeof(m_fracMultipliers[0]);
321 if ((m_eTEstimate >= *m_maxEtThresholdFrac) or m_eTEstimateOverflow or
322 m_HAD_eTEstimateOverflow) {
324 m_fracCondition = (1 << (n_multipliers - 1)) - 1;
328 if (m_EM_eTEstimateOverflow) {
333 m_hadEstimateShifted = BitLeftShift(m_HAD_eTEstimate, 3,
ENERGY_WIDTH);
335 for (;
i < n_multipliers;
i++) {
338 m_emEtXMultiplier[
i] = multWithOverflow(
339 *(m_fracMultipliers[
i]), m_EM_eTEstimate, overflow,
ENERGY_WIDTH);
340 m_emEtXMultiplierOverflow[
i] = (
int)overflow;
342 if ((m_hadEstimateShifted < m_emEtXMultiplier[
i]) or
343 m_emEtXMultiplierOverflow[
i]) {
352 int n_thresholds =
sizeof(m_bdtThresholds) /
sizeof(m_bdtThresholds[0]);
354 int toShiftRight = m_bdt.getScorePrecision() -
PARAM_WIDTH;
357 m_bdtScoreShifted = (m_bdtScore >> toShiftRight);
359 if ((m_eTEstimate >= *m_maxEtThreshold) or m_eTEstimateOverflow or
360 m_eTEstimate < *m_bdtMinEtThreshold) {
362 m_bdtCondition = (1 << (n_thresholds - 1)) - 1;
367 for (;
i < n_thresholds;
i++) {
368 if (m_bdtScoreShifted < *(m_bdtThresholds[
i])) {
381 unsigned int centralET = m_towers[4];
385 for (
unsigned int bphi = 0; bphi < 3; bphi++) {
386 int flatIndex = flatTowerIndex(
beta, bphi);
388 if ((
beta == 1) && (bphi == 1)) {
393 if (
beta == 2 || (
beta == 1 && bphi == 2)) {
394 if (centralET <= m_towers[flatIndex]) {
400 else if (
beta == 0 || (
beta == 1 && bphi == 0)) {
401 if (centralET < m_towers[flatIndex]) {
408 if (m_eTEstimate < *m_etThreshold) {
416 if (m_eTEstimateOverflow) {
419 if (m_eTEstimate < *m_etThreshold) {