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