ATLAS Offline Software
Loading...
Searching...
No Matches
EFexTauAlgorithm.cxx
Go to the documentation of this file.
1// Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
2
9
10#include "EFexTauAlgorithm.h"
11#include "CaloEvent/CaloCellContainer.h"
12#include <string>
13#include "TH2F.h"
14#include "TVector2.h"
15#include "xAODTrigger/EmTauRoIAuxContainer.h"
16#include "xAODTrigger/EmTauRoI.h"
17
18namespace
19{
20 const static SG::AuxElement::Decorator<float> decR3ClusterET("R3ClusterET");
21 const static SG::AuxElement::Decorator<float> decR3OreClusterET("R3_Ore_ClusterET");
22 const static SG::AuxElement::Decorator<float> decR3BCClusterET("R3_BC_ClusterET");
23 const static SG::AuxElement::Decorator<float> decR3BCClusterIso("R3_BC_ClusterIso");
24 const static SG::AuxElement::Decorator<float> decR3OreClusterIso("R3_Ore_ClusterIso");
25 const static SG::AuxElement::Decorator<bool> decR3OreIsoPass12("R3_Ore_ClusterIso_12pass");
26 const static SG::AuxElement::Decorator<bool> decR3OreIsoPass20("R3_Ore_ClusterIso_20pass");
27 const static SG::AuxElement::Decorator<bool> decR3BCIsoPass12("R3_BC_ClusterIso_12pass");
28 const static SG::AuxElement::Decorator<bool> decR3BCIsoPass20("R3_BC_ClusterIso_20pass");
29 const static SG::AuxElement::Decorator<float> decR3ClusterIso("R3ClusterIso");
30} // namespace
31
32LVL1::EFexTauAlgorithm::EFexTauAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
33 : AthReentrantAlgorithm(name, pSvcLocator)
34{
35 declareProperty("InputSuperCellContainer", m_inputCellContainerKey = "SCell");
36 declareProperty("InputTileCellContainer", m_inputTileCellContainerKey = "AllCalo");
37 declareProperty("InputTriggerTowerContainer", m_inputTriggerTowerContainerKey = "xAODTriggerTowers");
38 declareProperty("OutputClusterName", m_outputClusterName = "SCluster");
39
40 declareProperty("UseTileCells", m_use_tileCells = false, "Use Tile cells instead of TriggerTowers");
41
42 declareProperty("CleanCellContainerSkim", m_useProvenanceSkim = false);
43 declareProperty("CleanCellContainer", m_useProvenance = true);
44 declareProperty("QualBitMask", m_qualBitMask = 0x40);
45
46 declareProperty("NominalDigitizationValue", m_nominalDigitization = 25.);
47 declareProperty("NominalNoiseThreshold", m_nominalNoise_thresh = 100.);
48 declareProperty("TimeThreshold", m_timeThr = 200);
49}
50
51LVL1::EFexTauAlgorithm::~EFexTauAlgorithm() = default;
52
53StatusCode
54LVL1::EFexTauAlgorithm::initialize()
55{
56 ATH_MSG_DEBUG("initialize");
57 ATH_CHECK(m_inputCellContainerKey.initialize());
58 ATH_CHECK(m_inputTileCellContainerKey.initialize(m_use_tileCells));
59 ATH_CHECK(m_inputTriggerTowerContainerKey.initialize(!m_use_tileCells));
60
61 ATH_CHECK(m_outputClusterName.initialize());
62 return StatusCode::SUCCESS;
63}
64
65float LVL1::EFexTauAlgorithm::CaloCellET(const CaloCell *const &inputCell, float digitScale, float digitThreshold) const
66{
67 if (inputCell == nullptr)
68 return 0.;
69 // Check that timing is correct
70 if (m_useProvenance)
71 {
72 bool correctProv = (inputCell->provenance() & m_qualBitMask);
73 if (!correctProv)
74 return 0.;
75 }
76 // Calculates the ET (before digitization)
77 float inputCell_energy = inputCell->energy();
78 float inputCell_eta = inputCell->eta();
79 float inputCell_ET = cosh(inputCell_eta) ? inputCell_energy / cosh(inputCell_eta) : inputCell_energy;
80 // Check to see if negative ET values are allowed
81 bool allowNegs = false;
82 if (digitScale < 0.)
83 {
84 digitScale = std::abs(digitScale);
85 allowNegs = true;
86 }
87 if (inputCell_ET == 0)
88 return 0.;
89 else if (digitScale == 0)
90 return inputCell_ET;
91 if (allowNegs || inputCell_ET > 0.)
92 {
93 // Split up ET into magnitude & whether it's positive or negative
94 float posOrNeg = inputCell_ET ? inputCell_ET / std::abs(inputCell_ET) : inputCell_ET;
95 inputCell_ET = std::abs(inputCell_ET);
96 // If no digitisation, return ET following noise cut
97 if (digitScale == 0)
98 {
99 if (inputCell_ET > digitThreshold)
100 return inputCell_ET * posOrNeg;
101 else
102 return 0.;
103 }
104 // Apply digitization & then noise cut
105 else
106 {
107 float divET = digitScale ? inputCell_ET / digitScale : inputCell_ET;
108 int roundET = divET;
109 float result = digitScale * roundET;
110 if (digitThreshold == 0)
111 return result * posOrNeg;
112 else if (result >= digitThreshold)
113 return result * posOrNeg;
114 else
115 return 0;
116 }
117 }
118 else
119 return 0.;
120}
121
122StatusCode
123LVL1::EFexTauAlgorithm::execute(const EventContext& ctx) const
124{
125
126 // supercells are used by both methods
127 auto scellsHandle = SG::makeHandle(m_inputCellContainerKey, ctx);
128 if (!scellsHandle.isValid())
129 {
130 ATH_MSG_ERROR("Failed to retrieve " << m_inputCellContainerKey.key());
131 return StatusCode::FAILURE;
132 }
133 CaloConstCellContainer scells(SG::VIEW_ELEMENTS);
134 if (m_useProvenanceSkim)
135 {
136 for (const CaloCell *scell : *scellsHandle)
137 if (scell->provenance() & m_qualBitMask)
138 scells.push_back(scell);
139 }
140 else
141 scells.assign(scellsHandle->begin(), scellsHandle->end());
142
143 CaloConstCellContainer tileCellCont(SG::VIEW_ELEMENTS);
144 const xAOD::TriggerTowerContainer *TTs{nullptr};
145 if (m_use_tileCells)
146 {
147 auto tileCellHandle = SG::makeHandle(m_inputTileCellContainerKey, ctx);
148 if (!tileCellHandle.isValid())
149 {
150 ATH_MSG_ERROR("Failed to retrieve " << m_inputTileCellContainerKey.key());
151 return StatusCode::FAILURE;
152 }
153 tileCellCont.assign(tileCellHandle->begin(), tileCellHandle->end());
154 }
155 else
156 {
157 auto triggerTowerHandle = SG::makeHandle(m_inputTriggerTowerContainerKey, ctx);
158 if (!triggerTowerHandle.isValid())
159 {
160 ATH_MSG_ERROR("Failed to retrieve " << m_inputTriggerTowerContainerKey.key());
161 return StatusCode::FAILURE;
162 }
163 TTs = triggerTowerHandle.cptr();
164 }
165 auto clustersForTau = std::make_unique<xAOD::EmTauRoIContainer>();
166 auto auxClustersForTau = std::make_unique<xAOD::EmTauRoIAuxContainer>();
167 clustersForTau->setStore(auxClustersForTau.get());
168
169 // prepare all supercells vector in whole ATLAS detector
170 std::vector<const CaloCell *> allSuperCells(scells.begin(), scells.end());
171
172 // clear all TH2 histograms for supercell map
173 TH2F supercellMapEM0("SupercellMapEM0", "Supercell map of EM0", 98, -4.9, 4.9, 64, 0, 2 * M_PI);
174 TH2F supercellMapEM1("SupercellMapEM1", "Supercell map of EM1", 392, -4.9, 4.9, 64, 0, 2 * M_PI);
175 TH2F supercellMapEM2("SupercellMapEM2", "Supercell map of EM2", 392, -4.9, 4.9, 64, 0, 2 * M_PI);
176 TH2F supercellMapEM1_coarse("SupercellMapEM1_coarse", "Supercell map of EM1 coarse", 196, -4.9, 4.9, 64, 0, 2 * M_PI);
177 TH2F supercellMapEM2_coarse("SupercellMapEM2_coarse", "Supercell map of EM2 coarse", 196, -4.9, 4.9, 64, 0, 2 * M_PI);
178 TH2F supercellMapEM3("SupercellMapEM3", "Supercell map of EM3", 98, -4.9, 4.9, 64, 0, 2 * M_PI);
179 TH2F supercellMapHAD("SupercellMapHAD", "Supercell map of HAD", 98, -4.9, 4.9, 64, 0, 2 * M_PI);
180 TH2F supercellMapTWR("SupercellMapTWR", "Supercell map of TWR", 98, -4.9, 4.9, 64, 0, 2 * M_PI);
181
182 int currentSampling = 0;
183 float currentEta = 0;
184 float currentPhi = 0;
185 float currentCellEt = 0;
186
187 // fill energy in all TH2 histograms for supercell map
188 for (auto scell : allSuperCells)
189 {
190 currentSampling = scell->caloDDE()->getSampling();
191 currentEta = scell->eta();
192 currentPhi = TVector2::Phi_0_2pi(scell->phi());
193 currentCellEt = CaloCellET(scell, m_nominalDigitization, m_nominalNoise_thresh);
194
195 // Store maps per layer
196 if (currentSampling == 0 || currentSampling == 4)
197 {
198 supercellMapEM0.Fill(currentEta, currentPhi, currentCellEt);
199 }
200 else if (currentSampling == 1 || currentSampling == 5)
201 {
202 supercellMapEM1.Fill(currentEta, currentPhi, currentCellEt);
203 supercellMapEM1_coarse.Fill(currentEta, currentPhi, currentCellEt);
204 }
205 else if (currentSampling == 2 || currentSampling == 6)
206 {
207 supercellMapEM2.Fill(currentEta, currentPhi, currentCellEt);
208 supercellMapEM2_coarse.Fill(currentEta, currentPhi, currentCellEt);
209 }
210 else if (currentSampling == 3 || currentSampling == 7)
211 {
212 supercellMapEM3.Fill(currentEta, currentPhi, currentCellEt);
213 }
214 else
215 {
216 supercellMapHAD.Fill(currentEta, currentPhi, currentCellEt);
217 }
218
219 // Store a map sum of all layers
220 supercellMapTWR.Fill(currentEta, currentPhi, currentCellEt);
221 }
222
223 // Need to also loop over Run-I towers to get central region hadronic energy
224 for (const xAOD::TriggerTower *tt : *TTs)
225 {
226
227 // Only use towers within 1.5, and in Tile
228 if (tt->sampling() != 1 || std::abs(tt->eta()) > 1.5)
229 {
230 continue;
231 }
232
233 // Conversion into ET
234 float cpET = tt->cpET() * 500.; // EM energy scale: 1 unit corresponds to 500 MeV
235 if (cpET < 0.)
236 {
237 cpET = 0;
238 }
239
240 // Fill hadronic maps
241 supercellMapHAD.Fill(tt->eta(), TVector2::Phi_0_2pi(tt->phi()), cpET);
242 supercellMapTWR.Fill(tt->eta(), TVector2::Phi_0_2pi(tt->phi()), cpET);
243 }
244
245 // Find local maxima
246 std::vector<TLorentzVector> localMaxima;
247 localMaxima.reserve(200);
248
249 // X is eta, Y is phi
250 for (int i = 0; i < supercellMapTWR.GetNbinsX(); i++)
251 {
252 for (int j = 0; j < supercellMapTWR.GetNbinsY(); j++)
253 {
254 // Start by filtering out 'useless towers' (ie: anything less than a GeV)
255 double towerET = supercellMapTWR.GetBinContent(i + 1, j + 1);
256 if (towerET < 1000.)
257 continue;
258
259 // Check if the current tower has the largest ET in this 3x3 window
260
261 // Need to be careful with wrap-around in phi, unfortunately.
262 std::vector<double> binsAbove;
263 binsAbove.reserve(4);
264 std::vector<double> binsBelow;
265 binsBelow.reserve(4);
266
267 // Handle wrap-around in phi
268 int aboveInPhi = j + 1;
269 if (j == supercellMapTWR.GetNbinsY() - 1)
270 aboveInPhi = 0;
271 int belowInPhi = j - 1;
272 if (j == 0)
273 belowInPhi = supercellMapTWR.GetNbinsY() - 1;
274
275 // The convention here is arbitrary, but needs to be mirrored
276 // Take the cells in the next row in phi, and the one cell above in eta
277 binsAbove.push_back(supercellMapTWR.GetBinContent(i, aboveInPhi + 1));
278 binsAbove.push_back(supercellMapTWR.GetBinContent(i + 1, aboveInPhi + 1));
279 binsAbove.push_back(supercellMapTWR.GetBinContent(i + 2, aboveInPhi + 1));
280 binsAbove.push_back(supercellMapTWR.GetBinContent(i + 2, j + 1));
281
282 // Inversely so for the bins below
283 binsBelow.push_back(supercellMapTWR.GetBinContent(i, belowInPhi + 1));
284 binsBelow.push_back(supercellMapTWR.GetBinContent(i + 1, belowInPhi + 1));
285 binsBelow.push_back(supercellMapTWR.GetBinContent(i + 2, belowInPhi + 1));
286 binsBelow.push_back(supercellMapTWR.GetBinContent(i, j + 1));
287
288 bool isMax = true;
289
290 // Check if it is a local maximum
291 for (unsigned int k = 0; k < binsAbove.size(); k++)
292 {
293 if (towerET < binsAbove[k])
294 isMax = false;
295 }
296 for (unsigned int k = 0; k < binsBelow.size(); k++)
297 {
298 if (towerET <= binsBelow[k])
299 isMax = false;
300 }
301
302 if (isMax)
303 {
304 TLorentzVector myMaximum;
305 myMaximum.SetPtEtaPhiM(towerET, supercellMapTWR.GetXaxis()->GetBinCenter(i + 1), supercellMapTWR.GetYaxis()->GetBinCenter(j + 1), 0);
306 localMaxima.push_back(myMaximum);
307 }
308 }
309 }
310
311 if (msgLvl(MSG::DEBUG))
312 {
313 for (int i = 1; i < supercellMapEM1_coarse.GetNbinsX() + 1; i++)
314 {
315 for (int j = 1; j < supercellMapEM1_coarse.GetNbinsY() + 1; j++)
316 {
317 ATH_MSG_DEBUG("supercellMapEM1_coarse.GetBinContent(" << i << "," << j << ") " << supercellMapEM1_coarse.GetBinContent(i, j));
318 ATH_MSG_DEBUG("supercellMapEM2_coarse.GetBinContent(" << i << "," << j << ") " << supercellMapEM2_coarse.GetBinContent(i, j));
319 ATH_MSG_DEBUG("supercellMapEM2_coarse.GetXaxis()->GetBinCenter(" << i << ") " << supercellMapEM2_coarse.GetXaxis()->GetBinCenter(i));
320 ATH_MSG_DEBUG("supercellMapEM2_coarse.GetYaxis()->GetBinCenter(" << j << ") " << supercellMapEM2_coarse.GetYaxis()->GetBinCenter(j));
321 }
322 }
323 }
324
325 // Now loop over local maxima, decide what to do
326 for (auto myMaximum : localMaxima)
327 {
328 // Check eta bounds
329 if (std::abs(myMaximum.Eta()) > 2.5)
330 {
331 continue;
332 }
333 // Cluster coordinates
334 // Get eta coordinate for coarse layers EM0, EM3, HAD, and TWR
335 int i = supercellMapTWR.GetXaxis()->FindFixBin(myMaximum.Eta());
336 // Careful, ROOT conventions are -pi,pi
337 int j = supercellMapTWR.GetYaxis()->FindFixBin(TVector2::Phi_0_2pi(myMaximum.Phi()));
338 // Get eta coordinate for fine layers EM1 and EM2
339 int i_fine_start = ((i - 1) * 4) + 1;
340 int i_offset = 0;
341 double i_max = 0;
342 for (unsigned int i_off_cand = 0; i_off_cand < 4; i_off_cand++)
343 {
344 int i_et = supercellMapEM2.GetBinContent(i_fine_start + i_off_cand, j);
345 if (i_et > i_max)
346 {
347 i_max = i_et;
348 i_offset = i_off_cand;
349 }
350 }
351 int i_fine = i_fine_start + i_offset;
352
353 // Prepare Phi Wrap-around
354 int aboveInPhi = j + 1;
355 if (j == supercellMapTWR.GetNbinsY())
356 {
357 aboveInPhi = 1;
358 }
359 int belowInPhi = j - 1;
360 if (j == 1)
361 {
362 belowInPhi = supercellMapTWR.GetNbinsY();
363 }
364
365 // Start calculating total energy
366 // Use Fixed 2x2 cluster, 4 possibilities
367 std::vector<double> allET;
368 allET.reserve(4);
369
370 double ET;
371 // Up and right
372 ET = supercellMapTWR.GetBinContent(i, j);
373 ET += supercellMapTWR.GetBinContent(i + 1, j);
374 ET += supercellMapTWR.GetBinContent(i, aboveInPhi);
375 ET += supercellMapTWR.GetBinContent(i + 1, aboveInPhi);
376 allET.push_back(ET);
377
378 // Up and left
379 ET = supercellMapTWR.GetBinContent(i, j);
380 ET += supercellMapTWR.GetBinContent(i - 1, j);
381 ET += supercellMapTWR.GetBinContent(i, aboveInPhi);
382 ET += supercellMapTWR.GetBinContent(i - 1, aboveInPhi);
383 allET.push_back(ET);
384
385 // Down and left
386 ET = supercellMapTWR.GetBinContent(i, j);
387 ET += supercellMapTWR.GetBinContent(i - 1, j);
388 ET += supercellMapTWR.GetBinContent(i, belowInPhi);
389 ET += supercellMapTWR.GetBinContent(i - 1, belowInPhi);
390 allET.push_back(ET);
391
392 // Down and right
393 ET = supercellMapTWR.GetBinContent(i, j);
394 ET += supercellMapTWR.GetBinContent(i + 1, j);
395 ET += supercellMapTWR.GetBinContent(i, belowInPhi);
396 ET += supercellMapTWR.GetBinContent(i + 1, belowInPhi);
397 allET.push_back(ET);
398
399 // Pick largest resulting sum
400 double eFEXOldCluster = 0;
401 for (unsigned int k = 0; k < allET.size(); k++)
402 {
403 if (allET.at(k) > eFEXOldCluster)
404 eFEXOldCluster = allET.at(k);
405 }
406
407 // Calculate Oregon algorithm reconstructed ET
408 // Determine if the Oregon shapes, which are asymmetric in phi, should be oriented up or down
409 bool sumAboveInPhi = true;
410
411 // Sum the cells above in phi for EM1 and EM2
412 double abovePhiCellET = supercellMapEM2.GetBinContent(i_fine, aboveInPhi);
413
414 // Sum the cells below in phi for EM1 and EM2
415 double belowPhiCellET = supercellMapEM2.GetBinContent(i_fine, belowInPhi);
416
417 // If more energy deposited below in phi, then orient shapes downward
418 if (belowPhiCellET > abovePhiCellET)
419 sumAboveInPhi = false;
420
421 // Hold the coordinate of the non-central phi row
422 int offPhiCoordinate = sumAboveInPhi ? aboveInPhi : belowInPhi;
423
424 // Construct 3x2 EM0 Oregon layer energy
425 double em0OregonET = 0;
426 em0OregonET += supercellMapEM0.GetBinContent(i, j);
427 em0OregonET += supercellMapEM0.GetBinContent(i - 1, j);
428 em0OregonET += supercellMapEM0.GetBinContent(i + 1, j);
429 em0OregonET += supercellMapEM0.GetBinContent(i, offPhiCoordinate);
430 em0OregonET += supercellMapEM0.GetBinContent(i - 1, offPhiCoordinate);
431 em0OregonET += supercellMapEM0.GetBinContent(i + 1, offPhiCoordinate);
432
433 // Construct 5x2 EM1 Oregon layer energy
434 double em1OregonET = 0;
435 em1OregonET += supercellMapEM1.GetBinContent(i_fine, j);
436 em1OregonET += supercellMapEM1.GetBinContent(i_fine - 1, j);
437 em1OregonET += supercellMapEM1.GetBinContent(i_fine - 2, j);
438 em1OregonET += supercellMapEM1.GetBinContent(i_fine + 1, j);
439 em1OregonET += supercellMapEM1.GetBinContent(i_fine + 2, j);
440 em1OregonET += supercellMapEM1.GetBinContent(i_fine, offPhiCoordinate);
441 em1OregonET += supercellMapEM1.GetBinContent(i_fine - 1, offPhiCoordinate);
442 em1OregonET += supercellMapEM1.GetBinContent(i_fine - 2, offPhiCoordinate);
443 em1OregonET += supercellMapEM1.GetBinContent(i_fine + 1, offPhiCoordinate);
444 em1OregonET += supercellMapEM1.GetBinContent(i_fine + 2, offPhiCoordinate);
445
446 // Construct 5x2 EM2 Oregon layer energy
447 double em2OregonET = 0;
448 em2OregonET += supercellMapEM2.GetBinContent(i_fine, j);
449 em2OregonET += supercellMapEM2.GetBinContent(i_fine - 1, j);
450 em2OregonET += supercellMapEM2.GetBinContent(i_fine - 2, j);
451 em2OregonET += supercellMapEM2.GetBinContent(i_fine + 1, j);
452 em2OregonET += supercellMapEM2.GetBinContent(i_fine + 2, j);
453 em2OregonET += supercellMapEM2.GetBinContent(i_fine, offPhiCoordinate);
454 em2OregonET += supercellMapEM2.GetBinContent(i_fine - 1, offPhiCoordinate);
455 em2OregonET += supercellMapEM2.GetBinContent(i_fine - 2, offPhiCoordinate);
456 em2OregonET += supercellMapEM2.GetBinContent(i_fine + 1, offPhiCoordinate);
457 em2OregonET += supercellMapEM2.GetBinContent(i_fine + 2, offPhiCoordinate);
458
459 // Construct 3x2 EM3 Oregon layer energy
460 double em3OregonET = 0;
461 em3OregonET += supercellMapEM3.GetBinContent(i, j);
462 em3OregonET += supercellMapEM3.GetBinContent(i - 1, j);
463 em3OregonET += supercellMapEM3.GetBinContent(i + 1, j);
464 em3OregonET += supercellMapEM3.GetBinContent(i, offPhiCoordinate);
465 em3OregonET += supercellMapEM3.GetBinContent(i - 1, offPhiCoordinate);
466 em3OregonET += supercellMapEM3.GetBinContent(i + 1, offPhiCoordinate);
467
468 // Construct 3x2 HAD Oregon layer energy
469 double hadOregonET = 0;
470 hadOregonET += supercellMapHAD.GetBinContent(i, j);
471 hadOregonET += supercellMapHAD.GetBinContent(i - 1, j);
472 hadOregonET += supercellMapHAD.GetBinContent(i + 1, j);
473 hadOregonET += supercellMapHAD.GetBinContent(i, offPhiCoordinate);
474 hadOregonET += supercellMapHAD.GetBinContent(i - 1, offPhiCoordinate);
475 hadOregonET += supercellMapHAD.GetBinContent(i + 1, offPhiCoordinate);
476
477 // Add individual layer energies to get the full Oregon reconstructed energy
478 double eFEX_OregonET = em0OregonET + em1OregonET + em2OregonET + em3OregonET + hadOregonET;
479
480 // Calculate Oregon algorithm isolation value
481 // Construct inner 3x2 energy
482 double oreIsoInnerET = 0;
483 oreIsoInnerET += supercellMapEM2.GetBinContent(i_fine, j);
484 oreIsoInnerET += supercellMapEM2.GetBinContent(i_fine - 1, j);
485 oreIsoInnerET += supercellMapEM2.GetBinContent(i_fine + 1, j);
486 oreIsoInnerET += supercellMapEM2.GetBinContent(i_fine, offPhiCoordinate);
487 oreIsoInnerET += supercellMapEM2.GetBinContent(i_fine - 1, offPhiCoordinate);
488 oreIsoInnerET += supercellMapEM2.GetBinContent(i_fine + 1, offPhiCoordinate);
489
490 // Construct outer 9x2 energy using inner sum plus extra cells
491 double oreIsoOuterET = oreIsoInnerET;
492 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine - 2, j);
493 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine - 3, j);
494 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine - 4, j);
495 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine + 2, j);
496 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine + 3, j);
497 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine + 4, j);
498 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine - 2, offPhiCoordinate);
499 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine - 3, offPhiCoordinate);
500 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine - 4, offPhiCoordinate);
501 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine + 2, offPhiCoordinate);
502 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine + 3, offPhiCoordinate);
503 oreIsoOuterET += supercellMapEM2.GetBinContent(i_fine + 4, offPhiCoordinate);
504
505 // Calculate isolation value as the ratio of inner over outer energies
506 double eFEX_OregonIso = oreIsoOuterET ? oreIsoInnerET / oreIsoOuterET : oreIsoInnerET;
507
508 // Calculation of isolation cut values discussed here
509 // https://indico.cern.ch/event/867020/contributions/3726146/attachments/2003208/3344698/L1CALOJoint03132020.pdf
510
511 // Set boolean for whether event passes 12 GeV isolation cut, hardcoding isolation threshold for now
512 bool eFEX_OregonIso_12pass = true;
513 if (10000. < eFEX_OregonET && 15000. > eFEX_OregonET && eFEX_OregonIso < 0.69)
514 {
515 eFEX_OregonIso_12pass = false;
516 }
517
518 // Set boolean for whether event passes 20 GeV isolation cut, hardcoding isolation threshold for now
519 bool eFEX_OregonIso_20pass = true;
520 if (20000. < eFEX_OregonET && 25000. > eFEX_OregonET && eFEX_OregonIso < 0.61)
521 {
522 eFEX_OregonIso_20pass = false;
523 }
524
525 // Code Oregon cluster later
526 // Will need to write a library of shape summation functions
527 // EM0: 1x2, 4 possibilities
528 // EM1: 5x2
529 // EM2: 5x2
530 // EM3: 3x2
531 // HAD: 3x2
532
533 //Code TLV bigCluster algorithm
534
535 //make vectors with EM0/3/HAD energies in each of 3x3 bins
536 //Make a function since this is disgusting but vim does it very fast
537 std::vector<double> E_EM0;
538 E_EM0.reserve(9);
539 E_EM0.push_back(supercellMapEM0.GetBinContent(i, j));
540 E_EM0.push_back(supercellMapEM0.GetBinContent(i - 1, j));
541 E_EM0.push_back(supercellMapEM0.GetBinContent(i - 1, aboveInPhi));
542 E_EM0.push_back(supercellMapEM0.GetBinContent(i, aboveInPhi));
543 E_EM0.push_back(supercellMapEM0.GetBinContent(i + 1, aboveInPhi));
544 E_EM0.push_back(supercellMapEM0.GetBinContent(i + 1, j));
545 E_EM0.push_back(supercellMapEM0.GetBinContent(i + 1, belowInPhi));
546 E_EM0.push_back(supercellMapEM0.GetBinContent(i, belowInPhi));
547 E_EM0.push_back(supercellMapEM0.GetBinContent(i - 1, belowInPhi));
548 std::vector<double> E_EM3;
549 E_EM3.reserve(9);
550 E_EM3.push_back(supercellMapEM3.GetBinContent(i, j));
551 E_EM3.push_back(supercellMapEM3.GetBinContent(i - 1, j));
552 E_EM3.push_back(supercellMapEM3.GetBinContent(i - 1, aboveInPhi));
553 E_EM3.push_back(supercellMapEM3.GetBinContent(i, aboveInPhi));
554 E_EM3.push_back(supercellMapEM3.GetBinContent(i + 1, aboveInPhi));
555 E_EM3.push_back(supercellMapEM3.GetBinContent(i + 1, j));
556 E_EM3.push_back(supercellMapEM3.GetBinContent(i + 1, belowInPhi));
557 E_EM3.push_back(supercellMapEM3.GetBinContent(i, belowInPhi));
558 E_EM3.push_back(supercellMapEM3.GetBinContent(i - 1, belowInPhi));
559 std::vector<double> E_HAD;
560 E_HAD.reserve(9);
561 E_HAD.push_back(supercellMapHAD.GetBinContent(i, j));
562 E_HAD.push_back(supercellMapHAD.GetBinContent(i - 1, j));
563 E_HAD.push_back(supercellMapHAD.GetBinContent(i - 1, aboveInPhi));
564 E_HAD.push_back(supercellMapHAD.GetBinContent(i, aboveInPhi));
565 E_HAD.push_back(supercellMapHAD.GetBinContent(i + 1, aboveInPhi));
566 E_HAD.push_back(supercellMapHAD.GetBinContent(i + 1, j));
567 E_HAD.push_back(supercellMapHAD.GetBinContent(i + 1, belowInPhi));
568 E_HAD.push_back(supercellMapHAD.GetBinContent(i, belowInPhi));
569 E_HAD.push_back(supercellMapHAD.GetBinContent(i - 1, belowInPhi));
570 std::vector<double> E_EM12_central;
571 E_EM12_central.reserve(5);
572 std::vector<double> E_EM12_above;
573 E_EM12_above.reserve(6);
574 std::vector<double> E_EM12_below;
575 E_EM12_below.reserve(6);
576 float seedEta = supercellMapTWR.GetXaxis()->GetBinCenter(i);
577 int EM2seedBin = supercellMapEM2_coarse.GetXaxis()->FindBin(seedEta - 0.025);
578
579 // Make a vector with the 5 possible energies in the central phi row
580 E_EM12_central.push_back(supercellMapEM2_coarse.GetBinContent(EM2seedBin - 2, j) + supercellMapEM2_coarse.GetBinContent(EM2seedBin - 1, j) + supercellMapEM1_coarse.GetBinContent(EM2seedBin - 2, j) + supercellMapEM1_coarse.GetBinContent(EM2seedBin - 1, j));
581 E_EM12_central.push_back(supercellMapEM2_coarse.GetBinContent(EM2seedBin - 1, j) + supercellMapEM2_coarse.GetBinContent(EM2seedBin - 0, j) + supercellMapEM1_coarse.GetBinContent(EM2seedBin - 1, j) + supercellMapEM1_coarse.GetBinContent(EM2seedBin - 0, j));
582 E_EM12_central.push_back(supercellMapEM2_coarse.GetBinContent(EM2seedBin - 0, j) + supercellMapEM2_coarse.GetBinContent(EM2seedBin + 1, j) + supercellMapEM1_coarse.GetBinContent(EM2seedBin - 0, j) + supercellMapEM1_coarse.GetBinContent(EM2seedBin + 1, j));
583 E_EM12_central.push_back(supercellMapEM2_coarse.GetBinContent(EM2seedBin + 1, j) + supercellMapEM2_coarse.GetBinContent(EM2seedBin + 2, j) + supercellMapEM1_coarse.GetBinContent(EM2seedBin + 1, j) + supercellMapEM1_coarse.GetBinContent(EM2seedBin + 2, j));
584 E_EM12_central.push_back(supercellMapEM2_coarse.GetBinContent(EM2seedBin + 2, j) + supercellMapEM2_coarse.GetBinContent(EM2seedBin + 3, j) + supercellMapEM1_coarse.GetBinContent(EM2seedBin + 2, j) + supercellMapEM1_coarse.GetBinContent(EM2seedBin + 3, j));
585 //For upper and lower phi rows take the bin with highest energy
586 for (int k = -2; k <= 3; k++)
587 {
588 E_EM12_above.push_back(supercellMapEM2_coarse.GetBinContent(EM2seedBin + k, aboveInPhi) + supercellMapEM1_coarse.GetBinContent(EM2seedBin + k, aboveInPhi));
589 E_EM12_below.push_back(supercellMapEM2_coarse.GetBinContent(EM2seedBin + k, belowInPhi) + supercellMapEM1_coarse.GetBinContent(EM2seedBin + k, belowInPhi));
590 }
591
592 sort(E_EM0.begin(), E_EM0.end());
593 sort(E_EM3.begin(), E_EM3.end());
594 sort(E_HAD.begin(), E_HAD.end());
595 sort(E_EM12_central.begin(), E_EM12_central.end());
596 sort(E_EM12_above.begin(), E_EM12_above.end());
597 sort(E_EM12_below.begin(), E_EM12_below.end());
598
599 double eFEX_BC = 0;
600
601 //Three hottest bins in EM0
602 eFEX_BC = E_EM0.at(E_EM0.size() - 1) + E_EM0.at(E_EM0.size() - 2) + E_EM0.at(E_EM0.size() - 3);
603 //Two hottest bins in EM3
604 eFEX_BC += E_EM3.at(E_EM3.size() - 1) + E_EM3.at(E_EM3.size() - 2);
605 //Three hottest bins in HAD
606 eFEX_BC += E_HAD.at(E_HAD.size() - 1) + E_HAD.at(E_HAD.size() - 2) + E_HAD.at(E_HAD.size() - 3);
607 eFEX_BC += E_EM12_central.at(E_EM12_central.size() - 1);
608 eFEX_BC += E_EM12_above.at(E_EM12_above.size() - 1);
609 eFEX_BC += E_EM12_below.at(E_EM12_below.size() - 1);
610
611 //Bigcluster isolation
612 float nomeFEX_BCiso = 0;
613 nomeFEX_BCiso += E_EM12_central.at(E_EM12_central.size() - 1);
614 nomeFEX_BCiso += E_EM12_above.at(E_EM12_above.size() - 1);
615 nomeFEX_BCiso += E_EM12_below.at(E_EM12_below.size() - 1);
616 nomeFEX_BCiso += E_EM0.at(E_EM0.size() - 1);
617 nomeFEX_BCiso += E_EM0.at(E_EM0.size() - 2);
618 nomeFEX_BCiso += E_EM0.at(E_EM0.size() - 3);
619 nomeFEX_BCiso += E_EM3.at(E_EM3.size() - 1);
620 nomeFEX_BCiso += E_EM3.at(E_EM3.size() - 2);
621
622 float denBCiso = 0;
623 denBCiso += supercellMapTWR.GetBinContent(i, j);
624 denBCiso += supercellMapTWR.GetBinContent(i - 1, j);
625 denBCiso += supercellMapTWR.GetBinContent(i + 1, j);
626 denBCiso += supercellMapTWR.GetBinContent(i, j + 1);
627 denBCiso += supercellMapTWR.GetBinContent(i - 1, j + 1);
628 denBCiso += supercellMapTWR.GetBinContent(i + 1, j + 1);
629 denBCiso += supercellMapTWR.GetBinContent(i, j - 1);
630 denBCiso += supercellMapTWR.GetBinContent(i - 1, j - 1);
631 denBCiso += supercellMapTWR.GetBinContent(i + 1, j - 1);
632
633 float eFEX_BCiso = denBCiso ? nomeFEX_BCiso / denBCiso : nomeFEX_BCiso;
634
635 // Set boolean for whether event passes 12 GeV isolation cut, hardcoding isolation threshold for now
636 bool eFEX_BCiso_12pass = true;
637 if (10000. < eFEX_BC && 15000. > eFEX_BC && eFEX_BCiso < 0.38)
638 {
639 eFEX_BCiso_12pass = false;
640 }
641
642 // Set boolean for whether event passes 20 GeV isolation cut, hardcoding isolation threshold for now
643 bool eFEX_BCiso_20pass = true;
644 if (20000. < eFEX_BC && 25000. > eFEX_BC && eFEX_BCiso < 0.18)
645 {
646 eFEX_BCiso_20pass = false;
647 }
648
649 // Calculate an EM2-based isolation
650 // Center in EM2, offset to the right in eta:
651 int em2i = supercellMapEM2.GetXaxis()->FindFixBin(myMaximum.Eta() + 0.05);
652 // Careful, ROOT conventions are -pi,pi
653 int em2j = supercellMapEM2.GetYaxis()->FindFixBin(TVector2::Phi_0_2pi(myMaximum.Phi()));
654
655 float maximumET = supercellMapEM2.GetBinContent(em2i, em2j);
656 int maximumCell = em2i;
657
658 // Find maximum in EM2
659 for (int k = -2; k < 2; k++)
660 {
661 float ETvalue = supercellMapEM2.GetBinContent(em2i + k, em2j);
662 if (ETvalue > maximumET)
663 {
664 maximumET = ETvalue;
665 maximumCell = em2i + k;
666 }
667 }
668
669 // Find highest pT nearest cell, but we don't care which it is
670 // Note that we don't need to worry about phi wrap-around
671 float nextET = supercellMapEM2.GetBinContent(maximumCell + 1, em2j);
672 nextET = ((supercellMapEM2.GetBinContent(maximumCell - 1, em2j) > nextET) ? supercellMapEM2.GetBinContent(maximumCell - 1, em2j) : nextET);
673 nextET = ((supercellMapEM2.GetBinContent(maximumCell, em2j + 1) > nextET) ? supercellMapEM2.GetBinContent(maximumCell, em2j + 1) : nextET);
674 nextET = ((supercellMapEM2.GetBinContent(maximumCell, em2j - 1) > nextET) ? supercellMapEM2.GetBinContent(maximumCell, em2j - 1) : nextET);
675
676 float numerator = maximumET + nextET;
677
678 // And now, run the full sum: avoid phi wrapping by converting back and forth
679 float denominator = 0;
680 float phicenter = supercellMapEM2.GetYaxis()->GetBinCenter(em2j);
681 for (int eta = -6; eta < 6; eta++)
682 for (int phi = -1; phi < 2; phi++)
683 denominator += supercellMapEM2.GetBinContent(em2i + eta, supercellMapEM2.GetYaxis()->FindFixBin(TVector2::Phi_0_2pi(phicenter + (phi * 0.1))));
684
685 // build cluster for tau
686 //xAOD::EmTauRoI_v2* clForTau = new xAOD::EmTauRoI_v2();
687 xAOD::EmTauRoI *clForTau = new xAOD::EmTauRoI();
688 clustersForTau->push_back(clForTau);
689 clForTau->setEta(myMaximum.Eta());
690 clForTau->setPhi(myMaximum.Phi());
691
692 decR3ClusterET(*clForTau) = eFEXOldCluster;
693 decR3OreClusterET(*clForTau) = eFEX_OregonET;
694 decR3BCClusterET(*clForTau) = eFEX_BC;
695 decR3BCClusterIso(*clForTau) = denBCiso > 0 ? eFEX_BCiso : 0;
696 decR3OreClusterIso(*clForTau) = oreIsoOuterET > 0 ? eFEX_OregonIso : 0;
697
698 decR3OreIsoPass12(*clForTau) = eFEX_OregonIso_12pass;
699 decR3OreIsoPass20(*clForTau) = eFEX_OregonIso_20pass;
700 decR3BCIsoPass12(*clForTau) = eFEX_BCiso_12pass;
701 decR3BCIsoPass20(*clForTau) = eFEX_BCiso_20pass;
702 decR3ClusterIso(*clForTau) = denominator > 0 ? numerator / denominator : 0;
703 }
704
705 auto writeHandle = SG::makeHandle(m_outputClusterName, ctx);
706 ATH_CHECK(writeHandle.record(std::move(clustersForTau), std::move(auxClustersForTau)));
707 return StatusCode::SUCCESS;
708}
M_PI
#define M_PI
Definition ActiveFraction.h:11
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition DVL_algorithms.h:554
AthCommonReentrantAlgorithm< Gaudi::Algorithm >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition AthCommonDataStore.h:145
AthCommonReentrantAlgorithm< Gaudi::Algorithm >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition AthCommonMsg.h:30
AthReentrantAlgorithm
An algorithm that can be simultaneously executed in multiple threads.
Definition AthReentrantAlgorithm.h:74
CaloCell
Data object for each calorimeter readout cell.
Definition CaloCell.h:57
CaloCell::energy
double energy() const
get energy (data member)
Definition CaloCell.h:327
CaloCell::provenance
uint16_t provenance() const
get provenance (data member)
Definition CaloCell.h:354
CaloCell::eta
virtual double eta() const override final
get eta (through CaloDetDescrElement)
Definition CaloCell.h:382
CaloConstCellContainer
CaloCellContainer that can accept const cell pointers.
Definition CaloConstCellContainer.h:45
ConstDataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
ConstDataVector::begin
iterator begin() noexcept
Return an iterator pointing at the beginning of the collection.
ConstDataVector::end
iterator end() noexcept
Return an iterator pointing past the end of the collection.
ConstDataVector::assign
void assign(InputIterator first, InputIterator last)
Assign from iterators.
LVL1::EFexTauAlgorithm::m_useProvenanceSkim
bool m_useProvenanceSkim
clear up container from bad BC by making a new container (Denis, old way)
Definition EFexTauAlgorithm.h:57
LVL1::EFexTauAlgorithm::m_useProvenance
bool m_useProvenance
clear up container from bad BC by skipping scells
Definition EFexTauAlgorithm.h:58
LVL1::EFexTauAlgorithm::initialize
StatusCode initialize() override
Definition EFexTauAlgorithm.cxx:54
LVL1::EFexTauAlgorithm::m_qualBitMask
int m_qualBitMask
Configurable quality bitmask.
Definition EFexTauAlgorithm.h:59
LVL1::EFexTauAlgorithm::m_use_tileCells
bool m_use_tileCells
properties
Definition EFexTauAlgorithm.h:55
LVL1::EFexTauAlgorithm::EFexTauAlgorithm
EFexTauAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Definition EFexTauAlgorithm.cxx:32
LVL1::EFexTauAlgorithm::m_nominalDigitization
float m_nominalDigitization
value of nominal digitisation
Definition EFexTauAlgorithm.h:61
LVL1::EFexTauAlgorithm::m_inputCellContainerKey
SG::ReadHandleKey< CaloCellContainer > m_inputCellContainerKey
input / output
Definition EFexTauAlgorithm.h:47
LVL1::EFexTauAlgorithm::m_inputTileCellContainerKey
SG::ReadHandleKey< CaloCellContainer > m_inputTileCellContainerKey
Tile cell input container.
Definition EFexTauAlgorithm.h:48
LVL1::EFexTauAlgorithm::execute
StatusCode execute(const EventContext &ctx) const override
Definition EFexTauAlgorithm.cxx:123
LVL1::EFexTauAlgorithm::m_outputClusterName
SG::WriteHandleKey< xAOD::EmTauRoIContainer > m_outputClusterName
Definition EFexTauAlgorithm.h:50
LVL1::EFexTauAlgorithm::m_nominalNoise_thresh
float m_nominalNoise_thresh
noise threshold
Definition EFexTauAlgorithm.h:62
LVL1::EFexTauAlgorithm::m_inputTriggerTowerContainerKey
SG::ReadHandleKey< xAOD::TriggerTowerContainer > m_inputTriggerTowerContainerKey
TriggerTowers (if needed)
Definition EFexTauAlgorithm.h:49
LVL1::EFexTauAlgorithm::CaloCellET
float CaloCellET(const CaloCell *const &inputCell, float digitScale, float digitThreshold) const
member functions
Definition EFexTauAlgorithm.cxx:65
SG::AuxElement::Decorator
SG::Decorator< T, ALLOC > Decorator
Definition AuxElement.h:575
xAOD::EmTauRoI_v2::setEta
void setEta(float v)
Set the pseudorapidity of the em/tau candidate.
xAOD::EmTauRoI_v2::setPhi
void setPhi(float v)
Set the azimuthal angle of the em/tau candidate.
SG::VIEW_ELEMENTS
@ VIEW_ELEMENTS
this data object is a view, it does not own its elmts
Definition OwnershipPolicy.h:18
SG::makeHandle
SG::ReadCondHandle< T > makeHandle(const SG::ReadCondHandleKey< T > &key, const EventContext &ctx=Gaudi::Hive::currentContext())
Definition ReadCondHandle.h:270
xAOD::TriggerTowerContainer
TriggerTowerContainer_v2 TriggerTowerContainer
Define the latest version of the TriggerTower container.
Definition TriggerTowerContainer.h:15
xAOD::TriggerTower
TriggerTower_v2 TriggerTower
Define the latest version of the TriggerTower class.
Definition Event/xAOD/xAODTrigL1Calo/xAODTrigL1Calo/TriggerTower.h:16
xAOD::EmTauRoI
EmTauRoI_v2 EmTauRoI
Definition EmTauRoI.h:16
result
Generated on for ATLAS Offline Software by doxygen 1.14.0