ATLAS Offline Software
PseudoTopReco.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3  */
4 
5 // definitions taken from the MAMbo code here: https://twiki.cern.ch/twiki/bin/view/Main/MAMbo
6 
7 // Local include(s):
8 #include "TopEventReconstructionTools/PseudoTopReco.h"
9 #include "TopEvent/Event.h"
10 #include "TopEvent/EventTools.h"
11 #include "TopEvent/PseudoTopResultContainer.h"
12 #include "TopConfiguration/TopConfig.h"
13 #include "PathResolver/PathResolver.h"
14 #include "xAODCore/AuxContainerBase.h"
15 
16 #include <algorithm>
17 
18 namespace top {
19  PseudoTopReco::PseudoTopReco(const std::string& name) :
20  asg::AsgTool(name),
21  m_config(nullptr),
22  // commented out variables are unused. experts please check and remove
23  // m_bTagCutValue(9999.9),
24  m_leptonType("SetMe") {
25  declareProperty("config", m_config, "Set the configuration");
26  declareProperty("LeptonType", m_leptonType = "kUndefined", "Define the lepton type");
27  }
28 
30  StatusCode PseudoTopReco::initialize() {
31  // Have you set the config??
32  if (m_config == nullptr) {
33  ATH_MSG_ERROR("Please set the top::TopConfig");
34  return StatusCode::FAILURE;
35  }
36 
37  m_config->setPseudoTop();
38 
39  // Figure out the b tagging working point
40  if (m_config->bTagWP().size() != 1) {
41  ATH_MSG_INFO(
42  m_config->bTagWP().size() <<
43  " b-tagging WP - cannot pick b-jets. Please select only 1 WP if you want to use the PseudoTop reconstruction");
44  }
45 
46 
47  ATH_MSG_INFO("++++++++++++++++++++++++++++++");
48  ATH_MSG_INFO(" Using Lepton \t" << m_leptonType);
49  ATH_MSG_INFO("++++++++++++++++++++++++++++++");
50 
52  return StatusCode::SUCCESS;
53  }
54 
56  StatusCode PseudoTopReco::execute(const top::Event& event) {
57  if ((!event.m_isLoose &&
58  evtStore()->contains<xAOD::PseudoTopResultContainer>(m_config->sgKeyPseudoTop(event.m_hashValue))) ||
59  (event.m_isLoose &&
60  evtStore()->contains<xAOD::PseudoTopResultContainer>(m_config->sgKeyPseudoTopLoose(event.m_hashValue)))) {
61  return StatusCode::SUCCESS;
62  }
63 
64  // Create the partonHistory xAOD object
65  xAOD::AuxContainerBase* pseudoTopAuxCont = new xAOD::AuxContainerBase {};
66  xAOD::PseudoTopResultContainer* pseudoTop = new xAOD::PseudoTopResultContainer {};
67  pseudoTop->setStore(pseudoTopAuxCont);
68 
69  xAOD::PseudoTopResult* PseudoTopResult = new xAOD::PseudoTopResult {};
70  pseudoTop->push_back(PseudoTopResult);
71 
72  PseudoTopResult->IniVar(true);
73 
74  // fill either reco particles or truth particles into four vectors
75  // store them in GeV
76  // check if particle/jet is b-tagged
77  // check if electron or muon channel
78  // store output so that it can be written later to the event saver
79 
80  m_bJets.clear();
81  m_lightJets.clear();
82 
83  // get the MET x and y components for the reconstruction
84  m_nu_px = event.m_met->mpx() / 1.e3;
85  m_nu_py = event.m_met->mpy() / 1.e3;
86  m_met_et = sqrt(m_nu_px * m_nu_px + m_nu_py * m_nu_py);
87 
88  SetJetInfo(event);
89  SetChargedLeptonInfo(event);
90 
91  if (m_bJets.size() == 2 && m_lightJets.size() >= 2) {
92  ReconstructLeptonicW();
93  RunReconstruction();
94 
95  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_top_lep_pt") = m_top_lep.Pt();
96  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_top_lep_eta") = m_top_lep.Eta();
97  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_top_lep_phi") = m_top_lep.Phi();
98  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_top_lep_m") = m_top_lep.M();
99 
100  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_top_had_pt") = m_top_had.Pt();
101  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_top_had_eta") = m_top_had.Eta();
102  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_top_had_phi") = m_top_had.Phi();
103  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_top_had_m") = m_top_had.M();
104 
105  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_ttbar_pt") = m_ttbar.Pt();
106  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_ttbar_eta") = m_ttbar.Eta();
107  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_ttbar_phi") = m_ttbar.Phi();
108  PseudoTopResult->auxdecor< float >("PseudoTop_Reco_ttbar_m") = m_ttbar.M();
109  } else {
110  // std::cout<< "NOT ENOUGH BJETS!!!" << std::endl;
111  }
112 
113  // Save to StoreGate / TStore
114  std::string outputSGKey("SetMe");
115  if (!event.m_isLoose) {
116  outputSGKey = m_config->sgKeyPseudoTop(event.m_hashValue);
117  } else {
118  outputSGKey = m_config->sgKeyPseudoTopLoose(event.m_hashValue);
119  }
120  std::string outputSGKeyAux = outputSGKey + "Aux.";
121 
122  StatusCode save = evtStore()->tds()->record(pseudoTop, outputSGKey);
123  StatusCode saveAux = evtStore()->tds()->record(pseudoTopAuxCont, outputSGKeyAux);
124  if (!save || !saveAux) {
125  return StatusCode::FAILURE;
126  }
127 
129  return StatusCode::SUCCESS;
130  }
131 
132  StatusCode PseudoTopReco::execute(const top::ParticleLevelEvent& plEvent) {
133  if (evtStore()->contains<xAOD::PseudoTopResultContainer>(m_config->sgKeyPseudoTop(0))) {
134  return StatusCode::SUCCESS;
135  }
136 
137  // Create the pseudoTopHistory xAOD object
138  xAOD::AuxContainerBase* pseudoTopAuxCont = new xAOD::AuxContainerBase {};
139  xAOD::PseudoTopResultContainer* pseudoTop = new xAOD::PseudoTopResultContainer {};
140  pseudoTop->setStore(pseudoTopAuxCont);
141 
142  xAOD::PseudoTopResult* PseudoTopResult = new xAOD::PseudoTopResult {};
143  pseudoTop->push_back(PseudoTopResult);
144 
145  PseudoTopResult->IniVar(false);
146 
147 
148  m_bJets.clear();
149  m_lightJets.clear();
150 
151  // get the MET x and y components for the reconstruction
152  m_nu_px = plEvent.m_met->mpx() / 1.e3;
153  m_nu_py = plEvent.m_met->mpy() / 1.e3;
154  m_met_et = sqrt(m_nu_px * m_nu_px + m_nu_py * m_nu_py);
155 
156  SetJetInfo(plEvent);
157  SetChargedLeptonInfo(plEvent);
158 
159  if (m_bJets.size() == 2 && m_lightJets.size() >= 2) {
160  ReconstructLeptonicW();
161  RunReconstruction();
162 
163  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_top_lep_pt") = m_top_lep.Pt();
164  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_top_lep_eta") = m_top_lep.Eta();
165  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_top_lep_phi") = m_top_lep.Phi();
166  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_top_lep_m") = m_top_lep.M();
167 
168  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_top_had_pt") = m_top_had.Pt();
169  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_top_had_eta") = m_top_had.Eta();
170  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_top_had_phi") = m_top_had.Phi();
171  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_top_had_m") = m_top_had.M();
172 
173  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_ttbar_pt") = m_ttbar.Pt();
174  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_ttbar_eta") = m_ttbar.Eta();
175  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_ttbar_phi") = m_ttbar.Phi();
176  PseudoTopResult->auxdecor< float >("PseudoTop_Particle_ttbar_m") = m_ttbar.M();
177  } else {
178  // std::cout << "NOT ENOUGH BJETS!!!" << std::endl;
179  }
180 
181 
182  // Save to StoreGate / TStore
183  std::string outputSGKey = m_config->sgKeyPseudoTop(0);
184  std::string outputSGKeyAux = outputSGKey + "Aux.";
185 
186  StatusCode save = evtStore()->tds()->record(pseudoTop, outputSGKey);
187  StatusCode saveAux = evtStore()->tds()->record(pseudoTopAuxCont, outputSGKeyAux);
188  if (!save || !saveAux) {
189  return StatusCode::FAILURE;
190  }
191 
193  return StatusCode::SUCCESS;
194  }
195 
196  bool PseudoTopReco::SetChargedLeptonInfo(const top::Event& event) {
197  if (m_leptonType == "kElectron") {
198  m_lepton.SetPtEtaPhiE(event.m_electrons.at(0)->pt() / 1.e3, event.m_electrons.at(0)->eta(), event.m_electrons.at(
199  0)->phi(), event.m_electrons.at(0)->e() / 1.e3);
200  } else if (m_leptonType == "kMuon") {
201  m_lepton.SetPtEtaPhiE(event.m_muons.at(0)->pt() / 1.e3, event.m_muons.at(0)->eta(), event.m_muons.at(
202  0)->phi(), event.m_muons.at(0)->e() / 1.e3);
203  } else {
204  ATH_MSG_ERROR(" Please supply a valid LeptonType : kElectron or kMuon");
205  return false;
206  }
207 
208  return true;
209  }
210 
211  bool PseudoTopReco::SetChargedLeptonInfo(const top::ParticleLevelEvent& plEvent) {
212  if (m_leptonType == "kElectron") {
213  m_lepton.SetPtEtaPhiE(plEvent.m_electrons->at(0)->pt() / 1.e3, plEvent.m_electrons->at(
214  0)->eta(), plEvent.m_electrons->at(0)->phi(), plEvent.m_electrons->at(0)->e() / 1.e3);
215  } else if (m_leptonType == "kMuon") {
216  m_lepton.SetPtEtaPhiE(plEvent.m_muons->at(0)->pt() / 1.e3, plEvent.m_muons->at(0)->eta(), plEvent.m_muons->at(
217  0)->phi(), plEvent.m_muons->at(0)->e() / 1.e3);
218  } else {
219  ATH_MSG_ERROR(" Please supply a valid LeptonType : kElectron or kMuon");
220  return false;
221  }
222 
223  return true;
224  }
225 
226  bool PseudoTopReco::SetJetInfo(const top::Event& event) {
227  int bCounter = 0;
228 
229  for (const auto* const jetPtr : event.m_jets) {
230  TLorentzVector helpVec(0, 0, 0, 0);
231  helpVec.SetPtEtaPhiE(jetPtr->pt() / 1.e3, jetPtr->eta(), jetPtr->phi(), jetPtr->e() / 1.e3);
232 
233  // if more than two b-jets are available, consider third leading b-tagged jet as light jet!
234  std::string out = m_config->bTagWP()[0];
235 
236  const bool hasbTagFlag = jetPtr->isAvailable<char>("isbtagged_" + out);
237 
238  if (hasbTagFlag) {
239  if (jetPtr->auxdataConst<char>("isbtagged_" + out) && bCounter < 2) {
240  m_bJets.push_back(helpVec);
241 
242  bCounter++;
243  } else {
244  m_lightJets.push_back(helpVec);
245  }
246  }
247  } // end loop over jets
248 
249 
250  return true;
251  }
252 
253  bool PseudoTopReco::SetJetInfo(const top::ParticleLevelEvent& plEvent) {
254  // if more than two b-jets are available, consider third leading b-tagged jet as light jet!
255 
256  int bCounter = 0;
257 
258  for (const auto *jetPtr : *plEvent.m_jets) {
259  TLorentzVector helpVec(0, 0, 0, 0);
260  helpVec.SetPtEtaPhiE(jetPtr->pt() / 1.e3, jetPtr->eta(), jetPtr->phi(), jetPtr->e() / 1.e3);
261 
262  int nGhosts = jetPtr->auxdata<int>("GhostBHadronsFinalCount");
263 
264  if (nGhosts >= 1 && bCounter < 2) {
265  m_bJets.push_back(helpVec);
266  bCounter++;
267  } else m_lightJets.push_back(helpVec);
268  }
269 
270  return true;
271  }
272 
273  bool PseudoTopReco::RunReconstruction() {
274  double bestwmass = 999.e9;
275  TLorentzVector Whelp(0, 0, 0, 0);
276 
277  // get first the b-jet closest to the lepton in dR
278  // define this as TLorentzVector m_b_lep;
279  // take now remaining b-jet as m_b_had
280  double dRb1 = m_bJets[0].DeltaR(m_lepton);
281  double dRb2 = m_bJets[1].DeltaR(m_lepton);
282 
283  if (dRb1 < dRb2) {
284  m_b_lep = m_bJets[0];
285  m_b_had = m_bJets[1];
286  } else {
287  m_b_lep = m_bJets[1];
288  m_b_had = m_bJets[0];
289  }
290 
291  // then loop over all light jets
292  // take combination which is closest to mW (PDG)
293  for (unsigned int iJet = 0; iJet < m_lightJets.size(); ++iJet) {
294  for (unsigned int kJet = 0; kJet < m_lightJets.size(); ++kJet) {
295  if (iJet == kJet) continue;
296 
297  Whelp = m_lightJets[iJet] + m_lightJets[kJet];
298 
299  if (fabs(mWPDG - Whelp.M()) < fabs(mWPDG - bestwmass)) {
300  m_W_had = Whelp;
301  bestwmass = Whelp.M();
302  }
303  }
304  }
305 
306  // define now the top and ttbar four-vectors
307  m_top_lep = m_W_lep + m_b_lep;
308  m_top_had = m_W_had + m_b_had;
309  m_ttbar = m_top_lep + m_top_had;
310 
311  return true;
312  }
313 
314  bool PseudoTopReco::ReconstructLeptonicW() {
315  double v_pz = -KinemEdge;
316  double delta = -KinemEdge;
317 
318  // lepton already made by reco or particle levels:
319  double l_px = m_lepton.Px();
320  double l_py = m_lepton.Py();
321  double l_pz = m_lepton.Pz();
322  double l_m = m_lepton.M();
323  double l_E = m_lepton.E();
324 
325  double mdiff = 0.5 * (mWPDG * mWPDG - l_m * l_m);
326  double pT_vl = m_nu_px * l_px + m_nu_py * l_py;
327 
328  double a = l_E * l_E - l_pz * l_pz;
329  double b = -2. * l_pz * (mdiff + pT_vl);
330  double c = m_met_et * m_met_et * l_E * l_E - mdiff * mdiff - pT_vl * pT_vl - 2. * mdiff * pT_vl;
331 
332  delta = b * b - 4. * a * c;
333 
334  if (delta <= 0.) {
335  // take only the real part
336  v_pz = -0.5 * b / a;
337  } else {
338  double v_pz_1 = 0.5 * (-b - sqrt(delta)) / a;
339  double v_pz_2 = 0.5 * (-b + sqrt(delta)) / a;
340 
341  v_pz = (fabs(v_pz_1) > fabs(v_pz_2)) ? v_pz_2 : v_pz_1;
342  }
343 
344  double v_E = sqrt(m_met_et * m_met_et + v_pz * v_pz);
345 
346  m_neutrino.SetPxPyPzE(m_nu_px, m_nu_py, v_pz, v_E);
347 
348  m_W_lep = m_neutrino + m_lepton;
349 
350  return true;
351  }
352 
354  StatusCode PseudoTopReco::finalize() {
356  return StatusCode::SUCCESS;
357  }
358 }
top::ParticleLevelEvent::m_muons
const xAOD::TruthParticleContainer * m_muons
Pointer to truth level muons.
Definition: ParticleLevelEvent.h:45
top::PseudoTopReco::m_ttbar
TLorentzVector m_ttbar
Definition: PseudoTopReco.h:85
top
TopConfig A simple configuration that is NOT a singleton.
Definition: AnalysisTrackingHelper.cxx:58
top::PseudoTopReco::m_nu_px
double m_nu_px
Definition: PseudoTopReco.h:87
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
PseudoTopResultContainer.h
top::PseudoTopReco::m_nu_py
double m_nu_py
Definition: PseudoTopReco.h:87
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
top::PseudoTopReco::m_met_et
double m_met_et
Definition: PseudoTopReco.h:87
AuxContainerBase.h
top::ParticleLevelEvent::m_jets
const xAOD::JetContainer * m_jets
Pointer to truth level jets.
Definition: ParticleLevelEvent.h:54
xAOD::AuxContainerBase
Common base class for the auxiliary containers.
Definition: AuxContainerBase.h:59
asg
Definition: DataHandleTestTool.h:28
python.AthDsoLogger.out
out
Definition: AthDsoLogger.py:71
top::PseudoTopReco::execute
virtual StatusCode execute(const top::Event &)
Function executing the tool.
Definition: PseudoTopReco.cxx:56
SG::AuxElement::auxdecor
Decorator< T, ALLOC >::reference_type auxdecor(const std::string &name) const
Fetch an aux decoration, as a non-const reference.
top::PseudoTopReco::SetJetInfo
bool SetJetInfo(const top::Event &event)
Definition: PseudoTopReco.cxx:226
EventTools.h
A few functions for doing operations on particles / events. Currently holds code for dR,...
top::PseudoTopReco::m_W_had
TLorentzVector m_W_had
Definition: PseudoTopReco.h:80
AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore
ServiceHandle< StoreGateSvc > & evtStore()
The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:85
top::PseudoTopReco::m_bJets
std::vector< TLorentzVector > m_bJets
Definition: PseudoTopReco.h:75
xAOD::PseudoTopResult
Interface class.
Definition: PseudoTopResult.h:18
checkTP.save
def save(self, fileName="./columbo.out")
Definition: checkTP.py:178
xAOD::PseudoTopResult::IniVar
void IniVar(bool)
Definition: PseudoTopResult.cxx:14
top::PseudoTopReco::initialize
virtual StatusCode initialize()
Function initialising the tool.
Definition: PseudoTopReco.cxx:30
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
event
POOL::TEvent event(POOL::TEvent::kClassAccess)
top::PseudoTopReco::mWPDG
const double mWPDG
Definition: PseudoTopReco.h:91
top::PseudoTopReco::KinemEdge
const double KinemEdge
Definition: PseudoTopReco.h:90
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
top::PseudoTopReco::SetChargedLeptonInfo
bool SetChargedLeptonInfo(const top::Event &event)
Definition: PseudoTopReco.cxx:196
top::PseudoTopReco::m_W_lep
TLorentzVector m_W_lep
Definition: PseudoTopReco.h:79
xAOD::MissingET_v1::mpx
float mpx() const
Returns .
top::PseudoTopReco::m_lepton
TLorentzVector m_lepton
Definition: PseudoTopReco.h:77
PseudoTopReco.h
top::ParticleLevelEvent
Definition: ParticleLevelEvent.h:24
Event.h
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
top::PseudoTopReco::m_b_had
TLorentzVector m_b_had
Definition: PseudoTopReco.h:82
top::PseudoTopReco::m_leptonType
std::string m_leptonType
Definition: PseudoTopReco.h:73
PathResolver.h
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
top::PseudoTopReco::PseudoTopReco
PseudoTopReco(const std::string &name)
Constructor.
Definition: PseudoTopReco.cxx:19
top::PseudoTopReco::m_b_lep
TLorentzVector m_b_lep
Definition: PseudoTopReco.h:81
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
TopConfig.h
top::PseudoTopReco::finalize
virtual StatusCode finalize()
Function finalizing the tool.
Definition: PseudoTopReco.cxx:354
top::PseudoTopReco::m_top_lep
TLorentzVector m_top_lep
Definition: PseudoTopReco.h:83
top::ParticleLevelEvent::m_electrons
const xAOD::TruthParticleContainer * m_electrons
Pointer to truth level electrons.
Definition: ParticleLevelEvent.h:42
top::PseudoTopReco::m_config
std::shared_ptr< top::TopConfig > m_config
Definition: PseudoTopReco.h:69
a
TList * a
Definition: liststreamerinfos.cxx:10
top::PseudoTopReco::ReconstructLeptonicW
bool ReconstructLeptonicW()
Definition: PseudoTopReco.cxx:314
top::PseudoTopReco::m_top_had
TLorentzVector m_top_had
Definition: PseudoTopReco.h:84
top::Event
Very simple class to hold event data after reading from a file.
Definition: Event.h:49
DataVector::at
const T * at(size_type n) const
Access an element, as an rvalue.
top::PseudoTopReco::m_lightJets
std::vector< TLorentzVector > m_lightJets
Definition: PseudoTopReco.h:76
xAOD::MissingET_v1::mpy
float mpy() const
Returns .
top::ParticleLevelEvent::m_met
const xAOD::MissingET * m_met
Definition: ParticleLevelEvent.h:69
python.compressB64.c
def c
Definition: compressB64.py:93
top::PseudoTopReco::RunReconstruction
bool RunReconstruction()
Definition: PseudoTopReco.cxx:273
top::PseudoTopReco::m_neutrino
TLorentzVector m_neutrino
Definition: PseudoTopReco.h:78
Generated on Tue Jul 2 2024 21:24:53 for ATLAS Offline Software by   doxygen 1.8.18