ATLAS Offline Software
TauElectronVetoVariables.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 #ifndef XAOD_ANALYSIS
6 
8 
9 #include "xAODTau/TauJet.h"
13 
14 #include "GaudiKernel/SystemOfUnits.h"
15 
16 #include <algorithm>
17 #include <cmath>
18 #include <unordered_map>
19 #include <array>
20 
21 #include "TVector2.h"
22 
23 using Gaudi::Units::GeV;
24 
27 }
28 
29 
30 
32  ATH_CHECK( m_caloExtensionTool.retrieve() );
33  if (!m_ParticleCacheKey.key().empty()) {
34  ATH_CHECK(m_ParticleCacheKey.initialize());
35  } else {
36  m_useOldCalo = true;
37  }
38  return StatusCode::SUCCESS;
39 }
40 
41 
42 
44  if (pTau.nTracks() < 1) {
45  return StatusCode::SUCCESS;
46  }
47  ATH_MSG_DEBUG("in execute()");
48  float detPhiTrk = 0.;
49  float detEtaTrk = 0.;
50  float clEtaTrk = 0.;
51  float distEtaTrk = 0.;
52  float energy_3phi[101] = {0.0};
53  float eta[101] = {0.0};
54  int max1 = 0;
55  int max2 = 0;
56  int max = 0;
57  int n = 0;
58  float Emax1 = 0.;
59  float Emax2 = 0.;
60  float etamaxcut = 0.158;
61  float phimaxcut = 0.1;
62  float signum_eta = 0.;
63  signum_eta = pTau.track(0)->eta() / std::abs(pTau.track(0)->eta());
64  float sumETCellsLAr = 0.;
65  float eta0cut = 0.075;
66  float eta1cut = 0.0475;
67  float eta2cut = 0.075;
68  float eta3cut = 1.5;
69  float phi0cut = 0.3;
70  float phi1cut = 0.3;
71  float phi2cut = 0.075;
72  float phi3cut = 0.075;
73  float etareg = 0.;
74  float etacase1 = 1.8;
75  float etagran1 = 0.00315;
76  float etagran2 = 0.00415;
77  float sumETCellsHad1 = 0.;
78  float etahadcut = 0.2;
79  float phihadcut = 0.2;
80  const CaloCell *pCell;
81  int trackIndex = -1;
82 
83  //---------------------------------------------------------------------
84  // Calculate eta, phi impact point of leading track at calorimeter layers EM 0,1,2,3
85  //---------------------------------------------------------------------
86  Trk::TrackParametersIdHelper parsIdHelper;
87  constexpr size_t numberOfEM_Layers{4};
88  constexpr double invalidCoordinate{-11111.};
89  constexpr double invalidCoordinateThreshold{-11110.};
90  std::array<double, numberOfEM_Layers> extrapolatedEta{};
91  std::array<double, numberOfEM_Layers> extrapolatedPhi{};
92  extrapolatedEta.fill(invalidCoordinate);
93  extrapolatedPhi.fill(invalidCoordinate);
94 
95  /*get the CaloExtension object*/
96  const Trk::CaloExtension * caloExtension = nullptr;
97  std::unique_ptr<Trk::CaloExtension> uniqueExtension ;
98  const xAOD::TrackParticle *orgTrack = pTau.track(0)->track();
99  trackIndex = orgTrack->index();
100  if (m_useOldCalo) {
101  /* If CaloExtensionBuilder is unavailable, use the calo extension tool */
102  ATH_MSG_VERBOSE("Using the CaloExtensionTool");
103  uniqueExtension = m_caloExtensionTool->caloExtension(
104  Gaudi::Hive::currentContext(), *orgTrack);
105  caloExtension = uniqueExtension.get();
106  } else {
107  /*get the CaloExtension object*/
108  ATH_MSG_VERBOSE("Using the CaloExtensionBuilder Cache");
110  caloExtension = (*particleCache)[trackIndex];
111  ATH_MSG_VERBOSE("Getting element " << trackIndex << " from the particleCache");
112  if( not caloExtension ){
113  ATH_MSG_VERBOSE("Cache does not contain a calo extension -> Calculating with the a CaloExtensionTool" );
114  uniqueExtension = m_caloExtensionTool->caloExtension(
115  Gaudi::Hive::currentContext(), *orgTrack);
116  caloExtension = uniqueExtension.get();
117  }
118  }
119  if( not caloExtension){
120  ATH_MSG_DEBUG("extrapolation of leading track to calo surfaces failed : caloExtension is nullptr" );
121  return StatusCode::RECOVERABLE;
122  }
123  const std::vector<Trk::CurvilinearParameters>& clParametersVector = caloExtension->caloLayerIntersections();
124  if(clParametersVector.empty() ){
125  ATH_MSG_DEBUG("extrapolation of leading track to calo surfaces failed : caloLayerIntersection is empty" );
126  return StatusCode::RECOVERABLE;
127  }
128  // loop over calo layers
129  for( const Trk::CurvilinearParameters& cur : clParametersVector ){
130  // only use entry layer
131  if( !parsIdHelper.isEntryToVolume(cur.cIdentifier()) ) continue;
132  CaloSampling::CaloSample sample = parsIdHelper.caloSample(cur.cIdentifier());
133  int index = -1;
135  else if( sample == CaloSampling::EME1 || sample == CaloSampling::EMB1 ) index = 1;
136  else if( sample == CaloSampling::EME2 || sample == CaloSampling::EMB2 ) index = 2;
137  else if( sample == CaloSampling::EME3 || sample == CaloSampling::EMB3 ) index = 3;
138  if( index < 0 ) continue;
139  extrapolatedEta[index] = cur.position().eta();
140  extrapolatedPhi[index] = cur.position().phi();
141  }
142  for (size_t i = 0; i < numberOfEM_Layers; ++i) {
143  if ( extrapolatedEta[i] < invalidCoordinateThreshold || extrapolatedPhi[i] < invalidCoordinateThreshold ){
144  ATH_MSG_DEBUG("extrapolation of leading track to calo surfaces failed for sampling : " << i );
145  return StatusCode::SUCCESS;
146  }
147  }
148 
149  // Loop through jets, get links to clusters
150  std::bitset<200000> cellSeen{};
151  const std::unordered_map<int, int> samplingLookup{
152  {4,0}, {5,1}, {6,2}, {7,3}, {8,12},
153  {15, 12}, {16,13}, {17,14}, {18,12}, {19, 13}, {20,14}
154  };
155  const auto notFound{samplingLookup.end()};
156 
157  std::vector<xAOD::CaloVertexedTopoCluster> vertexedClusterList = pTau.vertexedClusters();
158  for (const xAOD::CaloVertexedTopoCluster& vertexedCluster : vertexedClusterList){
159 
160  const xAOD::CaloCluster& cluster = vertexedCluster.clust();
161  auto cell_links = cluster.getCellLinks();
162  if (cell_links == nullptr) {
163  ATH_MSG_DEBUG("NO Cell links found for cluster with pT " << cluster.pt());
164  continue;
165  }
166  CaloClusterCellLink::const_iterator pCellIter = cluster.getCellLinks()->begin();
167  CaloClusterCellLink::const_iterator pCellIterE = cluster.getCellLinks()->end();
168  for (; pCellIter != pCellIterE; ++pCellIter) {
169  double cellEta{}, cellPhi{}, cellET{};
170  pCell = *pCellIter;
171  if (cellSeen.test(pCell->caloDDE()->calo_hash())) continue;
172  else cellSeen.set(pCell->caloDDE()->calo_hash());
173  if (m_doVertexCorrection && pTau.vertex()!=nullptr) {
174  CaloVertexedCell vxCell (*pCell, pTau.vertex()->position());
175  cellPhi = vxCell.phi();
176  cellEta = vxCell.eta();
177  cellET = vxCell.et();
178  } else {
179  cellPhi = pCell->phi();
180  cellEta = pCell->eta();
181  cellET = pCell->et();
182  }
183  int sampling = pCell->caloDDE()->getSampling();
184  if (const auto & pElement {samplingLookup.find(sampling)};pElement != notFound){
185  sampling = pElement->second;
186  }
187  int i = 2;
188  if (sampling < 4) i = sampling;
189  if (sampling == 12 || sampling == 13 || sampling == 14) i = 3;
190  detPhiTrk = TVector2::Phi_mpi_pi(cellPhi-extrapolatedPhi[i]);
191  detEtaTrk = std::abs( cellEta - extrapolatedEta[i] );
192  clEtaTrk = extrapolatedEta[i];
193  distEtaTrk = cellEta - extrapolatedEta[i];
194  if ((sampling == 0 && detEtaTrk < eta0cut && detPhiTrk < phi0cut) ||
195  (sampling == 1 && detEtaTrk < eta1cut && detPhiTrk < phi1cut) ||
196  (sampling == 2 && detEtaTrk < eta2cut && detPhiTrk < phi2cut) ||
197  (sampling == 3 && detEtaTrk < eta3cut && detPhiTrk < phi3cut)) {
198  sumETCellsLAr += cellET;
199  }
200  if (sampling == 12 && detEtaTrk < etahadcut && detPhiTrk < phihadcut) sumETCellsHad1 += cellET;
201  if (std::abs(cellEta) > 0.8 && std::abs(cellEta) <= 1.2 && (sampling == 13 || sampling == 14) && detEtaTrk < etahadcut && detPhiTrk < phihadcut) {
202  sumETCellsHad1 += cellET;
203  }
204  if (std::abs(pTau.track(0)->eta()) <= 1.7) {
205  if (sampling == 1 && detEtaTrk < etamaxcut && detPhiTrk <= phimaxcut) {
206  if ((std::abs(cellEta) < 1.37 || std::abs(cellEta) > 1.52) && std::abs(cellEta) < 1.85) {
207  if (std::abs(clEtaTrk) <= etacase1 && std::abs(cellEta) <= etacase1) {
208  n = 50 + int(distEtaTrk / etagran1);
209  }
210  if (std::abs(clEtaTrk) <= etacase1 && std::abs(cellEta) > etacase1) {
211  n = 50 + int(signum_eta * ((etacase1 - std::abs(clEtaTrk)) / etagran1 + (-etacase1 + std::abs(cellEta)) / etagran2));
212  }
213  energy_3phi[n] = energy_3phi[n] + cellET / GeV;
214  eta[n] = signum_eta * (clEtaTrk - cellEta);
215  } else {
216  energy_3phi[n] = 0;
217  eta[n] = 0;
218  }
219  }
220  } else {
221  energy_3phi[n] = 0;
222  eta[n] = 0;
223  }
224  if (std::abs(cellEta) <= etacase1) {
225  etareg = 0.00315;
226  } else {
227  etareg = 0.00415;
228  }
229  } //end cell loop
230  }// end jet constituent loop
231  for (int m1 = 0; m1 < 101; m1++) {
232  if ((energy_3phi[m1] > Emax1)) {
233  Emax1 = energy_3phi[m1];
234  max1 = m1;
235  }
236  }
237  for (int m2 = 1; m2 < 100; m2++) {
238  if (m2 == max1) continue;
239  if ((energy_3phi[m2] > Emax2) && (energy_3phi[m2] > energy_3phi[m2 - 1]) && (energy_3phi[m2] > energy_3phi[m2 + 1])) {
240  Emax2 = energy_3phi[m2];
241  max2 = m2;
242  }
243  }
244  if (std::abs(eta[max1]) >= etareg) {
245  max = max1;
246  } else {
247  max = max2;
248  }
249 
250  const xAOD::TrackParticle* leadTrack = pTau.track(0)->track();
251 
252  uint8_t TRTHits{0};
253  if ( !leadTrack->summaryValue( TRTHits, xAOD::SummaryType::numberOfTRTHits ) ){
254  ATH_MSG_DEBUG("retrieval of track summary value failed. Not filling electron veto variables for this one prong candidate");
255  return StatusCode::SUCCESS;
256  }
257 
258  uint8_t TRTHTHits{0};
259  if ( !leadTrack->summaryValue( TRTHTHits, xAOD::SummaryType::numberOfTRTHighThresholdHits ) ){
260  ATH_MSG_DEBUG("retrieval of track summary value failed. Not filling electron veto variables for this one prong candidate");
261  return StatusCode::SUCCESS;
262  }
263 
264  uint8_t TRTOutliers{0};
265  if ( !leadTrack->summaryValue( TRTOutliers, xAOD::SummaryType::numberOfTRTOutliers ) ){
266  ATH_MSG_DEBUG("retrieval of track summary value failed. Not filling electron veto variables for this one prong candidate");
267  return StatusCode::SUCCESS;
268  }
269 
270  uint8_t TRTHTOutliers{0};
271  if ( !leadTrack->summaryValue( TRTHTOutliers, xAOD::SummaryType::numberOfTRTHighThresholdOutliers ) ) {
272  ATH_MSG_DEBUG("retrieval of track summary value failed. Not filling electron veto variables for this one prong candidate");
273  return StatusCode::SUCCESS;
274  }
275 
276  float TRTratio = -9999.0;
277  if (TRTHits + TRTOutliers != 0) {
278  TRTratio = float( TRTHTHits + TRTHTOutliers) / float( TRTHits + TRTOutliers );
279  } else {
280  TRTratio = 0.0;
281  }
282 
285  pTau.setDetail(xAOD::TauJetParameters::hadLeakEt , static_cast<float>( sumETCellsHad1 / ( pTau.track(0)->pt() ) ) );
286  pTau.setDetail(xAOD::TauJetParameters::sumEMCellEtOverLeadTrkPt , static_cast<float>( ( sumETCellsLAr / ( pTau.track(0)->pt() ) ) ) );
287  return StatusCode::SUCCESS;
288 }
289 
290 #endif
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:18
CaloCell::phi
virtual double phi() const override final
get phi (through CaloDetDescrElement)
Definition: CaloCell.h:369
CaloVertexedCell::eta
virtual double eta() const final
The pseudorapidity of the particle.
Definition: CaloVertexedCell.h:65
constants.EMB1
int EMB1
Definition: Calorimeter/CaloClusterCorrection/python/constants.py:53
TileDCSDataPlotter.max1
max1
Definition: TileDCSDataPlotter.py:886
xAOD::TauJetParameters::hadLeakEt
@ hadLeakEt
Definition: TauDefs.h:232
xAOD::uint8_t
uint8_t
Definition: Muon_v1.cxx:558
TauElectronVetoVariables::m_useOldCalo
Gaudi::Property< bool > m_useOldCalo
Definition: TauElectronVetoVariables.h:44
Trk::CaloExtension
Tracking class to hold the extrapolation through calorimeter Layers Both the caloEntryLayerIntersecti...
Definition: CaloExtension.h:19
beamspotman.cur
def cur
Definition: beamspotman.py:669
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:67
index
Definition: index.py:1
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
xAOD::TauTrack_v1::eta
virtual double eta() const
The pseudorapidity ( ) of the particle.
xAOD::TrackParticle_v1::summaryValue
bool summaryValue(uint8_t &value, const SummaryType &information) const
Accessor for TrackSummary values.
Definition: TrackParticle_v1.cxx:737
xAOD::TauJet_v3::nTracks
size_t nTracks(TauJetParameters::TauTrackFlag flag=TauJetParameters::TauTrackFlag::classifiedCharged) const
Definition: TauJet_v3.cxx:488
TauRecToolBase
The base class for all tau tools.
Definition: TauRecToolBase.h:21
Trk::TrackParametersIdHelper::caloSample
CaloSampling::CaloSample caloSample(TrackParametersIdentifier id) const
CaloSample encoded in id, returns CaloSampling::Unknown if id is not valid
Definition: TrackParametersIdHelper.h:91
xAOD::numberOfTRTHits
@ numberOfTRTHits
number of TRT hits [unit8_t].
Definition: TrackingPrimitives.h:276
xAOD::TauJetParameters::secMaxStripEt
@ secMaxStripEt
migrate only seedTrk_ variables which are used in reco and ID and without prefix
Definition: TauDefs.h:230
xAOD::Vertex_v1::position
const Amg::Vector3D & position() const
Returns the 3-pos.
CaloVertexedCell::phi
virtual double phi() const final
The aximuthal angle of the particle.
Definition: CaloVertexedCell.h:68
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Phi_mpi_pi
__HOSTDEV__ double Phi_mpi_pi(double)
Definition: GeoRegion.cxx:7
Trk::TrackParametersIdHelper::isEntryToVolume
bool isEntryToVolume(TrackParametersIdentifier id) const
returns true if the id belongs to the volume entrance
Definition: TrackParametersIdHelper.h:81
xAOD::numberOfTRTHighThresholdHits
@ numberOfTRTHighThresholdHits
number of TRT hits which pass the high threshold (only xenon counted) [unit8_t].
Definition: TrackingPrimitives.h:279
Trk::TrackParametersIdHelper
helper class to encode and decode a TrackParametersIdentifier
Definition: TrackParametersIdHelper.h:18
xAOD::numberOfTRTHighThresholdOutliers
@ numberOfTRTHighThresholdOutliers
number of TRT high threshold outliers (only xenon counted) [unit8_t].
Definition: TrackingPrimitives.h:282
xAOD::CaloCluster_v1
Description of a calorimeter cluster.
Definition: CaloCluster_v1.h:62
TauElectronVetoVariables::m_doVertexCorrection
Gaudi::Property< bool > m_doVertexCorrection
Definition: TauElectronVetoVariables.h:43
CaloDetDescrElement::calo_hash
IdentifierHash calo_hash() const
cell calo hash
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:412
python.changerun.m1
m1
Definition: changerun.py:30
TrackParametersIdHelper.h
TauElectronVetoVariables::m_ParticleCacheKey
SG::ReadHandleKey< CaloExtensionCollection > m_ParticleCacheKey
Definition: TauElectronVetoVariables.h:47
constants.EMB2
int EMB2
Definition: Calorimeter/CaloClusterCorrection/python/constants.py:54
FullCPAlgorithmsTest_eljob.sample
sample
Definition: FullCPAlgorithmsTest_eljob.py:116
TauElectronVetoVariables::execute
virtual StatusCode execute(xAOD::TauJet &pTau) const override
Execute - called for each tau candidate.
Definition: TauElectronVetoVariables.cxx:43
lumiFormat.i
int i
Definition: lumiFormat.py:85
CaloSampling::CaloSample
CaloSample
Definition: Calorimeter/CaloGeoHelpers/CaloGeoHelpers/CaloSampling.h:22
beamspotman.n
n
Definition: beamspotman.py:729
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
CaloCell::caloDDE
const CaloDetDescrElement * caloDDE() const
get pointer to CaloDetDescrElement (data member)
Definition: CaloCell.h:315
xAOD::TauJet_v3
Class describing a tau jet.
Definition: TauJet_v3.h:41
xAOD::TauJet_v3::track
const TauTrack * track(size_t i, TauJetParameters::TauTrackFlag flag=TauJetParameters::TauTrackFlag::classifiedCharged, int *container_index=0) const
Get the pointer to a given tauTrack associated with this tau /*container index needed by trackNonCons...
Definition: TauJet_v3.cxx:422
CaloCell::et
virtual double et() const override final
get et
Definition: CaloCell.h:417
constants.EME1
int EME1
Definition: Calorimeter/CaloClusterCorrection/python/constants.py:55
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
TauElectronVetoVariables::TauElectronVetoVariables
TauElectronVetoVariables(const std::string &name)
Definition: TauElectronVetoVariables.cxx:25
Trk::CurvilinearParametersT
Definition: CurvilinearParametersT.h:48
SG::AuxElement::index
size_t index() const
Return the index of this element within its container.
xAOD::CaloCluster_v1::getCellLinks
const CaloClusterCellLink * getCellLinks() const
Get a pointer to the CaloClusterCellLink object (const version)
Definition: CaloCluster_v1.cxx:859
xAOD::CaloCluster_v1::pt
virtual double pt() const
The transverse momentum ( ) of the particle (negative for negative-energy clusters)
Definition: CaloCluster_v1.cxx:247
xAOD::TauTrack_v1::pt
virtual double pt() const
The transverse momentum ( ) of the particle.
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
CaloCell_ID_FCS::EME3
@ EME3
Definition: FastCaloSim_CaloCell_ID.h:26
TauElectronVetoVariables.h
xAOD::TauJet_v3::vertexedClusters
std::vector< xAOD::CaloVertexedTopoCluster > vertexedClusters() const
Definition: TauJet_v3.cxx:586
xAOD::TauJetParameters::TRT_NHT_OVER_NLT
@ TRT_NHT_OVER_NLT
TRT hits high threshold over low threshold.
Definition: TauDefs.h:260
xAOD::numberOfTRTOutliers
@ numberOfTRTOutliers
number of TRT outliers [unit8_t].
Definition: TrackingPrimitives.h:277
IParticleCaloExtensionTool.h
DeMoScan.index
string index
Definition: DeMoScan.py:362
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
CaloCell
Data object for each calorimeter readout cell.
Definition: CaloCell.h:57
CaloVertexedCell
Evaluate cell kinematics with a different vertex.
Definition: CaloVertexedCell.h:37
TauElectronVetoVariables::initialize
virtual StatusCode initialize() override
Tool initializer.
Definition: TauElectronVetoVariables.cxx:31
CaloDetDescrElement::getSampling
CaloCell_ID::CaloSample getSampling() const
cell sampling
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:395
xAOD::TauJet_v3::vertex
const Vertex * vertex() const
xAOD::TauJetParameters::sumEMCellEtOverLeadTrkPt
@ sumEMCellEtOverLeadTrkPt
Definition: TauDefs.h:231
CaloCell_ID_FCS::PreSamplerE
@ PreSamplerE
Definition: FastCaloSim_CaloCell_ID.h:23
CaloCell_ID_FCS::PreSamplerB
@ PreSamplerB
Definition: FastCaloSim_CaloCell_ID.h:19
TauJet.h
xAOD::TauJet_v3::setDetail
void setDetail(TauJetParameters::Detail detail, int value)
Definition: TauJet_v3.cxx:309
P4EEtaPhiMBase::et
virtual double et() const
transverse energy defined to be e*sin(theta)
Definition: P4EEtaPhiMBase.cxx:106
xAOD::TauTrack_v1::track
const TrackParticle * track() const
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
xAOD::CaloVertexedTopoCluster
Evaluate cluster kinematics with a different vertex / signal state.
Definition: Event/xAOD/xAODCaloEvent/xAODCaloEvent/CaloVertexedTopoCluster.h:38
python.SystemOfUnits.m2
float m2
Definition: SystemOfUnits.py:107
CaloCell_ID_FCS::EMB3
@ EMB3
Definition: FastCaloSim_CaloCell_ID.h:22
constants.EME2
int EME2
Definition: Calorimeter/CaloClusterCorrection/python/constants.py:56
CaloVertexedCell.h
Evaluate cell kinematics with a different vertex.
CaloCell::eta
virtual double eta() const override final
get eta (through CaloDetDescrElement)
Definition: CaloCell.h:376
python.LArMinBiasAlgConfig.float
float
Definition: LArMinBiasAlgConfig.py:65
TauElectronVetoVariables::m_caloExtensionTool
ToolHandle< Trk::IParticleCaloExtensionTool > m_caloExtensionTool
Definition: TauElectronVetoVariables.h:45