ATLAS Offline Software
HadronOriginClassifier.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
6 #include "StoreGate/ReadHandle.h"
7 
8 namespace {
10  struct Sample {
12 
14  Sample(int low, int high, GEN_id gen) :
15  low(low), high(high), gen(gen) {}
16 
18  Sample(int id, GEN_id gen) :
19  Sample(id, id, gen) {}
20 
21  int low{};
22  int high{};
23  GEN_id gen{GEN_id::Pythia6};
24  };
25 }
26 
27 namespace DerivationFramework{
28 
29  HadronOriginClassifier::HadronOriginClassifier(const std::string& t, const std::string& n, const IInterface* p):
30  AthAlgTool(t,n,p)
31  {
32  }
33 
35 
37  ATH_MSG_INFO("Initialize " );
38  ATH_MSG_INFO("DSID " << m_DSID );
39 
40  ATH_CHECK(m_mcName.initialize());
41 
42  static const std::vector<Sample> samples = {
43  // all Herwig++/Herwig7 showered samples
44  {346346, 346348, GEN_id::HerwigPP},
45  {410003, GEN_id::HerwigPP}, {410008, GEN_id::HerwigPP}, //aMC@NLO+Hpp
46  {410004, GEN_id::HerwigPP}, {410163, GEN_id::HerwigPP}, //Powheg+Hpp
47  {410232, 410233, GEN_id::HerwigPP}, //first attempt for Powheg+H7 / aMC@NLO+H7
48  {410525, 410530, GEN_id::HerwigPP}, //New Powheg+H7 samples
49  {407037, 407040, GEN_id::HerwigPP}, //Powheg+Hpp MET/HT sliced
50  {410536, 410537, GEN_id::HerwigPP}, {410245, GEN_id::HerwigPP}, //aMC@NLO+H++ , ttbb
51  {410557, 410559, GEN_id::HerwigPP}, // new Powheg+H7, mc16
52  {411082, 411090, GEN_id::HerwigPP}, //Powheg+H7 HF-filtered
53  {407354, 407356, GEN_id::HerwigPP}, //Powheg+H7 ttbar HT-filtered
54  {411233, 411234, GEN_id::HerwigPP}, //Powheg+H7.1.3 ttbar
55  {411316, GEN_id::HerwigPP}, //Powheg+H7 allhad ttbar
56  {411329, 411334, GEN_id::HerwigPP}, //Powheg+H7.1.3 ttbar HF-filtered
57  {411335, 411337, GEN_id::HerwigPP}, //Powheg+H7.1.3 ttbar HT-filtered
58  {412116, 412117, GEN_id::HerwigPP}, //amc@NLO+H7.1.3 ttbar
59  {504329, GEN_id::HerwigPP}, {504333, GEN_id::HerwigPP}, {504341, GEN_id::HerwigPP}, //amc@NLO+H7.2.1 refined ttZ
60  {601239, 601240, GEN_id::HerwigPP},
61  {601668, GEN_id::HerwigPP},
62  {603905, 603906, GEN_id::HerwigPP}, // ttbb Powheg+H7 dilep, ljet, allhad
63  {504337, GEN_id::HerwigPP}, {504345, GEN_id::HerwigPP}, // aMC@NLO+H7 ttZ
64  {526034, GEN_id::HerwigPP}, // aMC@NLO+H7 4tops
65  {600666, 600667, GEN_id::HerwigPP}, // Powheg+H7 ttbar H7UE
66  {601414, 601415, GEN_id::HerwigPP}, // Powheg+H7 ttbar A14
67  {602635, 602635, GEN_id::HerwigPP}, // Powheg+H7 ttH PDF4LHC21
68  {602846, 602849, GEN_id::HerwigPP}, // Powheg+H7 ttW NNPDF30NLO EW
69 
70  // all Pythia8 showered samples
71  {304014, GEN_id::Pythia8}, // amc@NLO+P8 3top
72  {346229, 346234, GEN_id::Pythia8}, // amc@NLO+P8 tHjb
73  {410006, GEN_id::Pythia8}, //Powheg+P8 old main31
74  {410081, GEN_id::Pythia8}, //amc@NLO+P8 ttV
75  {410500, GEN_id::Pythia8}, //Powheg+P8 new main31, hdamp=mt
76  {410501, 410508, GEN_id::Pythia8}, //Powheg+P8 new main31, hdamp=1.5m // Boosted samples are included 410507 410508
77  {410511, 410524, GEN_id::Pythia8}, //Powheg+P8 new main31, hdamp=1.5mt, radiation systematics
78  {410531, 410535, GEN_id::Pythia8}, //Powheg+P8 allhad samples
79  {346343, 346345, GEN_id::Pythia8}, //Powheg+P8 ttH
80  {412123, GEN_id::Pythia8}, // MG+P8 ttW
81  {410155, GEN_id::Pythia8}, // aMC@NlO+P8 ttW
82  {410159, 410160, GEN_id::Pythia8}, //aMC@NLO+P8, old settings
83  {410218, 410220, GEN_id::Pythia8}, // aMC@NlO+P8 ttZ
84  {410276, 410278, GEN_id::Pythia8}, // aMC@NlO+P8 ttZ_lowMass
85  {410225, 410227, GEN_id::Pythia8}, {410274, 410275, GEN_id::Pythia8}, //aMC@NLO+P8, new settings
86  {410568, 410569, GEN_id::Pythia8}, // nonallhad boosted c-filtered
87  {410244, GEN_id::Pythia8}, //aMC@NLO+P8, ttbb (old)
88  {410441, 410442, GEN_id::Pythia8}, //new aMC@NLO+P8 mc16, new shower starting scale
89  {410464, 410466, GEN_id::Pythia8}, //new aMC@NLO+P8 mc16, new shower starting scale, no shower weights
90  {410470, 410472, GEN_id::Pythia8}, {410480, 410482, GEN_id::Pythia8}, //new Powheg+P8 mc16
91  {410452, GEN_id::Pythia8}, //new aMC@NLO+P8 FxFx mc16
92  {411073, 411081, GEN_id::Pythia8}, //Powheg+P8 HF-filtered
93  {412043, 412044, GEN_id::Pythia8}, {500326, GEN_id::Pythia8}, //aMC@NLO+P8 4top
94  {412066, 412074, GEN_id::Pythia8}, //aMC@NLO+P8 HF-filtered
95  {411068, 411070, GEN_id::Pythia8}, //Powheg+P8 ttbb
96  {410265, 410267, GEN_id::Pythia8}, //aMC@NLO+P8 ttbb
97  {411178, 411180, GEN_id::Pythia8}, {411275, GEN_id::Pythia8}, //Powheg+P8 ttbb OTF production - ATLMCPROD-7240
98  {501720, GEN_id::Pythia8}, // aMC@NLO+P8 FxFx ttW
99  {600791, 600792, GEN_id::Pythia8}, //Powheg+P8 ttbb - ATLMCPROD-9179
100  {600737, 600738, GEN_id::Pythia8}, //Powheg+P8 ttbb - ATLMCPROD-9179
101  {601226, 601228, GEN_id::Pythia8}, // Powheg+P8 ttbb bornzerodamp cut 5, ATLMCPROD-9694
102  {407342, 407344, GEN_id::Pythia8}, {411391, GEN_id::Pythia8}, //Powheg+P8 ttbar HT-filtered
103  {407345, 407347, GEN_id::Pythia8}, //Powheg+P8 ttbar MET-filtered
104  {407348, 407350, GEN_id::Pythia8}, //aMC@NLO+P8 ttbar HT-filtered
105  {504330, 504332, GEN_id::Pythia8}, {504334, 504336, GEN_id::Pythia8}, {504338, GEN_id::Pythia8}, {504342, 504344, GEN_id::Pythia8}, {504346, GEN_id::Pythia8}, //aMC@NLO+P8 refined ttZ
106  {601491, 601492, GEN_id::Pythia8}, //Pow+Py8 ttbar pTHard variations - ATLMCPROD-10168
107  {601495, 601498, GEN_id::Pythia8}, //Pow+Py8 ttbar pTHard variations - ATLMCPROD-10168
108  {601229, 601230, GEN_id::Pythia8}, // mc23 ttbar dilep, singlelep
109  {601237, GEN_id::Pythia8}, // mc23 ttbar allhad
110  {601398, 601399, GEN_id::Pythia8}, // mc23 ttbar dilep, singlelep hdamp517p5
111  {601491, GEN_id::Pythia8}, {601495, GEN_id::Pythia8}, {601497, GEN_id::Pythia8}, // mc23 ttbar pThard variations, dilep, singlelep, allhad
112  {601783, 601784, GEN_id::Pythia8}, // Powheg+P8 ttbb bornzerodamp cut 5 pThard variations - ATLMCPROD-10527
113  {603003, 603004, GEN_id::Pythia8}, // Powheg+P8 ttbb nominal and pthard1 allhad
114  {603190, 603193, GEN_id::Pythia8}, // Powheg+P8 ttbb nominal and pthard1 dilep, ljet
115  {604482, 604483, GEN_id::Pythia8}, // mc23 Powheg+P8 ttbar recoilToTop
116  {411369, 411374, GEN_id::Pythia8}, // PP8 ttbar Var1
117  {411290, GEN_id::Pythia8}, // PP8 ttbar rb1p05
118  {508792, GEN_id::Pythia8}, // aMC@NLO+P8 ttbar smeftsim
119  {510203, GEN_id::Pythia8}, // aMC@NLO+P8 4tops
120  {522035, 522038, GEN_id::Pythia8}, // aMC@NLO+P8 ttZqq
121  {523243, GEN_id::Pythia8}, // aMC@NLO+P8 4tops
122  {601356, 601357, GEN_id::Pythia8}, // PP8 ttbar Trec
123  {602067, 602072, GEN_id::Pythia8}, // PP8 ttH pthard1/2
124  {602637, GEN_id::Pythia8}, // PP8 ttH PDF4LHC21
125  {602852, 602853, GEN_id::Pythia8}, // PP8 ttH PDF4LHC21 pthard1
126  {602886, 602889, GEN_id::Pythia8}, // PP8 ttW NNPDF23
127  {603851, 603854, GEN_id::Pythia8}, // PP8 ggH
128  {604224, GEN_id::Pythia8}, // PP8 ttH allhad HTop
129 
130  // all Sherpa showered samples
131  {410186, 410189, GEN_id::Sherpa}, //Sherpa 2.2.0
132  {410249, 410252, GEN_id::Sherpa}, //Sherpa 2.2.1
133  {410342, 410347, GEN_id::Sherpa}, //Sherpa 2.2.1 sys
134  {410350, 410355, GEN_id::Sherpa}, //Sherpa 2.2.1 sys
135  {410357, 410359, GEN_id::Sherpa}, //Sherpa 2.2.1 sys
136  {410361, 410367, GEN_id::Sherpa}, //Sherpa 2.2.1 sys
137  {410281, 410283, GEN_id::Sherpa}, //Sherpa BFilter
138  {410051, GEN_id::Sherpa}, //Sherpa ttbb (ICHEP sample)
139  {410323, 410325, GEN_id::Sherpa}, {410369, GEN_id::Sherpa}, //New Sherpa 2.2.1 ttbb
140  {364345, 364348, GEN_id::Sherpa}, //Sherpa 2.2.4 (test)
141  {410424, 410427, GEN_id::Sherpa}, //Sherpa 2.2.4
142  {410661, 410664, GEN_id::Sherpa}, //Sherpa 2.2.4 ttbb
143  {421152, 421158, GEN_id::Sherpa}, //Sherpa2.2.8 ttbar
144  {413023, GEN_id::Sherpa}, // sherpa 2.2.1 ttZ
145  {700000, GEN_id::Sherpa}, // Sherpa 2.2.8 ttW
146  {700168, GEN_id::Sherpa}, // Sherpa 2.2.10 ttW
147  {700205, GEN_id::Sherpa}, // Sherpa 2.2.10 ttW EWK
148  {700309, GEN_id::Sherpa}, // Sherpa 2.2.11 ttZ
149  {700051, 700054, GEN_id::Sherpa}, //Sherpa2.2.8 ttbb
150  {700121, 700124, GEN_id::Sherpa}, //Sherpa2.2.10 ttbar
151  {700164, 700167, GEN_id::Sherpa}, //Sherpa2.2.10 ttbb
152  {700807, 700809, GEN_id::Sherpa}, //Sherpa2.2.14 ttbar
153  {700659, 700662, GEN_id::Sherpa}, // Sherpa 2.2.12 ttbar maxHTavrgTopPT
154  {700712, GEN_id::Sherpa}, // Sherpa 2.2.14 4tops muQHT2
155  {700986, 700997, GEN_id::Sherpa}, // Sherpa 2.2.14 ttW
156  {701251, 701259, GEN_id::Sherpa}, // Sherpa ttW
157 
158  };
159 
160  // Linear search for sample and assign properties:
161  for (const auto& s : samples) {
162  if (m_DSID>=s.low && m_DSID<=s.high) {
163  m_GenUsed = s.gen;
164  return StatusCode::SUCCESS;
165  }
166  }
167 
168  // the default is Pythia6, so no need to list the Pythia6 showered samples
169  // these are:
170  // 410000-410002
171  // 410007, 410009, 410120-410121
172  // 301528-301532
173  // 303722-303726
174  // 407009-407012
175  // 407029-407036
176  // 410120
177  // 426090-426097
178  // 429007
180 
181  return StatusCode::SUCCESS;
182  }
183 
184  /*
185  --------------------------------------------------------------------------------------------------------------------------------------
186  ------------------------------------------------------------- Hadron Map -------------------------------------------------------------
187  --------------------------------------------------------------------------------------------------------------------------------------
188  */
189 
190  // Define the function GetOriginMap that determines the origin of the hadrons.
191  std::map<const xAOD::TruthParticle*, DerivationFramework::HadronOriginClassifier::HF_id> HadronOriginClassifier::GetOriginMap() const {
192  // Create a set of maps to store the information about the hadrons and the partons
193  std::map<const xAOD::TruthParticle*, int> mainHadronMap; // Map with main hadrons and their flavor.
194  std::map<const xAOD::TruthParticle*, HF_id> partonsOrigin; // Map with partons and their category (from top, W, H, MPI, FSR, extra).
195  std::map<const xAOD::TruthParticle*, const xAOD::TruthParticle*> hadronsPartons; // Map with hadrons and their matched parton.
196  std::map<const xAOD::TruthParticle*, HF_id> hadronsOrigin; // Map with hadrons and their category (from top, W, H, MPI, FSR, extra)
197  // Fill the maps mainHadronMap and partonsOrigin
198  buildPartonsHadronsMaps(mainHadronMap, partonsOrigin);
199  // Create two maps to know which partons and hadrons have already been matched.
200  std::vector<const xAOD::TruthParticle*> matched_partons;
201  std::vector<const xAOD::TruthParticle*> matched_hadrons;
202  // Use a while to go through the HF hadrons in mainHadronMap and partons in partonsOrigin.
203  while (matched_partons.size()<partonsOrigin.size() && matched_hadrons.size()<mainHadronMap.size()){
204  // Create a float variable to store the DeltaR between a parton and the closest hadron.
205  float dR=999.;
206  // Create two pointers for TruthParticle type to go through the partons and hadrons.
207  const xAOD::TruthParticle* hadron=nullptr;
208  const xAOD::TruthParticle* parton=nullptr;
209  // Use a for to go through the partonsOrigin.
210  for(std::map<const xAOD::TruthParticle*, HF_id>::iterator itr = partonsOrigin.begin(); itr!=partonsOrigin.end(); ++itr){
211  // Check if the parton has already been matched to an hadron.
212  if(std::find(matched_partons.begin(), matched_partons.end(), (*itr).first) != matched_partons.end()) continue;
213  // Extract the pt of the parton.
214  TVector3 v, vtmp;
215  if ((*itr).first->pt()>0.)
216  v.SetPtEtaPhi((*itr).first->pt(),(*itr).first->eta(),(*itr).first->phi());
217  else // Protection against FPE from eta and phi calculation
218  v.SetXYZ(0.,0.,(*itr).first->pz());
219  // Use a for to go through the HF hadrons in mainHadronMap.
220  for(std::map<const xAOD::TruthParticle*, int>::iterator it = mainHadronMap.begin(); it!=mainHadronMap.end(); ++it){
221  // Check if the hadron has already been matched to a parton.
222  if(std::find(matched_hadrons.begin(), matched_hadrons.end(), (*it).first) != matched_hadrons.end()) continue;
223  // Check if the hadron's flavour matches the one of the parton.
224  if((*it).second != (*itr).first->absPdgId()) continue;
225  // Extract the pt of the hadron.
226  vtmp.SetPtEtaPhi((*it).first->pt(),(*it).first->eta(),(*it).first->phi());
227  // Compute Delta R between hadron and parton and store in dR if it is smaller than the current value.
228  // Also store the parton and hadron in the pointers that have been previous created.
229  if(vtmp.DeltaR(v) < dR){
230  dR = vtmp.DeltaR(v);
231  hadron = (*it).first;
232  parton = (*itr).first;
233  }
234  }//loop hadrons
235  }//loop partons
236  // Add the matched part-hadron pair in the corresponding maps.
237  matched_partons.push_back(parton);
238  matched_hadrons.push_back(hadron);
239  hadronsPartons[ hadron ] = parton;
240  }
241 
242  // Use a for to go through the HF hadrons in mainHadronMap.
243  for(std::map<const xAOD::TruthParticle*, int>::iterator it = mainHadronMap.begin(); it!=mainHadronMap.end(); ++it){
244  // Extract the current hadron.
245  const xAOD::TruthParticle* hadron = (*it).first;
246  // Check if the hadron has been matched to a parton.
247  // If it has been matched to any hadron, use it to determine the origin.
248  // Otherwise, the hadron is considered extra.
249  if(hadronsPartons.find(hadron)!=hadronsPartons.end()){
250  hadronsOrigin[hadron] = partonsOrigin[ hadronsPartons[hadron] ];
251  } else{
252  hadronsOrigin[hadron] = extrajet;
253  }
254  }
255  return hadronsOrigin;
256  }
257 
258  // Define the function buildPartonsHadronsMaps that determines the flavour of the hadrons and the origin of the partons.
259  void HadronOriginClassifier::buildPartonsHadronsMaps(std::map<const xAOD::TruthParticle*,int>& mainHadronMap, std::map<const xAOD::TruthParticle*,HF_id>& partonsOrigin) const {
260  // Extract the TruthParticles container.
261  const EventContext& ctx = Gaudi::Hive::currentContext();
262  SG::ReadHandle<xAOD::TruthEventContainer> xTruthEventContainer(m_mcName, ctx);
263  if (!xTruthEventContainer.isValid()) {
264  ATH_MSG_WARNING("Could not retrieve " <<m_mcName);
265  }
266 
267  // Create a container with TruthParticles to store the hadrons that has already been saved.
268  std::set<const xAOD::TruthParticle*> usedHadron;
269  for ( const auto* truthevent : *xTruthEventContainer ) {
270  // Use a for to go through the TruthParticles.
271  for(unsigned int i = 0; i < truthevent->nTruthParticles(); i++){
272  // Extract the i-th particle.
273  const xAOD::TruthParticle* part = truthevent->truthParticle(i);
274  if(!part) continue;
275  // Simulated particles are not considered.
277  // Create a set of boolean variables to indicate the type of particle.
278  bool isbquark = false; // The particle is a b-quark.
279  bool iscquark = false; // The particle is a c-quark.
280  bool isHFhadron = false; // The particle is a HF hadron.
281  // Extract the pdgid of the particle and use it to determine the type of particle.
282  int pdgid = part->absPdgId();
283  if( MC::isBottom(pdgid) ){
284  isbquark=true;
285  }
286  else if( MC::isCharm(pdgid) ){
287  iscquark=true;
288  }
290  isHFhadron=true;
291  }
292  else{
293  continue;
294  }
295  // For HF quarks (b or c), check their category.
296  // The category is determined looking for the parents.
297  if(isbquark){
298  // In this case, the parton is a b-quark.
299  // Check the category of the b-quark.
301  partonsOrigin[ part ] = b_from_W;
302  }
303  else if(isDirectlyFromTop(part)){
304  partonsOrigin[ part ] = b_from_top;
305  }
306  else if((IsHerwigPP()||IsSherpa())&&isDirectlyFSR(part)){
307  partonsOrigin[ part ] = b_FSR;
308  }
309  else if(IsPythia8()&&isDirectlyFSRPythia8(part)){
310  partonsOrigin[ part ] = b_FSR;
311  }
312  else if(IsPythia6()&&isDirectlyFSRPythia6(part)){
313  partonsOrigin[ part ] = b_FSR;
314  }
315  else if(IsPythia6()&&isDirectlyMPIPythia6(part)){
316  partonsOrigin[ part ] = b_MPI;
317  }
318  else if(IsPythia8()&&isDirectlyMPIPythia8(part)){
319  partonsOrigin[ part ] = b_MPI;
320  }
321  else if(IsSherpa()&&isDirectlyMPISherpa(part)){
322  partonsOrigin[ part ] = b_MPI;
323  }
324  }
325  if(iscquark){
326  // In this case, the parton is a c-quark.
327  // Check the category of the b-quark.
329  partonsOrigin[ part ] = c_from_W;
330  }
331  else if(isDirectlyFromTop(part)){
332  partonsOrigin[ part ] = c_from_top;
333  }
334  else if((IsHerwigPP()&&IsSherpa())&&isDirectlyFSR(part)){
335  partonsOrigin[ part ] = c_FSR;
336  }
337  else if(IsPythia8()&&isDirectlyFSRPythia8(part)){
338  partonsOrigin[ part ] = c_FSR;
339  }
340  else if(IsPythia6()&&isDirectlyFSRPythia6(part)){
341  partonsOrigin[ part ] = c_FSR;
342  }
343  else if(IsPythia6()&&isDirectlyMPIPythia6(part)){
344  partonsOrigin[ part ] = c_MPI;
345  }
346  else if(IsPythia8()&&isDirectlyMPIPythia8(part)){
347  partonsOrigin[ part ] = c_MPI;
348  }
349  else if(IsSherpa()&&isDirectlyMPISherpa(part)){
350  partonsOrigin[ part ] = c_MPI;
351  }
352  }
353  // The HF hadrons are stored in the map mainHadronMap if they are not repeated.
354  if(isHFhadron && !isCHadronFromB(part)){
355  // In this case, the particle is a HF hadron but not a C-Hadron from a B-hadron.
356  // If the hadron is not in usedHadron, then add it in mainHadronMap with fillHadronMap function.
357  if(usedHadron.insert(part).second) {
358  fillHadronMap(usedHadron, mainHadronMap,part,part);
359  }
360  }
361  }//loop on particles
362  }//loop on truthevent container
363  }
364 
365  /*
366  ---------------------------------------------------------------------------------------------------------------------------------------
367  ------------------------------------------------------------ Particle Type ------------------------------------------------------------
368  ---------------------------------------------------------------------------------------------------------------------------------------
369  */
370  bool HadronOriginClassifier::isCHadronFromB(const xAOD::TruthParticle* part, std::shared_ptr<std::set<const xAOD::TruthParticle*>> checked ) const{
371  if(!MC::isCharmHadron(part)) return false;
372  if (!checked) checked = std::make_shared<std::set<const xAOD::TruthParticle*>>();
373  checked ->insert(part);
374 
375  for(unsigned int i=0; i<part->nParents(); ++i){
376  const xAOD::TruthParticle* parent = part->parent(i);
377  if(!parent) continue;
378  if(checked->count(parent)) continue;
379  checked->insert(parent);
380  if( MC::isBottomHadron(parent) ){
381  return true;
382  }
384  if(isCHadronFromB(parent))return true;
385  }
386  }
387 
388  return false;
389  }
390 
391  // Define the function fillHadronMap that fills the map of hadrons with their flavour.
392  void HadronOriginClassifier::fillHadronMap(std::set<const xAOD::TruthParticle*>& usedHadron, std::map<const xAOD::TruthParticle*,int>& mainHadronMap, const xAOD::TruthParticle* mainhad, const xAOD::TruthParticle* ihad, bool decayed) const {
393  // Fist, check that the consdired hadron has a non-null pointer
394  if (!ihad) return;
395  usedHadron.insert(ihad);
396  // Create two variables to indicate the flavour of the parents and childrens particles that will be considered.
397  // Create a boolean to indicate if the particles considered are from the final state.
398  int parent_flav,child_flav;
399  bool isFinal = true;
400  // Check if the considered hadron has children.
401  if(!ihad->nChildren()) return;
402  // Use a for to go through the children.
403  for(unsigned int j=0; j<ihad->nChildren(); ++j){
404  // Extract the j-th children.
405  const xAOD::TruthParticle* child = ihad->child(j);
406  if(!child) continue;
407  if(decayed){
408  fillHadronMap(usedHadron, mainHadronMap,mainhad,child,true);
409  isFinal=false;
410  }
411  else{
412  child_flav = std::abs(MC::leadingQuark(child));
413  if(child_flav!=4 && child_flav!=5) continue;
414  parent_flav = std::abs(MC::leadingQuark(mainhad));
415  if(child_flav!=parent_flav) continue;
416  fillHadronMap(usedHadron, mainHadronMap,mainhad,child);
417  isFinal=false;
418  }
419  }
420 
421  if(isFinal && !decayed){
422  mainHadronMap[mainhad]=std::abs(MC::leadingQuark(mainhad));
423  for(unsigned int j=0; j<ihad->nChildren(); ++j){
424  const xAOD::TruthParticle* child = ihad->child(j);
425  if(!child) continue;
426  fillHadronMap(usedHadron, mainHadronMap,mainhad,child,true);
427  }
428  }
429  }
430 
431  /*
432  ---------------------------------------------------------------------------------------------------------------------------------------
433  ----------------------------------------------------------- Particle Origin -----------------------------------------------------------
434  ---------------------------------------------------------------------------------------------------------------------------------------
435  */
436 
437  // Define the function isFromTop that indicates if a particle comes from top.
438 
440  // Find the first parent of the considered particle that is different from the particle.
441  const xAOD::TruthParticle* initpart = findInitial(part);
442  // Check if this parent comes from the top with function isDirectlyFromTop.
443  return isDirectlyFromTop(initpart);
444  }
445 
446  // Define the function isDirectlyFromTop that indicates if a particle comes from the direct decay of top.
448  // First, make sure the consdired particle has a non-null pointer and it has parents.
449  // Otherwise, return false.
450  if(!part || !part->nParents()) return false;
451  // Go through the parents of the particle.
452  for(unsigned int i=0; i<part->nParents(); ++i){
453  // Extract the i-th parent.
454  const xAOD::TruthParticle* parent = part->parent(i);
455  if(!parent) continue;
456  // If the i-th parent is a top, then return true
457  if( MC::isTop(parent) ) return true;
458  }
459  // If a top is no the parent, then return false.
460  return false;
461  }
462 
463  // Define the function isFromWTop that indicates if a particle comes from the decay chain t->Wb.
464 
466  // Find the first parent of the considered particle that is different from the particle.
467  const xAOD::TruthParticle* initpart = findInitial(part);
468  return isDirectlyFromWTop(initpart);
469  }
470 
471  // Define the function isDirectlyFromWTop that indicates if a particle comes from the direct decay of a W from a top.
473  // First, make sure the consdired particle has a non-null pointer and it has parents.
474  // Otherwise, return false.
475  if(!part || !part->nParents()) return false;
476  // Use a for to go though the parents.
477  for(unsigned int i=0; i<part->nParents(); ++i){
478  // Get the i-th parent.
479  const xAOD::TruthParticle* parent = part->parent(i);
480  if(!parent) continue;
481  if( MC::isW(parent)){
482  if( isFromTop(parent) ) return true;
483  }
484  }
485  // In this case, none of the parents of the particle is a W from top.
486  // Hence, return false.
487  return false;
488  }
489 
491  if(!part->nParents()) return false;
492  for(unsigned int i=0; i<part->nParents(); ++i){
493  const xAOD::TruthParticle* parent = part->parent(i);
494  if(!parent) continue;
495  if( MC::isPhoton(parent) || parent->absPdgId()<MC::BQUARK ) return true;
496  }
497  return false;
498  }
499 
501  const xAOD::TruthParticle* initpart = findInitial(part);
502  return isDirectlyFromGluonQuark(initpart);
503  }
504 
506  if(!part->nParents()) return false;
507  for(unsigned int i=0; i<part->nParents(); ++i){
508  const xAOD::TruthParticle* parent = part->parent(i);
509  if(!parent) continue;
510  if(!MC::isW(parent)) continue;
511  if(MC::isCharm(part)){
512  //trick to get at least 50% of PowhegPythia c from FSR
513  if(part->pdgId()==-(parent->pdgId())/6){
514  if( isFromGluonQuark(parent) ) return true;
515  }
516  }
517  else{
518  if( isFromGluonQuark(parent) ) return true;
519  }
520  }
521  return false;
522  }
523 
525  if(!part->nParents()) return false;
526  for(unsigned int i=0; i<part->nParents(); ++i){
527  const xAOD::TruthParticle* parent = part->parent(i);
528  if(!parent) continue;
530  if( isFromQuarkTop( parent ) ) return true;
531  }
532  }
533  return false;
534  }
535 
537  if(!part->nParents()) return false;
538  for(unsigned int i=0; i<part->nParents(); ++i){
539  const xAOD::TruthParticle* parent = part->parent(i);
540  if(!parent) continue;
541  if( parent->absPdgId()<MC::TQUARK ) {
542  if(isFromTop(parent)){
543  return true;
544  }
545  else if(isFromWTop(parent)){
546  return true;
547  }
548  }
549  }
550 
551  return false;
552  }
553 
555  const xAOD::TruthParticle* initpart = findInitial(part);
556  return isDirectlyFromQuarkTop(initpart);
557  }
558 
559  // Define the function isDirectlyFSRPythia8 that indicates if a particle comes from Final State Radiation in samples generated with Pythia8.
560 
562  // First, check if the particle has parents and return false if it does not.
563  if(!part->nParents()) return false;
564  // Use a for to go through the parents.
565  for(unsigned int i=0; i<part->nParents(); ++i){
566 
567  // Extract the i-th parent.
568 
569  const xAOD::TruthParticle* parent = part->parent(i);
570  if(!parent) continue;
572  if( isFromQuarkTopPythia8( parent ) ) return true;
573  }
574  }
575  // In this case, no parent from the particle is a gluon or a photon coming from a top
576  // Hence, the particle is not from FSR and false is not returned.
577  return false;
578  }
579 
580  // Define the function isDirectlyFromQuarkTopPythia8 that indicates if a particle comes from direct decay of the top in samples generated with Pythia8.
582  // First, make sure the consdired particle has a non-null pointer and it has parents.
583  // Otherwise, return false.
584  if(!part->nParents()) return false;
585  // Use a for to go through the parents.
586  for(unsigned int i=0; i<part->nParents(); ++i){
587  // Extract the i-th parent.
588  const xAOD::TruthParticle* parent = part->parent(i);
589  if(!parent) continue;
590  // Check if the parent is a quark different from the top.
591  if( parent->absPdgId()<MC::TQUARK ) {
592  // In this case, the parent is a quark different from top.
593  // Check if it comes from the decay chain of the t->Wb.
594  // If it is the case, return true.
595  if(isFromWTop(parent)){
596  return true;
597  }
598  }
599  }
600  // In this case, any of the parents of the particle comes from t->Wb chaing.
601  // Hence, the particle does not come from the top directly and false is returned.
602  return false;
603  }
604 
605  // Define the function isFromQuarkTopPythia8 that indicates if a particle comes from top in samples generated with Pythia8.
607  // Find the first parent of the considered particle that is different from the particle.
608  const xAOD::TruthParticle* initpart = findInitial(part);
609  // Check if this parent comes from the top with function isDirectlyFromQuarkTopPythia8.
610  return isDirectlyFromQuarkTopPythia8(initpart);
611  }
612 
614  if(!part->nParents()) return false;
615  for(unsigned int i=0; i<part->nParents(); ++i){
616  const xAOD::TruthParticle* parent = part->parent(i);
617  if(!parent) continue;
618  if( parent->absPdgId() == MC::PROTON && MC::isPhysical(part) ) return true;
619  }
620  return false;
621  }
622 
624  const xAOD::TruthParticle* initpart = findInitial(part);
625  return MC::Pythia8::isConditionC(initpart);
626  }
627 
629  if(!part->hasProdVtx()) return false;
630  const xAOD::TruthVertex* vertex = part->prodVtx();
631  return HepMC::status(vertex) == 2;
632  }
633 
634  /*
635  --------------------------------------------------------------------------------------------------------------------------------------
636  ---------------------------------------------------------- Particle Parents ----------------------------------------------------------
637  --------------------------------------------------------------------------------------------------------------------------------------
638  */
639 
640  // Define the function findInitial which finds the first parent of a particle that is not the particle itself.
641  const xAOD::TruthParticle* HadronOriginClassifier::findInitial(const xAOD::TruthParticle* part, std::shared_ptr<std::set<const xAOD::TruthParticle*>> checked ) const{
642  // If the particle has no parent, return the particle.
643  if(!part->nParents()) return part;
644  if (!checked) checked = std::make_shared<std::set<const xAOD::TruthParticle*>>();
645  // Use a for to go through the parents.
646  for(unsigned int i=0; i<part->nParents(); ++i){
647  // Extract the i-th parent.
648  const xAOD::TruthParticle* parent = part->parent(i);
649  if(!parent) continue;
650  if(checked->count(parent)) continue;
651  checked->insert(parent);
652  // If the parent has the same pdgId as the particle, then it means that the parent is the same as the considered particle.
653  // This happens if the particle irradiates for example.
654  // In this case, try to look for the first parent of i-th parent that is being considered.
655  // Repeat the process until you find a particle different from the considred one or that has no parent.
656 
657  if( part->pdgId() == parent->pdgId() ){
658  return findInitial(parent, checked);
659  }
660  }
661  // In this case, no parent different from the considered particle has been found.
662  // Hence, return the particle.
663  return part;
664  }
665 
666 }//namespace
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
DerivationFramework::HadronOriginClassifier::isDirectlyMPIPythia8
bool isDirectlyMPIPythia8(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:623
DerivationFramework::HadronOriginClassifier::m_DSID
Gaudi::Property< int > m_DSID
Definition: HadronOriginClassifier.h:90
MCTruthPartClassifier::hadron
@ hadron
Definition: TruthClassifiers.h:148
DerivationFramework::HadronOriginClassifier::isDirectlyMPISherpa
static bool isDirectlyMPISherpa(const xAOD::TruthParticle *part)
Definition: HadronOriginClassifier.cxx:628
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
DerivationFramework::HadronOriginClassifier::GEN_id::HerwigPP
@ HerwigPP
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:67
DerivationFramework::HadronOriginClassifier::b_FSR
@ b_FSR
Definition: HadronOriginClassifier.h:39
HadronOriginClassifier.h
DerivationFramework::HadronOriginClassifier::isFromQuarkTop
bool isFromQuarkTop(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:554
DerivationFramework::HadronOriginClassifier::GEN_id::Sherpa
@ Sherpa
MC::Pythia8::isConditionC
bool isConditionC(const T &p)
Definition: HepMCHelpers.h:27
DerivationFramework::HadronOriginClassifier::c_FSR
@ c_FSR
Definition: HadronOriginClassifier.h:39
DerivationFramework::HadronOriginClassifier::IsHerwigPP
bool IsHerwigPP() const
Definition: HadronOriginClassifier.h:82
skel.it
it
Definition: skel.GENtoEVGEN.py:407
DerivationFramework::HadronOriginClassifier::isFromQuarkTopPythia8
bool isFromQuarkTopPythia8(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:606
DerivationFramework::HadronOriginClassifier::GetOriginMap
std::map< const xAOD::TruthParticle *, HF_id > GetOriginMap() const
Definition: HadronOriginClassifier.cxx:191
DerivationFramework::HadronOriginClassifier::c_MPI
@ c_MPI
Definition: HadronOriginClassifier.h:38
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
DerivationFramework::HadronOriginClassifier::c_from_W
@ c_from_W
Definition: HadronOriginClassifier.h:40
DerivationFramework::HadronOriginClassifier::isFromTop
bool isFromTop(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:439
isBottomHadron
bool isBottomHadron(const T &p)
Definition: AtlasPID.h:901
master.gen
gen
Definition: master.py:32
DerivationFramework::HadronOriginClassifier::c_from_top
@ c_from_top
Definition: HadronOriginClassifier.h:41
isGluon
bool isGluon(const T &p)
Definition: AtlasPID.h:370
DerivationFramework::HadronOriginClassifier::m_mcName
SG::ReadHandleKey< xAOD::TruthEventContainer > m_mcName
Definition: HadronOriginClassifier.h:87
DerivationFramework::HadronOriginClassifier::isDirectlyFSR
bool isDirectlyFSR(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:524
DerivationFramework::HadronOriginClassifier::b_from_W
@ b_from_W
Definition: HadronOriginClassifier.h:40
DerivationFramework::HadronOriginClassifier::fillHadronMap
void fillHadronMap(std::set< const xAOD::TruthParticle * > &usedHadron, std::map< const xAOD::TruthParticle *, int > &mainHadronMap, const xAOD::TruthParticle *mainhad, const xAOD::TruthParticle *ihad, bool decayed=false) const
Definition: HadronOriginClassifier.cxx:392
DerivationFramework::HadronOriginClassifier::findInitial
const xAOD::TruthParticle * findInitial(const xAOD::TruthParticle *part, std::shared_ptr< std::set< const xAOD::TruthParticle * >> checked=nullptr) const
Definition: HadronOriginClassifier.cxx:641
DerivationFramework::HadronOriginClassifier::~HadronOriginClassifier
virtual ~HadronOriginClassifier()
Definition: HadronOriginClassifier.cxx:34
MC::isPhysical
bool isPhysical(const T &p)
Identify if the particle is physical, i.e. is stable or decayed.
Definition: HepMCHelpers.h:51
DerivationFramework::HadronOriginClassifier::IsSherpa
bool IsSherpa() const
Definition: HadronOriginClassifier.h:85
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
isBottom
bool isBottom(const T &p)
Definition: AtlasPID.h:179
HepMC::is_simulation_particle
bool is_simulation_particle(const T &p)
Method to establish if a particle (or barcode) was created during the simulation (TODO update to be s...
Definition: MagicNumbers.h:354
DerivationFramework::HadronOriginClassifier::b_MPI
@ b_MPI
Definition: HadronOriginClassifier.h:38
lumiFormat.i
int i
Definition: lumiFormat.py:85
DerivationFramework::HadronOriginClassifier::isFromGluonQuark
bool isFromGluonQuark(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:500
leadingQuark
int leadingQuark(const T &p)
Definition: AtlasPID.h:878
DerivationFramework::HadronOriginClassifier::isDirectlyFromQuarkTop
bool isDirectlyFromQuarkTop(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:536
beamspotman.n
n
Definition: beamspotman.py:729
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
DerivationFramework::HadronOriginClassifier::isFromWTop
bool isFromWTop(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:465
DerivationFramework::HadronOriginClassifier::GEN_id
GEN_id
Definition: HadronOriginClassifier.h:44
xAOD::TruthParticle_v1
Class describing a truth particle in the MC record.
Definition: TruthParticle_v1.h:37
DerivationFramework::HadronOriginClassifier::isDirectlyFromWTop
bool isDirectlyFromWTop(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:472
DerivationFramework::HadronOriginClassifier::isDirectlyFromGluonQuark
static bool isDirectlyFromGluonQuark(const xAOD::TruthParticle *part)
Definition: HadronOriginClassifier.cxx:490
DerivationFramework::HadronOriginClassifier::m_GenUsed
GEN_id m_GenUsed
Definition: HadronOriginClassifier.h:91
test_pyathena.parent
parent
Definition: test_pyathena.py:15
DerivationFramework::HadronOriginClassifier::GEN_id::Pythia8
@ Pythia8
DerivationFramework::HadronOriginClassifier::isDirectlyFSRPythia6
bool isDirectlyFSRPythia6(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:505
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
DerivationFramework
THE reconstruction tool.
Definition: ParticleSortingAlg.h:24
DerivationFramework::HadronOriginClassifier::HadronOriginClassifier
HadronOriginClassifier(const std::string &t, const std::string &n, const IInterface *p)
Definition: HadronOriginClassifier.cxx:29
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
xAOD::TruthVertex_v1
Class describing a truth vertex in the MC record.
Definition: TruthVertex_v1.h:37
xAOD::TruthParticle_v1::nChildren
size_t nChildren() const
Number of children of this particle.
Definition: TruthParticle_v1.cxx:135
id
SG::auxid_t id
Definition: Control/AthContainers/Root/debug.cxx:239
DerivationFramework::HadronOriginClassifier::IsPythia6
bool IsPythia6() const
Definition: HadronOriginClassifier.h:84
xAOD::TruthParticle_v1::child
const TruthParticle_v1 * child(size_t i) const
Retrieve the i-th mother (TruthParticle) of this TruthParticle.
Definition: TruthParticle_v1.cxx:144
DerivationFramework::HadronOriginClassifier::isDirectlyFSRPythia8
bool isDirectlyFSRPythia8(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:561
DerivationFramework::HadronOriginClassifier::isDirectlyMPIPythia6
static bool isDirectlyMPIPythia6(const xAOD::TruthParticle *part)
Definition: HadronOriginClassifier.cxx:613
python.PyAthena.v
v
Definition: PyAthena.py:154
DerivationFramework::HadronOriginClassifier::IsPythia8
bool IsPythia8() const
Definition: HadronOriginClassifier.h:83
isW
bool isW(const T &p)
Definition: AtlasPID.h:379
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
isTop
bool isTop(const T &p)
Definition: AtlasPID.h:182
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
xAOD::EgammaHelpers::isPhoton
bool isPhoton(const xAOD::Egamma *eg)
is the object a photon
Definition: EgammaxAODHelpers.cxx:21
isCharm
bool isCharm(const T &p)
Definition: AtlasPID.h:176
DerivationFramework::HadronOriginClassifier::isDirectlyFromTop
static bool isDirectlyFromTop(const xAOD::TruthParticle *part)
Definition: HadronOriginClassifier.cxx:447
isCharmHadron
bool isCharmHadron(const T &p)
Definition: AtlasPID.h:900
python.SystemOfUnits.s
float s
Definition: SystemOfUnits.py:147
DerivationFramework::HadronOriginClassifier::isDirectlyFromQuarkTopPythia8
bool isDirectlyFromQuarkTopPythia8(const xAOD::TruthParticle *part) const
Definition: HadronOriginClassifier.cxx:581
DerivationFramework::HadronOriginClassifier::extrajet
@ extrajet
Definition: HadronOriginClassifier.h:37
ReadHandle.h
Handle class for reading from StoreGate.
AthAlgTool
Definition: AthAlgTool.h:26
HepMC::status
int status(const T &p)
Definition: MagicNumbers.h:143
DerivationFramework::HadronOriginClassifier::buildPartonsHadronsMaps
void buildPartonsHadronsMaps(std::map< const xAOD::TruthParticle *, int > &mainHadronMap, std::map< const xAOD::TruthParticle *, HF_id > &partonsOrigin) const
Definition: HadronOriginClassifier.cxx:259
HepMCHelpers.h
DerivationFramework::HadronOriginClassifier::b_from_top
@ b_from_top
Definition: HadronOriginClassifier.h:41
DerivationFramework::HadronOriginClassifier::isCHadronFromB
bool isCHadronFromB(const xAOD::TruthParticle *part, std::shared_ptr< std::set< const xAOD::TruthParticle * >> checked=nullptr) const
Definition: HadronOriginClassifier.cxx:370
DerivationFramework::HadronOriginClassifier::GEN_id::Pythia6
@ Pythia6
DerivationFramework::HadronOriginClassifier::initialize
virtual StatusCode initialize() override
Definition: HadronOriginClassifier.cxx:36