ATLAS Offline Software
EFexTauAlgorithm.cxx
Go to the documentation of this file.
1 // Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
2 
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 
18 namespace
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 
32 LVL1::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 
51 LVL1::EFexTauAlgorithm::~EFexTauAlgorithm() = default;
52 
53 StatusCode
54 LVL1::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 
65 float 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 
122 StatusCode
123 LVL1::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 }
xAOD::EmTauRoI_v2::setEta
void setEta(float v)
Set the pseudorapidity of the em/tau candidate.
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
get_generator_info.result
result
Definition: get_generator_info.py:21
LVL1::EFexTauAlgorithm::m_inputTriggerTowerContainerKey
SG::ReadHandleKey< xAOD::TriggerTowerContainer > m_inputTriggerTowerContainerKey
TriggerTowers (if needed)
Definition: EFexTauAlgorithm.h:49
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
LVL1::EFexTauAlgorithm::m_nominalNoise_thresh
float m_nominalNoise_thresh
noise threshold
Definition: EFexTauAlgorithm.h:62
EmTauRoI.h
SG::VIEW_ELEMENTS
@ VIEW_ELEMENTS
this data object is a view, it does not own its elmts
Definition: OwnershipPolicy.h:18
TH2F
Definition: rootspy.cxx:420
ConstDataVector::end
iterator end() noexcept
Return an iterator pointing past the end of the collection.
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
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
EFexTauAlgorithm.h
AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
LVL1::EFexTauAlgorithm::m_useProvenance
bool m_useProvenance
clear up container from bad BC by skipping scells
Definition: EFexTauAlgorithm.h:58
M_PI
#define M_PI
Definition: ActiveFraction.h:11
CaloCell::provenance
uint16_t provenance() const
get provenance (data member)
Definition: CaloCell.h:338
CaloCell::energy
double energy() const
get energy (data member)
Definition: CaloCell.h:311
AthReentrantAlgorithm
An algorithm that can be simultaneously executed in multiple threads.
Definition: AthReentrantAlgorithm.h:83
SG::makeHandle
SG::ReadCondHandle< T > makeHandle(const SG::ReadCondHandleKey< T > &key, const EventContext &ctx=Gaudi::Hive::currentContext())
Definition: ReadCondHandle.h:270
xAOD::EmTauRoI
EmTauRoI_v2 EmTauRoI
Definition: EmTauRoI.h:16
LVL1::EFexTauAlgorithm::m_inputTileCellContainerKey
SG::ReadHandleKey< CaloCellContainer > m_inputTileCellContainerKey
Tile cell input container.
Definition: EFexTauAlgorithm.h:48
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
SG::Decorator
Helper class to provide type-safe access to aux data.
Definition: Decorator.h:58
lumiFormat.i
int i
Definition: lumiFormat.py:92
xAOD::EmTauRoI_v2
Class describing a LVL1 em/tau region of interest.
Definition: EmTauRoI_v2.h:35
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
LVL1::EFexTauAlgorithm::m_inputCellContainerKey
SG::ReadHandleKey< CaloCellContainer > m_inputCellContainerKey
input / output
Definition: EFexTauAlgorithm.h:47
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
xAOD::TriggerTower_v2
Description of TriggerTower_v2.
Definition: TriggerTower_v2.h:49
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
ReadTripsProbsFromCool.denominator
denominator
Definition: ReadTripsProbsFromCool.py:96
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
TH2F::GetBinContent
double GetBinContent(int) const
Definition: rootspy.cxx:425
TrigConf::name
Definition: HLTChainList.h:35
LVL1::EFexTauAlgorithm::m_outputClusterName
SG::WriteHandleKey< xAOD::EmTauRoIContainer > m_outputClusterName
Definition: EFexTauAlgorithm.h:50
LVL1::EFexTauAlgorithm::initialize
StatusCode initialize() override
Definition: EFexTauAlgorithm.cxx:54
LVL1::EFexTauAlgorithm::execute
StatusCode execute(const EventContext &ctx) const override
Definition: EFexTauAlgorithm.cxx:123
ConstDataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
CaloCellContainer.h
CaloCell
Data object for each calorimeter readout cell.
Definition: CaloCell.h:57
CaloConstCellContainer
CaloCellContainer that can accept const cell pointers.
Definition: CaloConstCellContainer.h:45
LVL1::EFexTauAlgorithm::m_nominalDigitization
float m_nominalDigitization
value of nominal digitisation
Definition: EFexTauAlgorithm.h:61
DEBUG
#define DEBUG
Definition: page_access.h:11
LVL1::EFexTauAlgorithm::CaloCellET
float CaloCellET(const CaloCell *const &inputCell, float digitScale, float digitThreshold) const
member functions
Definition: EFexTauAlgorithm.cxx:65
EmTauRoIAuxContainer.h
xAOD::EmTauRoI_v2::setPhi
void setPhi(float v)
Set the azimuthal angle of the em/tau candidate.
LVL1::EFexTauAlgorithm::m_qualBitMask
int m_qualBitMask
Configurable quality bitmask.
Definition: EFexTauAlgorithm.h:59
TileDCSDataPlotter.tt
tt
Definition: TileDCSDataPlotter.py:874
ConstDataVector::begin
iterator begin() noexcept
Return an iterator pointing at the beginning of the collection.
LVL1::EFexTauAlgorithm::m_timeThr
float m_timeThr
Definition: EFexTauAlgorithm.h:63
CaloCell::eta
virtual double eta() const override final
get eta (through CaloDetDescrElement)
Definition: CaloCell.h:366
LVL1::EFexTauAlgorithm::~EFexTauAlgorithm
virtual ~EFexTauAlgorithm()
fitman.k
k
Definition: fitman.py:528
Generated on Thu Jul 4 2024 21:14:30 for ATLAS Offline Software by   doxygen 1.8.18