ATLAS Offline Software
DiTauTruthMatchingTool.cxx
Go to the documentation of this file.
1 
10 // Local include(s)
11 #include "TauAnalysisTools/DiTauTruthMatchingTool.h"
12 
13 // Core include(s):
14 #include "AthLinks/ElementLink.h"
15 #include "TruthUtils/HepMCHelpers.h"
16 #include "AthContainers/ConstAccessor.h"
17 #include "AthContainers/Decorator.h"
18 
19 #include "MCTruthClassifier/MCTruthClassifier.h"
20 
21 using namespace TauAnalysisTools;
22 
23 //=================================PUBLIC-PART==================================
24 //______________________________________________________________________________
25 DiTauTruthMatchingTool::DiTauTruthMatchingTool( const std::string& name )
26  : BuildTruthTaus(name)
27  , m_accPtVis("pt_vis")
28  , m_accEtaVis("eta_vis")
29  , m_accPhiVis("phi_vis")
30  , m_accMVis("m_vis")
31 {
32 }
33 
34 //______________________________________________________________________________
35 DiTauTruthMatchingTool::~DiTauTruthMatchingTool( )
36 {
37 
38 }
39 
40 //______________________________________________________________________________
41 StatusCode DiTauTruthMatchingTool::initialize()
42 {
43  ATH_MSG_INFO( "Initializing DiTauTruthMatchingTool" );
44 
45  // configure BuildTruthTaus in truth matching mode, not truth tau building mode
46  DiTauTruthMatchingTool::BuildTruthTaus::setTruthMatchingMode();
47 
48  if (DiTauTruthMatchingTool::BuildTruthTaus::initialize().isFailure())
49  {
50  ATH_MSG_FATAL("Failed initializing BuildTruthTaus");
51  return StatusCode::FAILURE;
52  }
53  return StatusCode::SUCCESS;
54 }
55 
56 //______________________________________________________________________________
57 void DiTauTruthMatchingTool::getTruth(const xAOD::DiTauJet& xDiTau)
58 {
59  if (retrieveTruthTaus().isFailure())
60  return;
61 
62  if (findTruthTau(xDiTau).isFailure())
63  ATH_MSG_WARNING("There was a failure in finding the matched truth tau");
64 
65  return;
66 }
67 
68 //______________________________________________________________________________
69 void DiTauTruthMatchingTool::getTruth(const std::vector<const xAOD::DiTauJet*>& vDiTaus)
70 {
71  for (auto xDiTau : vDiTaus)
72  getTruth(*xDiTau);
73  return;
74 }
75 
76 
78 // Private Part //
80 
81 //______________________________________________________________________________
82 StatusCode DiTauTruthMatchingTool::findTruthTau(const xAOD::DiTauJet& xDiTau)
83 {
84  // check if decorations were already added to the first passed tau
85  if (!m_bIsTruthMatchedAvailable.isValid()) {
86  static const SG::ConstAccessor<char> accIsTruthMatched("IsTruthMatched");
87  m_bIsTruthMatchedAvailable.set (accIsTruthMatched.isAvailable(xDiTau));
88  }
89  if (*m_bIsTruthMatchedAvailable.ptr())
90  return StatusCode::SUCCESS;
91 
92  if (m_bTruthTauAvailable)
93  return checkTruthMatch(xDiTau, *m_truthTausEvent.m_xTruthTauContainerConst);
94  else
95  return checkTruthMatch(xDiTau, *m_truthTausEvent.m_xTruthTauContainer);
96 }
97 
98 //______________________________________________________________________________
99 StatusCode DiTauTruthMatchingTool::checkTruthMatch (const xAOD::DiTauJet& xDiTau, const xAOD::TruthParticleContainer& xTruthTauContainer) const
100 {
101  std::vector<const xAOD::TruthParticle*> vTruthMatch;
102  std::vector<TruthMatchedParticleType> vTruthMatchedParticleType;
103  std::vector<const xAOD::Jet*> vTruthJetMatch;
104 
105 
106  xAOD::TruthParticleContainer xRemainingTruthTaus = xTruthTauContainer;
107 
108  static const SG::Decorator<char> decIsTruthMatched("IsTruthMatched");
109  static const SG::Decorator<char> decIsTruthHadronic("IsTruthHadronic");
110  static const SG::ConstAccessor<int> accNSubjets("n_subjets");
111  static const SG::ConstAccessor<char> accIsTruthHadronic("IsTruthHadronic");
112 
113  int n_subjets = 0;
114  if(!(accNSubjets.isAvailable(xDiTau))){
115  // n_subjets decoration is not available, recalculation on the fly
116  while (xDiTau.subjetPt(n_subjets) > 0. )
117  {
118  n_subjets++;
119  }
120  } else {
121  n_subjets = accNSubjets(xDiTau);
122  }
123 
124  // set default values for each subjet
125  for (int i = 0; i < n_subjets; ++i)
126  {
127  const xAOD::TruthParticle* xTruthMatch = nullptr;
128  TruthMatchedParticleType eTruthMatchedParticleType = Unknown;
129 
130  vTruthMatch.push_back(xTruthMatch);
131  vTruthMatchedParticleType.push_back(eTruthMatchedParticleType);
132 
133  const xAOD::Jet* xTruthJetMatch = nullptr;
134  vTruthJetMatch.push_back(xTruthJetMatch);
135  }
136 
137  // truthmatching for subjets:
138  for (int i = 0; i < n_subjets; ++i)
139  {
140  TLorentzVector vSubjetTLV;
141  vSubjetTLV.SetPtEtaPhiE(xDiTau.subjetPt(i),
142  xDiTau.subjetEta(i),
143  xDiTau.subjetPhi(i),
144  xDiTau.subjetE(i));
145  if ( truthMatch(vSubjetTLV,
146  xRemainingTruthTaus,
147  vTruthMatch.at(i),
148  vTruthJetMatch.at(i),
149  vTruthMatchedParticleType.at(i)).isFailure() )
150  {
151  ATH_MSG_WARNING("There was a failure in matching truth taus with subjet " << i);
152  return StatusCode::FAILURE;
153  }
154  if (vTruthMatch.at(i) && vTruthMatchedParticleType.at(i) == TruthHadronicTau)
155  {
156  xRemainingTruthTaus.erase( std::find(xRemainingTruthTaus.begin(),
157  xRemainingTruthTaus.end(),
158  vTruthMatch.at(i)) );
159  }
160  }
161 
162  // create links for jets
163  std::vector< ElementLink < xAOD::JetContainer > > vTruthJetLinks;
164  for (int i = 0; i < n_subjets; ++i)
165  {
166  const xAOD::Jet* xTruthJetMatch = vTruthJetMatch.at(i);
167  if(xTruthJetMatch){
168  ElementLink < xAOD::JetContainer > lTruthParticleLink(xTruthJetMatch, *m_truthTausEvent.m_xTruthJetContainerConst);
169  vTruthJetLinks.push_back(lTruthParticleLink);
170  }
171  else
172  {
173  ElementLink < xAOD::JetContainer > lTruthParticleLink;
174  vTruthJetLinks.push_back(lTruthParticleLink);
175  }
176  }
177  static const SG::Decorator<std::vector<ElementLink<xAOD::JetContainer>>>
178  decTruthJetLinks ("truthJetLinks");
179  decTruthJetLinks(xDiTau) = std::move(vTruthJetLinks);
180 
181 
182  bool bTruthMatched = true;
183 
184  // create link to the original TruthParticle
185  std::vector< ElementLink < xAOD::TruthParticleContainer > > vTruthLinks;
186  for (int i = 0; i < n_subjets; ++i)
187  {
188  const xAOD::TruthParticle* xTruthMatch = vTruthMatch.at(i);
189  TruthMatchedParticleType eTruthMatchedParticleType = vTruthMatchedParticleType.at(i);
190  if (xTruthMatch)
191  {
192  if (eTruthMatchedParticleType == TruthHadronicTau)
193  {
194  ElementLink < xAOD::TruthParticleContainer > lTruthParticleLink(xTruthMatch, xTruthTauContainer);
195  vTruthLinks.push_back(lTruthParticleLink);
196  }
197  else if (eTruthMatchedParticleType == TruthMuon)
198  {
199  ElementLink <xAOD::TruthParticleContainer> lTruthParticleLink(xTruthMatch, *m_truthTausEvent.m_xTruthMuonContainerConst);
200  vTruthLinks.push_back(lTruthParticleLink);
201  }
202  else if (eTruthMatchedParticleType == TruthElectron)
203  {
204  ElementLink <xAOD::TruthParticleContainer> lTruthParticleLink(xTruthMatch, *m_truthTausEvent.m_xTruthElectronContainerConst);
205  vTruthLinks.push_back(lTruthParticleLink);
206  }
207  }
208  else
209  {
210  ElementLink < xAOD::TruthParticleContainer > lTruthParticleLink;
211  vTruthLinks.push_back(lTruthParticleLink);
212 
213  // ditau is not truth matched if one of the two leading subjets is not truth matched
214  if (i == 0 || i == 1) bTruthMatched = false;
215  }
216  }
217 
218  static const SG::Decorator<std::vector<ElementLink<xAOD::TruthParticleContainer>>>
219  decTruthParticleLinks ("truthParticleLinks");
220  decTruthParticleLinks(xDiTau) = vTruthLinks;
221  if (!m_bTruthTauAvailable)
222  {
223  static const SG::Decorator<std::vector<ElementLink<xAOD::TruthParticleContainer>>>
224  decTruthTaus ("TruthTaus");
225  decTruthTaus(xDiTau) = vTruthLinks;
226  }
227 
228  static const SG::Decorator<float> decTruthLeadPt("TruthVisLeadPt");
229  static const SG::Decorator<float> decTruthLeadEta("TruthVisLeadEta");
230  static const SG::Decorator<float> decTruthLeadPhi("TruthVisLeadPhi");
231  static const SG::Decorator<float> decTruthLeadM("TruthVisLeadM");
232  static const SG::Decorator<float> decTruthLeadPdgID("TruthLeadPdgID");
233  static const SG::Decorator<float> decTruthSubleadPt("TruthVisSubleadPt");
234  static const SG::Decorator<float> decTruthSubleadEta("TruthVisSubleadEta");
235  static const SG::Decorator<float> decTruthSubleadPhi("TruthVisSubleadPhi");
236  static const SG::Decorator<float> decTruthSubleadM("TruthVisSubleadM");
237  static const SG::Decorator<float> decTruthSubleadPdgID("TruthSubleadPdgID");
238  static const SG::Decorator<float> decTruthDeltaR("TruthVisDeltaR");
239  static const SG::Decorator<float> decTruthMass("TruthVisMass");
240 
241  // the ditau candidate should have at least 2 subjets to be truth matched
242  if ( n_subjets < 2) {
243  decIsTruthMatched(xDiTau) = (char)false;
244  decIsTruthHadronic(xDiTau) = (char)false;
245  decTruthLeadPt(xDiTau) = -1234.;
246  decTruthLeadEta(xDiTau) = -1234.;
247  decTruthLeadPhi(xDiTau) = -1234.;
248  decTruthLeadM(xDiTau) = -1234.;
249  decTruthLeadPdgID(xDiTau) = -1234.;
250  decTruthSubleadPt(xDiTau) = -1234.;
251  decTruthSubleadEta(xDiTau) = -1234.;
252  decTruthSubleadPhi(xDiTau) = -1234.;
253  decTruthSubleadM(xDiTau) = -1234.;
254  decTruthSubleadPdgID(xDiTau) = -1234.;
255  decTruthDeltaR(xDiTau) = -1234.;
256  decTruthMass(xDiTau) = -1234.;
257  return StatusCode::SUCCESS;
258  }
259 
260  decIsTruthMatched(xDiTau) = (char)bTruthMatched;
261  if (bTruthMatched)
262  {
263  // ditau is hadronic if two leading subjets are truth matched with hadronic decay
264  decIsTruthHadronic(xDiTau) = (char)(vTruthMatchedParticleType[0]==TruthHadronicTau
265  && vTruthMatchedParticleType[1]==TruthHadronicTau );
266  }
267  else
268  decIsTruthHadronic(xDiTau) = (char)false;
269 
270  if (accIsTruthHadronic(xDiTau))
271  {
272  TLorentzVector tlvTruthTau1;
273  TLorentzVector tlvTruthTau2;
274  tlvTruthTau1.SetPtEtaPhiM(m_accPtVis(*(*vTruthLinks.at(0))),
275  m_accEtaVis(*(*vTruthLinks.at(0))),
276  m_accPhiVis(*(*vTruthLinks.at(0))),
277  m_accMVis(*(*vTruthLinks.at(0))));
278  tlvTruthTau2.SetPtEtaPhiM(m_accPtVis(*(*vTruthLinks.at(1))),
279  m_accEtaVis(*(*vTruthLinks.at(1))),
280  m_accPhiVis(*(*vTruthLinks.at(1))),
281  m_accMVis(*(*vTruthLinks.at(1))));
282 
283  decTruthLeadPt(xDiTau) = std::max(tlvTruthTau1.Pt(), tlvTruthTau2.Pt());
284  decTruthLeadEta(xDiTau) = (tlvTruthTau1.Pt() > tlvTruthTau2.Pt()) ? tlvTruthTau1.Eta() : tlvTruthTau2.Eta();
285  decTruthLeadPhi(xDiTau) = (tlvTruthTau1.Pt() > tlvTruthTau2.Pt()) ? tlvTruthTau1.Phi() : tlvTruthTau2.Phi();
286  decTruthLeadM(xDiTau) = (tlvTruthTau1.Pt() > tlvTruthTau2.Pt()) ? tlvTruthTau1.M() : tlvTruthTau2.M();
287  decTruthLeadPdgID(xDiTau) = (*vTruthLinks.at(0))->pdgId();
288  decTruthSubleadPt(xDiTau) = std::min(tlvTruthTau1.Pt(), tlvTruthTau2.Pt());
289  decTruthSubleadEta(xDiTau) = (tlvTruthTau1.Pt() > tlvTruthTau2.Pt()) ? tlvTruthTau2.Eta() : tlvTruthTau1.Eta();
290  decTruthSubleadPhi(xDiTau) = (tlvTruthTau1.Pt() > tlvTruthTau2.Pt()) ? tlvTruthTau2.Phi() : tlvTruthTau1.Phi();
291  decTruthSubleadM(xDiTau) = (tlvTruthTau1.Pt() > tlvTruthTau2.Pt()) ? tlvTruthTau2.M() : tlvTruthTau1.M();
292  decTruthSubleadPdgID(xDiTau) = (*vTruthLinks.at(1))->pdgId();
293  decTruthDeltaR(xDiTau) = tlvTruthTau1.DeltaR(tlvTruthTau2);
294  decTruthMass(xDiTau) = (tlvTruthTau1 + tlvTruthTau2).M();
295  }
296  else {
297  // set to a default value
298  decTruthLeadPt(xDiTau) = -1234.;
299  decTruthLeadEta(xDiTau) = -1234.;
300  decTruthLeadPhi(xDiTau) = -1234.;
301  decTruthLeadM(xDiTau) = -1234.;
302  decTruthLeadPdgID(xDiTau) = -1234.;
303  decTruthSubleadPt(xDiTau) = -1234.;
304  decTruthSubleadEta(xDiTau) = -1234.;
305  decTruthSubleadPhi(xDiTau) = -1234.;
306  decTruthSubleadM(xDiTau) = -1234.;
307  decTruthSubleadPdgID(xDiTau) = -1234.;
308  decTruthDeltaR(xDiTau) = -1234.;
309  decTruthMass(xDiTau) = -1234.;
310  }
311 
312  return StatusCode::SUCCESS;
313 }
314 //______________________________________________________________________________
315 StatusCode DiTauTruthMatchingTool::truthMatch(const TLorentzVector& vSubjetTLV,
316  const xAOD::TruthParticleContainer& xTruthTauContainer,
317  const xAOD::TruthParticle* &xTruthMatch,
318  const xAOD::Jet* &xTruthJetMatch,
319  TruthMatchedParticleType &eTruthMatchedParticleType) const
320 {
321  for (auto xTruthTauIt : xTruthTauContainer)
322  {
323  TLorentzVector vTruthVisTLV;
324  vTruthVisTLV.SetPtEtaPhiM(m_accPtVis(*xTruthTauIt),
325  m_accEtaVis(*xTruthTauIt),
326  m_accPhiVis(*xTruthTauIt),
327  m_accMVis(*xTruthTauIt));
328  if (vSubjetTLV.DeltaR(vTruthVisTLV) <= m_dMaxDeltaR)
329  {
330  static const SG::ConstAccessor<char> accIsHadronicTau("IsHadronicTau");
331  if (static_cast<bool>(accIsHadronicTau(*xTruthTauIt)))
332  eTruthMatchedParticleType = TruthHadronicTau;
333  else
334  continue; // don't let leptonic taus steal truthmatch just by chance
335 
336  xTruthMatch = xTruthTauIt;
337  break;
338  }
339  }
340 
341  if (!xTruthMatch and m_truthTausEvent.m_xTruthMuonContainerConst)
342  {
343  double dPtMax = 0.;
344  for (auto xTruthMuonIt : *m_truthTausEvent.m_xTruthMuonContainerConst)
345  {
346  if (vSubjetTLV.DeltaR(xTruthMuonIt->p4()) <= m_dMaxDeltaR)
347  {
348  if (xTruthMuonIt->pt()<dPtMax)
349  continue;
350  eTruthMatchedParticleType = TruthMuon;
351 
352  xTruthMatch = xTruthMuonIt;
353  dPtMax = xTruthMuonIt->pt();
354  }
355  }
356  }
357 
358  if (!xTruthMatch and m_truthTausEvent.m_xTruthElectronContainerConst)
359  {
360  double dPtMax = 0.;
361  for (auto xTruthElectronIt : *m_truthTausEvent.m_xTruthElectronContainerConst)
362  {
363  if (vSubjetTLV.DeltaR(xTruthElectronIt->p4()) <= m_dMaxDeltaR)
364  {
365  if (xTruthElectronIt->pt()<dPtMax)
366  continue;
367  eTruthMatchedParticleType = TruthElectron;
368  xTruthMatch = xTruthElectronIt;
369  dPtMax = xTruthElectronIt->pt();
370  }
371  }
372  }
373 
374  if (m_truthTausEvent.m_xTruthJetContainerConst)
375  {
376  double dPtMax = 0.;
377  for (auto xTruthJetIt : *m_truthTausEvent.m_xTruthJetContainerConst)
378  {
379  if (vSubjetTLV.DeltaR(xTruthJetIt->p4()) <= m_dMaxDeltaR)
380  {
381  if (xTruthJetIt->pt()<dPtMax)
382  continue;
383  xTruthJetMatch = xTruthJetIt;
384  dPtMax = xTruthJetIt->pt();
385  }
386  }
387  }
388 
389  return StatusCode::SUCCESS;
390 }
TauAnalysisTools
Definition: TruthCollectionMakerTau.h:16
TauAnalysisTools::TruthElectron
@ TruthElectron
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/TauAnalysisTools/Enums.h:102
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
TauAnalysisTools::getTruth
const xAOD::TruthParticle * getTruth(const xAOD::TauJet &xTau)
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:195
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
TauAnalysisTools::TruthHadronicTau
@ TruthHadronicTau
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/TauAnalysisTools/Enums.h:99
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
initialize
void initialize()
Definition: run_EoverP.cxx:894
xAOD::char
char
Definition: TrigDecision_v1.cxx:38
xAOD::DiTauJet_v1::subjetPhi
float subjetPhi(unsigned int numSubjet) const
Definition: DiTauJet_v1.cxx:111
SG::ConstAccessor< char >
MCTruthClassifier.h
xAOD::DiTauJet_v1::subjetE
float subjetE(unsigned int numSubjet) const
Definition: DiTauJet_v1.cxx:121
TauAnalysisTools::BuildTruthTaus
Definition: BuildTruthTaus.h:39
SG::Decorator< char >
lumiFormat.i
int i
Definition: lumiFormat.py:85
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
xAOD::TruthParticle_v1
Class describing a truth particle in the MC record.
Definition: TruthParticle_v1.h:37
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
TauAnalysisTools::TruthMuon
@ TruthMuon
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/TauAnalysisTools/Enums.h:101
xAOD::DiTauJet_v1::subjetEta
float subjetEta(unsigned int numSubjet) const
Definition: DiTauJet_v1.cxx:101
TauAnalysisTools::TruthMatchedParticleType
TruthMatchedParticleType
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/TauAnalysisTools/Enums.h:97
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
runIDPVM.pdgId
pdgId
Definition: runIDPVM.py:91
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
xAOD::Jet_v1
Class describing a jet.
Definition: Jet_v1.h:57
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
DiTauTruthMatchingTool.h
Tau, lepton and jet truth matching for ditau jets.
xAOD::DiTauJet_v1
Definition: DiTauJet_v1.h:31
xAOD::TruthParticle_v1::pt
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition: TruthParticle_v1.cxx:161
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
xAOD::DiTauJet_v1::subjetPt
float subjetPt(unsigned int numSubjet) const
Definition: DiTauJet_v1.cxx:91
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
DataVector::erase
iterator erase(iterator position)
Remove element at a given position.
xAOD::Jet_v1::pt
virtual double pt() const
The transverse momentum ( ) of the particle.
Definition: Jet_v1.cxx:44
Decorator.h
Helper class to provide type-safe access to aux data.
HepMCHelpers.h
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
Generated on Tue Sep 2 2025 21:10:04 for ATLAS Offline Software by   doxygen 1.8.18