Loading [MathJax]/extensions/tex2jax.js
ATLAS Offline Software
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
MCTruthClassifierGen.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /*
6  * Implementation file that mainly contains the code logic
7  * dealing with Truth - record classification
8  * Contributors: Pierre-Antoine Delsart
9  * Andrii Verbytskyi <andrii.verbytskyi@mpp.mpg.de>
10  */
11 
15 #ifndef XAOD_ANALYSIS
16 #include "AtlasHepMC/GenEvent.h"
17 #include "AtlasHepMC/GenVertex.h"
18 #include "AtlasHepMC/GenParticle.h"
19 #endif
21 using namespace MCTruthPartClassifier;
22 using std::abs;
23 
24 #ifndef XAOD_ANALYSIS
25 std::pair<ParticleType, ParticleOrigin>
27  // Retrieve the links between HepMC and xAOD::TruthParticle
28  const EventContext& ctx = info ? info->eventContext : Gaudi::Hive::currentContext();
29  SG::ReadHandle<xAODTruthParticleLinkVector> truthParticleLinkVecReadHandle(m_truthLinkVecReadHandleKey, ctx);
30  if (!truthParticleLinkVecReadHandle.isValid()) {
31  ATH_MSG_WARNING(" Invalid ReadHandle for xAODTruthParticleLinkVector with key: " << truthParticleLinkVecReadHandle.key());
32  return std::make_pair(Unknown, NonDefined);
33  }
34  ElementLink<xAOD::TruthParticleContainer> tplink = truthParticleLinkVecReadHandle->find (theLink);
35  if (tplink.isValid()) {
36  return particleTruthClassifier (*tplink, info);
37  }
38  return std::make_pair(Unknown, NonDefined);
39 }
40 
41 std::pair<ParticleType, ParticleOrigin>
43  ParticleType partType = Unknown;
44  ParticleOrigin partOrig = NonDefined;
45 
46  if (!theGenPart) return std::make_pair(partType, partOrig);
47 
48  // Retrieve the links between HepMC and xAOD::TruthParticle
49  const EventContext& ctx = info ? info->eventContext : Gaudi::Hive::currentContext();
50 
51  SG::ReadHandle<xAODTruthParticleLinkVector> truthParticleLinkVecReadHandle(m_truthLinkVecReadHandleKey, ctx);
52  if (!truthParticleLinkVecReadHandle.isValid()) {
53  ATH_MSG_WARNING( " Invalid ReadHandle for xAODTruthParticleLinkVector with key: " << truthParticleLinkVecReadHandle.key());
54  return std::make_pair(partType, partOrig);
55  }
56  for (const auto *const entry : *truthParticleLinkVecReadHandle) {
57  if (entry->first.isValid() && entry->second.isValid() && HepMC::is_same_particle(entry->first,theGenPart)) {
58  const xAOD::TruthParticle* truthParticle = *entry->second;
59  if (!theGenPart || !truthParticle ||
60  theGenPart->pdg_id() != truthParticle->pdgId() ||
61  theGenPart->status() != truthParticle->status() ||
62  HepMC::is_same_particle(theGenPart,truthParticle)) {
64  "HepMC::GenParticle and xAOD::TruthParticle do not match");
65  return std::make_pair(partType, partOrig);
66  }
67  return particleTruthClassifier(truthParticle, info);
68  }
69  }
70  return std::make_pair(partType, partOrig);
71 }
72 #endif
73 
74 std::pair<ParticleType, ParticleOrigin>
77  ATH_MSG_DEBUG("Executing particleTruthClassifier");
78 
79  ParticleType partType = Unknown;
80  ParticleOrigin partOrig = NonDefined;
81  if (!thePart){
82  return std::make_pair(partType, partOrig);
83  }
84 
85  const EventContext& ctx = info ? info->eventContext : Gaudi::Hive::currentContext();
87  if (!info) { info = &tmpinfo; }
88  info->genPart = thePart;
89 
90  // retrieve collection and get a pointer
91  SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainerReadHandle(m_truthParticleContainerKey,ctx);
92  if (!truthParticleContainerReadHandle.isValid()) {
93  ATH_MSG_WARNING( " Invalid ReadHandle for xAOD::TruthParticleContainer with key: " << truthParticleContainerReadHandle.key());
94  return std::make_pair(partType, partOrig);
95  }
96 
97  ATH_MSG_DEBUG("xAODTruthParticleContainer with key " << truthParticleContainerReadHandle.key() << " has valid ReadHandle ");
98 
99  if (!MC::isStable(thePart) && !MC::isDecayed(thePart)) {
100  return std::make_pair(GenParticle, partOrig);
101  }
102  bool isPartHadr = MC::isHadron(thePart)&&!MC::isBeam(thePart);
103  if (MC::isDecayed(thePart) && (!MC::isTau(thePart) && !isPartHadr)) return std::make_pair(GenParticle, partOrig);
104 
105  // SUSY datasets: tau(satus==2)->tau(satus==2)
106  if (MC::isDecayed(thePart) && MC::isTau(thePart)) {
107  const xAOD::TruthVertex* endVert = thePart->decayVtx();
108  if (endVert != nullptr) {
109  int numOfDaught = endVert->nOutgoingParticles();
110  if (numOfDaught == 1 && MC::isTau(endVert->outgoingParticle(0))) {
111  return std::make_pair(GenParticle, partOrig);
112  }
113  }
114  }
115 
116  if (MC::isStable(thePart) && MC::isSUSY(thePart)) return std::make_pair(SUSYParticle, partOrig);
117 
118  if (MC::isStable(thePart) && MC::isBSM(thePart)) return std::make_pair(OtherBSMParticle, partOrig);
119 
120  if (MC::isDecayed(thePart) &&
121  (!MC::isElectron(thePart) && !MC::isMuon(thePart) &&
122  !MC::isTau(thePart) && !MC::isPhoton(thePart)) &&
123  !isPartHadr)
124  return std::make_pair(GenParticle, partOrig);
125 
126  if (abs(thePart->pdg_id()) > 1000000000) return std::make_pair(NuclFrag, partOrig);
127 
128  if ( !MC::isSMLepton(thePart) && !MC::isPhoton(thePart) && !isPartHadr) return std::make_pair(partType, partOrig);
129  // don't consider generator particles
130 
131  const xAOD::TruthVertex* partOriVert = thePart->hasProdVtx() ? thePart->prodVtx() : nullptr;
132 
133  const xAOD::TruthParticle* theMoth{};
134  if (partOriVert != nullptr) {
135  for (const auto& temp: partOriVert->particles_in()) {if (temp) theMoth = temp;}
136  }
137  int motherPDG = theMoth?theMoth->pdg_id():0;
138  info->setMotherProperties(theMoth);
139 
140  if (!partOriVert && HepMC::is_simulation_particle(thePart)) {
141  return std::make_pair(NonPrimary, partOrig);
142  }
143  if (!partOriVert && MC::isElectron(thePart)) {
144  // to define electron out come status
145  bool isPrompt = false;
146  partOrig = defOrigOfElectron(truthParticleContainerReadHandle.ptr(), thePart, isPrompt, info);
147  return std::make_pair(UnknownElectron, partOrig);
148  }
149  if (!partOriVert && MC::isMuon(thePart)) {
150  // to define electron out come status
151  bool isPrompt = false;
152  partOrig = defOrigOfMuon(truthParticleContainerReadHandle.ptr(), thePart, isPrompt, info);
153  return std::make_pair(UnknownMuon, partOrig);
154  }
155  if (!partOriVert && MC::isTau(thePart)) {
156  // to define electron out come status
157  partOrig = defOrigOfTau(truthParticleContainerReadHandle.ptr(), thePart, motherPDG, info);
158  return std::make_pair(UnknownTau, partOrig);
159  }
160  if (!partOriVert && MC::isPhoton(thePart)) {
161  // to define photon out come
162  bool isPrompt = false;
163  partOrig = defOrigOfPhoton(truthParticleContainerReadHandle.ptr(), thePart, isPrompt, info);
164  return std::make_pair(UnknownPhoton, partOrig);
165  }
166  if (!partOriVert && MC::isNeutrino(thePart)) {
167  // to define neutrino outcome
168  info->particleOutCome = NonInteract;
169  return std::make_pair(Neutrino, partOrig);
170  }
171 
172  if (thePart&& info && info->Mother() && HepMC::is_same_generator_particle(thePart,info->Mother()))
173  return std::make_pair(NonPrimary, partOrig);
174 
175  if (isPartHadr) return std::make_pair(Hadron, partOrig);
176 
177  if (partOriVert && motherPDG == 0 && partOriVert->nOutgoingParticles() == 1 &&
178  partOriVert->nIncomingParticles() == 0) {
179  if (MC::isElectron(thePart)) {
180  info->particleOutCome = defOutComeOfElectron(thePart);
181  return std::make_pair(IsoElectron, SingleElec);
182  }
183  if (MC::isMuon(thePart)) {
184  info->particleOutCome = defOutComeOfMuon(thePart);
185  return std::make_pair(IsoMuon, SingleMuon);
186  }
187  if (MC::isTau(thePart)) {
188  info->particleOutCome = defOutComeOfTau(thePart);
189  return std::make_pair(IsoTau, SingleTau);
190  }
191  if (MC::isPhoton(thePart)) {
192  info->particleOutCome = defOutComeOfPhoton(thePart);
193  return std::make_pair(IsoPhoton, SinglePhot);
194  }
195  }
196 
197  if (motherPDG == thePart->pdg_id() && theMoth && theMoth->status() == 3 && MC::isDecayed(thePart)) return std::make_pair(GenParticle, partOrig);
198 
199  if (MC::isElectron(thePart)) {
200  bool isPrompt = false;
201  partOrig = defOrigOfElectron(truthParticleContainerReadHandle.ptr(), thePart, isPrompt, info);
202  partType = defTypeOfElectron(partOrig, isPrompt);
203  } else if (MC::isMuon(thePart)) {
204  bool isPrompt = false;
205  partOrig = defOrigOfMuon(truthParticleContainerReadHandle.ptr(), thePart, isPrompt, info);
206  partType = defTypeOfMuon(partOrig, isPrompt);
207  } else if (MC::isTau(thePart)) {
208  partOrig = defOrigOfTau(truthParticleContainerReadHandle.ptr(), thePart, motherPDG, info);
209  partType = defTypeOfTau(partOrig);
210  } else if (MC::isPhoton(thePart)) {
211  bool isPrompt = false;
212  partOrig = defOrigOfPhoton(truthParticleContainerReadHandle.ptr(), thePart, isPrompt, info);
213  partType = defTypeOfPhoton(partOrig);
214  } else if (MC::isNeutrino(thePart)) {
215  bool isPrompt = false;
216  partOrig = defOrigOfNeutrino(truthParticleContainerReadHandle.ptr(), thePart, isPrompt, info);
217  partType = Neutrino;
218  }
219 
220  ATH_MSG_DEBUG("particleTruthClassifier succeeded ");
221  return std::make_pair(partType, partOrig);
222 }
223 
224 
225 
227  const xAOD::TruthParticle* thePart,
228  bool& isPrompt,
229  MCTruthPartClassifier::Info* infoin) const
230 {
232  ATH_MSG_DEBUG("Executing DefOrigOfElectron ");
233 
234  const xAOD::TruthParticle* thePriPart = MC::findMatching(mcTruthTES, thePart);
235  if (!thePriPart) return NonDefined;
236  if (!MC::isElectron(thePriPart)) return NonDefined;
237 
238  const xAOD::TruthVertex* partOriVert = thePriPart->hasProdVtx() ? thePriPart->prodVtx() : nullptr;
239 
240  //-- to define electron outcome status
241  if (info) info->particleOutCome = defOutComeOfElectron(thePriPart);
243  if (!info) { info = &tmpinfo; }
244 
245  if (!partOriVert) return NonDefined;
246 
247  int numOfParents = -1;
248  numOfParents = partOriVert->nIncomingParticles();
249  if (numOfParents > 1) ATH_MSG_DEBUG("DefOrigOfElectron:: electron has more than one mother ");
250 
251  const xAOD::TruthParticle* mother = MC::findMother(thePriPart);
252  if (info) info->setMotherProperties(mother);
253  if (!mother) {
254  return NonDefined;
255  }
256  int motherPDG = mother->pdgId();
257  if (info) info->setMotherProperties(mother);
258  const xAOD::TruthVertex* mothOriVert = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
259 
260  bool samePart = false;
261  for (const auto& theDaug: partOriVert->particles_out()) {
262  if (!theDaug) continue;
263  if (motherPDG == theDaug->pdgId() && info && info->Mother() && HepMC::is_same_generator_particle(theDaug, info->Mother())) samePart = true;
264  }
265 
266  // to resolve Sherpa loop
267  if (mothOriVert && HepMC::is_same_vertex(mothOriVert,partOriVert)) samePart = true;
268  //
269 
270  if ((MC::isMuon(motherPDG) || MC::isTau(motherPDG) || MC::isW(motherPDG)) && mothOriVert != nullptr && !samePart) {
271  int pPDG(0);
272  const xAOD::TruthParticle* MotherParent(nullptr);
273  do {
274  pPDG = 0;
275  MotherParent = MC::findMother(mother);
276  // to prevent Sherpa loop
277  const xAOD::TruthVertex* mother_prdVtx(nullptr);
278  const xAOD::TruthVertex* mother_endVtx(nullptr);
279  if (mother) {
280  mother_prdVtx = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
281  mother_endVtx = mother->decayVtx();
282  }
283  const xAOD::TruthVertex* parent_prdVtx(nullptr);
284  const xAOD::TruthVertex* parent_endVtx(nullptr);
285  if (MotherParent) {
286  parent_prdVtx = MotherParent->hasProdVtx() ? MotherParent->prodVtx() : nullptr;
287  parent_endVtx = MotherParent->decayVtx();
288  }
289  if (mother_endVtx == parent_prdVtx && mother_prdVtx == parent_endVtx) {
290  MotherParent = mother;
291  break;
292  }
293  //
294  if (MotherParent) pPDG = MotherParent->pdgId();
295  // to prevent Sherpa loop
296  if (mother == MotherParent) break;
297  if (MC::isMuon(pPDG) || MC::isTau(pPDG) || MC::isW(pPDG)) mother = MotherParent;
298 
299  } while ((MC::isMuon(pPDG) || MC::isTau(pPDG) || MC::isW(pPDG)));
300 
301  if (MC::isMuon(pPDG) || MC::isTau(pPDG) || MC::isW(pPDG) || MC::isZ(pPDG) || MC::isHiggs(pPDG) ||
302  MC::isMSSMHiggs(pPDG) || MC::isHeavyBoson(pPDG) || MC::isTop(pPDG) || // MSSM Higgs bosons, Heavy bosons( Z', Z'', W'+)
303  abs(pPDG) == 9900024 || abs(pPDG) == 9900012 || abs(pPDG) == 9900014 || abs(pPDG) == 9900016 || // Left-right symmetric model WBoson || Right-handed NU_E || Right-handed NU_MU || Right-handed NU_TAU
304  (abs(pPDG) < 2000040 && abs(pPDG) > 1000001)) // FIXME This is nearly MC::isSUSY(pPDG), but slightly stricter.
305  mother = MotherParent;
306  }
307 
308  motherPDG = mother->pdgId();
309  partOriVert = mother->decayVtx();
310  mothOriVert = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
311  numOfParents = partOriVert->nIncomingParticles();
312  int numOfDaug = partOriVert->nOutgoingParticles();
313 
314  if (info) info->setMotherProperties(mother);
315 
316  int NumOfPhot(0);
317  int NumOfEl(0);
318  int NumOfPos(0);
319  int NumOfNucFr(0);
320  int NumOfquark(0);
321  int NumOfgluon(0);
322  int NumOfElNeut(0);
323  int NumOfLQ(0);
324  int NumOfMuPl(0);
325  int NumOfMuMin(0);
326  int NumOfMuNeut(0);
327  int NumOfTau(0);
328  int NumOfTauNeut(0);
329  samePart = false;
330 
331  for (const auto& theDaug: partOriVert->particles_out()) {
332  if (!theDaug) continue;
333  int DaugType = theDaug->pdgId();
334  if (MC::isSMQuark(DaugType)) NumOfquark++;
335  else if (MC::isGluon(DaugType)) NumOfgluon++;
336  else if (abs(DaugType) == MC::NU_E) NumOfElNeut++;
337  else if (abs(DaugType) == MC::NU_MU) NumOfMuNeut++;
338  else if (MC::isPhoton(DaugType)) NumOfPhot++;
339  else if (DaugType == MC::ELECTRON) NumOfEl++;
340  else if (DaugType == MC::POSITRON) NumOfPos++;
341  else if (DaugType == MC::MUON) NumOfMuMin++;
342  else if (DaugType == -MC::MUON) NumOfMuPl++;
343  else if (MC::isTau(DaugType)) NumOfTau++;
344  else if (abs(DaugType) == MC::NU_TAU) NumOfTauNeut++;
345  else if (MC::isLeptoQuark(DaugType)) NumOfLQ++;
346  if (abs(DaugType) == abs(motherPDG) && theDaug && info && info->Mother() && HepMC::is_same_generator_particle(theDaug, info->Mother() )) samePart = true;
347  if (numOfParents == 1 &&
348  (MC::isPhoton(motherPDG) || MC::isElectron(motherPDG) || MC::isMuon(motherPDG) || abs(motherPDG) == MC::PIPLUS) &&
349  (DaugType > 1000000000 || DaugType == 0 || DaugType == MC::PROTON || DaugType == MC::NEUTRON ||
350  abs(DaugType) == MC::PIPLUS || abs(DaugType) == MC::PI0))
351  NumOfNucFr++;
352  }
353 
354  if (MC::isPhoton(motherPDG) && mothOriVert != nullptr) {
355  for (const auto& theMother: mothOriVert->particles_in()) {
356  if (!theMother) continue;
357  if (!info) continue;
358  info->photonMother = theMother;
359  }
360  }
361 
362  if ((MC::isPhoton(motherPDG) && numOfDaug == 2 && NumOfEl == 1 && NumOfPos == 1) || (MC::isPhoton(motherPDG) && numOfDaug == 1 && (NumOfEl == 1 || NumOfPos == 1))) return PhotonConv;
363 
364  // e,gamma,pi+Nuclear->NuclearFragments+nuclons+e
365  if ((numOfParents == 1 && (MC::isPhoton(motherPDG) || MC::isElectron(motherPDG) || MC::isTau(motherPDG))) && numOfDaug > 1 && NumOfNucFr != 0) return ElMagProc;
366 
367  if (numOfParents == 1 && abs(motherPDG) == MC::PIPLUS && numOfDaug > 2 && NumOfNucFr != 0) return ElMagProc;
368 
369  // nuclear photo fission
370  if (MC::isPhoton(motherPDG) && numOfDaug > 4 && NumOfNucFr != 0) return ElMagProc;
371 
372  // unknown process el(pos)->el+pos??
373  if (MC::isElectron(motherPDG) && numOfDaug == 2 && NumOfEl == 1 && NumOfPos == 1) return ElMagProc;
374 
375  // unknown process el->el+el??
376  if (motherPDG == MC::ELECTRON && numOfDaug == 2 && NumOfEl == 2 && NumOfPos == 0) return ElMagProc;
377 
378  // unknown process pos->pos+pos??
379  if (motherPDG == MC::POSITRON && numOfDaug == 2 && NumOfEl == 0 && NumOfPos == 2) return ElMagProc;
380 
381  // unknown process pos/el->pos/el??
382  if (MC::isElectron(motherPDG) && !MC::isDecayed(mother) && motherPDG == thePriPart->pdgId() && numOfDaug == 1 && !samePart) return ElMagProc;
383 
384  // pi->pi+e+/e-; mu->mu+e+/e- ;
385  // gamma+ atom->gamma(the same) + e (compton scattering)
386  if (numOfDaug == 2 && (NumOfEl == 1 || NumOfPos == 1) && !MC::isElectron(motherPDG) && samePart) return ElMagProc;
387 
388  if ((motherPDG == MC::PI0 && numOfDaug == 3 && NumOfPhot == 1 && NumOfEl == 1 && NumOfPos == 1) ||
389  (motherPDG == MC::PI0 && numOfDaug == 4 && NumOfPhot == 0 && NumOfEl == 2 && NumOfPos == 2))
390  return DalitzDec;
391 
392  // Quark weak decay
393  if (MC::isSMQuark(motherPDG) && numOfParents == 1 && numOfDaug == 3 && NumOfquark == 1 && NumOfElNeut == 1) return QuarkWeakDec;
394 
395  if (MC::isMuon(motherPDG) && NumOfNucFr != 0) return ElMagProc;
396 
397  if (MC::isTop(motherPDG)) return top;
398 
399  if (MC::isW(motherPDG) && mothOriVert != nullptr && mothOriVert->nIncomingParticles() != 0) {
400 
401  const xAOD::TruthVertex* prodVert = mothOriVert;
402  const xAOD::TruthParticle* ptrPart;
403  do {
404  ptrPart = prodVert->incomingParticle(0);
405  prodVert = ptrPart->hasProdVtx() ? ptrPart->prodVtx() : nullptr;
406  } while (MC::isW(ptrPart) && prodVert != nullptr);
407 
408  if (prodVert && prodVert->nIncomingParticles() == 1) {
409  if (abs(ptrPart->pdgId()) == 9900012) return NuREle; // Right-handed NU_E
410  if (abs(ptrPart->pdgId()) == 9900014) return NuRMu; // Right-handed NU_MU
411  if (abs(ptrPart->pdgId()) == 9900016) return NuRTau; // Right-handed NU_TAU
412  }
413  return WBoson;
414  }
415  if (MC::isW(motherPDG)) return WBoson;
416  if (MC::isZ(motherPDG)) return ZBoson;
417 
418  // MadGraphPythia ZWW*->lllnulnu
419  if (numOfParents == 1 && numOfDaug > 4 && (MC::isSMQuark(motherPDG) || MC::isGluon(motherPDG))) {
420 
421  const xAOD::TruthParticle* thePartToCheck = thePriPart;
422  const xAOD::TruthParticle* theMother = thePriPart->hasProdVtx() ? thePriPart->prodVtx()->incomingParticle(0) : nullptr;
423  if (theMother != nullptr && MC::isElectron(theMother) && MC::isDecayed(theMother)) thePartToCheck = theMother;
424 
425  bool isZboson = false;
426  bool isWboson = false;
427  bool skipnext = false;
428 
429  for (unsigned int ipOut = 0; ipOut + 1 < partOriVert->nOutgoingParticles(); ++ipOut) {
430  const xAOD::TruthParticle* theDaug = partOriVert->outgoingParticle(ipOut);
431  if (!theDaug) continue;
432  const xAOD::TruthParticle* theNextDaug = nullptr;
433  for (unsigned int ipOut1 = ipOut + 1; ipOut1 < partOriVert->nOutgoingParticles(); ipOut1++) {
434  theNextDaug = partOriVert->outgoingParticle(ipOut1);
435  if (theNextDaug != nullptr) break;
436  }
437  if (!theNextDaug) continue;
438  if (skipnext) {
439  skipnext = false;
440  continue;
441  }
442 
443  if (MC::isElectron(theDaug) && MC::isElectron(theNextDaug)) {
444  // Zboson
445  if (thePartToCheck == theDaug || thePartToCheck == theNextDaug) {
446  isZboson = true;
447  break;
448  }
449  skipnext = true;
450  } else if (MC::isElectron(theDaug) && abs(theNextDaug->pdgId()) == MC::NU_E) {
451  // WBoson
452  if (thePartToCheck == theDaug || thePartToCheck == theNextDaug) {
453  isWboson = true;
454  break;
455  }
456  skipnext = true;
457  }
458  }
459  if (isWboson) return WBoson;
460  if (isZboson) return ZBoson;
461  }
462  if (numOfParents == 2) {
463  const int pdg1 = partOriVert->incomingParticle(0)->pdgId();
464  const int pdg2 = partOriVert->incomingParticle(1)->pdgId();
465  //--Sherpa Z->ee
466  if ((numOfDaug - NumOfquark - NumOfgluon) == 2 && NumOfEl == 1 && NumOfPos == 1) return ZBoson;
467 
468  //--Sherpa W->enu ??
469  if ((numOfDaug - NumOfquark - NumOfgluon) == 2 && (NumOfEl == 1 || NumOfPos == 1) && NumOfElNeut == 1) return WBoson;
470 
471  //--Sherpa ZZ,ZW
472  if ((numOfDaug - NumOfquark - NumOfgluon) == 4 && (NumOfEl + NumOfPos + NumOfMuPl + NumOfMuMin + NumOfTau + NumOfElNeut + NumOfMuNeut + NumOfTauNeut == 4) &&
473  (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return DiBoson;
474 
475  //--Sherpa VVV -- Note, have to allow for prompt photon radiation or these get lost
476  if ((numOfDaug - NumOfquark - NumOfgluon - NumOfPhot) == 6 && (NumOfEl + NumOfPos + NumOfMuPl + NumOfMuMin + NumOfTau + NumOfElNeut + NumOfMuNeut + NumOfTauNeut == 6) &&
477  (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return MultiBoson;
478  }
479 
480  // New Sherpa Z->ee
481  if (partOriVert == mothOriVert && partOriVert != nullptr) {
482  int NumOfEleLoop = 0;
483  int NumOfLepLoop = 0;
484  int NumOfEleNeuLoop = 0;
485  for (const auto *const pout: partOriVert->particles_out()) {
486  if (!pout) continue;
487  for (const auto *const pin: partOriVert->particles_in()) {
488  if (!pin) continue;
489  if (!HepMC::is_same_particle(pout,pin)) continue;
490  if (MC::isElectron(pout)) NumOfEleLoop++;
491  if (std::abs(pin->pdgId()) == MC::NU_E) NumOfEleNeuLoop++;
492  if (MC::isSMLepton(pout)) NumOfLepLoop++;
493  }
494  }
495  if (NumOfEleLoop == 2 && NumOfEleNeuLoop == 0) return ZBoson;
496  if (NumOfEleLoop == 1 && NumOfEleNeuLoop == 1) return WBoson;
497  if ((NumOfEleLoop == 4 && NumOfEleNeuLoop == 0) || (NumOfEleLoop == 3 && NumOfEleNeuLoop == 1) ||
498  (NumOfEleLoop == 2 && NumOfEleNeuLoop == 2))
499  return DiBoson;
500  if (NumOfLepLoop == 4) return DiBoson;
501  }
502 
503  //-- McAtNLo
504 
505  if (MC::isHiggs(motherPDG)) return Higgs;
506 
507  if (MC::isMSSMHiggs(motherPDG)) return HiggsMSSM; // MSSM Higgs bosons
508 
509  if (MC::isHeavyBoson(motherPDG)) return HeavyBoson; // Heavy bosons( Z', Z'', W'+)
510 
511  if (MC::isMuon(motherPDG)) return Mu;
512  if (MC::isTau(motherPDG)) {
513  ParticleOrigin tauOrig = defOrigOfTau(mcTruthTES, mother, motherPDG, info);
514  ParticleType tautype = defTypeOfTau(tauOrig);
515  return (tautype == IsoTau)?tauOrig:TauLep;
516  }
517 
518  if (abs(motherPDG) == 9900024) return WBosonLRSM; // Left-right symmetric model WBoson
519  if (abs(motherPDG) == 9900012) return NuREle; // Right-handed NU_E
520  if (abs(motherPDG) == 9900014) return NuRMu; // Right-handed NU_MU
521  if (abs(motherPDG) == 9900016) return NuRTau; // Right-handed NU_TAU
522  if (MC::isLeptoQuark(motherPDG) || NumOfLQ != 0) return LQ;
523  if (MC::isSUSY(motherPDG)) return SUSY;
524  if (MC::isBSM(motherPDG)) return OtherBSM;
525  ParticleType pType = defTypeOfHadron(motherPDG);
526  if ((pType == BBbarMesonPart || pType == CCbarMesonPart) && mothOriVert != nullptr && MC::isHardScatteringVertex(mothOriVert)) isPrompt = true;
527  return convHadronTypeToOrig(pType, motherPDG);
528 }
529 
530 
532  const xAOD::TruthParticle* thePart,
533  bool& isPrompt,
534  MCTruthPartClassifier::Info* infoin) const
535 {
537  ATH_MSG_DEBUG("Executing DefOrigOfMuon ");
538 
539  const xAOD::TruthParticle* thePriPart = MC::findMatching(mcTruthTES, thePart);
540  if (!thePriPart) return NonDefined;
541  if (!MC::isMuon(thePriPart)) return NonDefined;
542 
543  const xAOD::TruthVertex* partOriVert = thePriPart->hasProdVtx() ? thePriPart->prodVtx() : nullptr;
544 
545  //-- to define muon outcome status
546  if (info) info->particleOutCome = defOutComeOfMuon(thePriPart);
547 
549  if (!info) { info = &tmpinfo; }
550  if (!partOriVert) return NonDefined;
551 
552  int numOfParents = partOriVert->nIncomingParticles();
553  if (numOfParents > 1) ATH_MSG_DEBUG("DefOrigOfMuon:: muon has more than one mother ");
554 
555  const xAOD::TruthParticle* mother = MC::findMother(thePriPart);
556  if (info) info->setMotherProperties(mother);
557  if (!mother) {
558  return NonDefined;
559  }
560 
561  const xAOD::TruthVertex* mothOriVert = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
562  int motherPDG = mother->pdgId();
563 
564  if ((MC::isTau(motherPDG)|| MC::isW(motherPDG)) && mothOriVert != nullptr) {
565  int pPDG(0);
566  const xAOD::TruthParticle* MotherParent(nullptr);
567  do {
568  //
569  pPDG = 0;
570  //
571  const xAOD::TruthVertex* mother_prdVtx(nullptr);
572  const xAOD::TruthVertex* mother_endVtx(nullptr);
573  MotherParent = MC::findMother(mother);
574  // to prevent Sherpa loop
575  mother_prdVtx = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
576  mother_endVtx = mother->decayVtx();
577  //
578  const xAOD::TruthVertex* parent_prdVtx(nullptr);
579  const xAOD::TruthVertex* parent_endVtx(nullptr);
580  if (MotherParent) {
581  parent_prdVtx = MotherParent->hasProdVtx() ? MotherParent->prodVtx() : nullptr;
582  parent_endVtx = MotherParent->decayVtx();
583  }
584  //
585  if (mother_endVtx == parent_prdVtx && mother_prdVtx == parent_endVtx) {
586  MotherParent = mother;
587  break;
588  }
589  //
590  // to prevent Sherpa loop
591  if (mother == MotherParent) {
592  break;
593  }
594 
595  if (MotherParent) {
596  pPDG = MotherParent->pdgId();
597  if (MC::isMuon(pPDG) || MC::isTau(pPDG) || MC::isW(pPDG)) {
598  mother = MotherParent;
599  if (info) info->setMotherProperties(mother);
600  }
601  }
602  } while ((MC::isMuon(pPDG) || MC::isTau(pPDG) || MC::isW(pPDG)));
603 
604  if (MC::isTau(pPDG) || MC::isW(pPDG) || MC::isZ(pPDG) || MC::isHiggs(pPDG) ||
605  MC::isMSSMHiggs(pPDG) || MC::isHeavyBoson(pPDG) || MC::isTop(pPDG) || // MSSM Higgs bosons, Heavy bosons( Z', Z'', W'+)
606  abs(pPDG) == 9900024 || abs(pPDG) == 9900012 || abs(pPDG) == 9900014 || abs(pPDG) == 9900016 || // Left-right symmetric model WBoson || Right-handed NU_E || Right-handed NU_MU || Right-handed NU_TAU
607  (abs(pPDG) < 2000040 && abs(pPDG) > 1000001)) { // FIXME This is nearly MC::isSUSY(pPDG), but slightly stricter.
608  if (info) info->setMotherProperties(mother);
609  }
610  }
611 
612  motherPDG = mother->pdgId();
613  mothOriVert = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
614  partOriVert = mother->decayVtx();
615  numOfParents = partOriVert->nIncomingParticles();
616  int numOfDaug = partOriVert->nOutgoingParticles();
617 
618  if (info) info->setMotherProperties(mother);
619  auto DP = DecayProducts(partOriVert);
620  const int NumOfPhot = DP.pd(MC::PHOTON);
621  const int NumOfEl = DP.pd(MC::ELECTRON);
622  const int NumOfPos = DP.pd(MC::POSITRON);
623  const int NumOfElNeut = DP.apd(MC::NU_E);
624  const int NumOfMuNeut = DP.apd(MC::NU_MU);
625  const int NumOfLQ = DP.apd(MC::LEPTOQUARK);
626  const int NumOfquark = DP.apd({MC::DQUARK,MC::UQUARK,MC::SQUARK,MC::CQUARK,MC::BQUARK,MC::TQUARK});
627  const int NumOfgluon = DP.apd(MC::GLUON);
628  const int NumOfMuPl = DP.pd(-MC::MUON);
629  const int NumOfMuMin = DP.pd(MC::MUON);
630  const int NumOfTau = DP.apd(MC::TAU);
631  const int NumOfTauNeut = DP.apd(MC::NU_TAU);
632  if (std::abs(motherPDG) == MC::PIPLUS && numOfDaug == 2 && NumOfMuNeut == 1) return PionDecay;
633  if (std::abs(motherPDG) == MC::KPLUS && numOfDaug == 2 && NumOfMuNeut == 1) return KaonDecay;
634  if (MC::isTau(motherPDG)) {
635  ParticleOrigin tauOrig = defOrigOfTau(mcTruthTES, mother, motherPDG, info);
636  ParticleType tautype = defTypeOfTau(tauOrig);
637  return (tautype == IsoTau)?tauOrig:TauLep;
638  }
639 
640  if (MC::isTop(motherPDG)) return top;
641  // Quark weak decay
642  if (MC::isSMQuark(motherPDG) && numOfParents == 1 && numOfDaug == 3 && NumOfquark == 1 && NumOfMuNeut == 1) return QuarkWeakDec;
643 
644  if (MC::isW(motherPDG) && mothOriVert != nullptr && mothOriVert->nIncomingParticles() != 0) {
645  const xAOD::TruthVertex* prodVert = mothOriVert;
646  const xAOD::TruthParticle* itrP;
647  do {
648  itrP = prodVert->incomingParticle(0);
649  prodVert = itrP->hasProdVtx() ? itrP->prodVtx() : nullptr;
650  } while (MC::isW(itrP) && prodVert != nullptr);
651 
652  if (prodVert && prodVert->nIncomingParticles() == 1) {
653  if (abs(itrP->pdgId()) == 9900012) return NuREle; // Right-handed NU_E
654  if (abs(itrP->pdgId()) == 9900014) return NuRMu; // Right-handed NU_MU
655  if (abs(itrP->pdgId()) == 9900016) return NuRTau; // Right-handed NU_TAU
656  }
657  return WBoson;
658  }
659  if (MC::isW(motherPDG)) return WBoson;
660  if (MC::isZ(motherPDG)) return ZBoson;
661  if (MC::isPhoton(motherPDG) && numOfDaug == 2 && NumOfMuMin == 1 && NumOfMuPl == 1) return PhotonConv;
662  //-- Exotics
663 
664  // MadGraphPythia ZWW*->lllnulnu
665 
666  if (numOfParents == 1 && numOfDaug > 4 && (MC::isSMQuark(motherPDG) || MC::isGluon(motherPDG))) {
667  bool isZboson = false;
668  bool isWboson = false;
669  bool skipnext = false;
670  for (unsigned int ipOut = 0; ipOut + 1 < partOriVert->nOutgoingParticles(); ipOut++) {
671  if (skipnext) {
672  skipnext = false;
673  continue;
674  }
675  const xAOD::TruthParticle* theDaug = partOriVert->outgoingParticle(ipOut);
676  if (!theDaug) continue;
677  const xAOD::TruthParticle* theNextDaug = nullptr;
678  for (unsigned int ipOut1 = ipOut + 1; ipOut1 < partOriVert->nOutgoingParticles(); ipOut1++) {
679  theNextDaug = partOriVert->outgoingParticle(ipOut1);
680  if (theNextDaug != nullptr) break;
681  }
682  if (!theNextDaug) continue;
683  if (MC::isMuon(theDaug) && MC::isMuon(theNextDaug)) {
684  // Zboson
685  if (thePriPart == theDaug || thePriPart == theNextDaug) {
686  isZboson = true;
687  break;
688  }
689  skipnext = true;
690  } else if (MC::isMuon(theDaug) && abs(theNextDaug->pdgId()) == MC::NU_MU) {
691  // WBoson
692  if (thePriPart == theDaug || thePriPart == theNextDaug) {
693  isWboson = true;
694  break;
695  }
696  skipnext = true;
697  }
698  }
699  if (isWboson) return WBoson;
700  if (isZboson) return ZBoson;
701  }
702  if (numOfParents == 2 ) {
703  const int pdg1 = partOriVert->incomingParticle(0)->pdgId();
704  const int pdg2 = partOriVert->incomingParticle(1)->pdgId();
705  //--Sherpa Z->mumu
706  if ((numOfDaug - NumOfquark - NumOfgluon) == 2 && NumOfMuPl == 1 && NumOfMuMin == 1) return ZBoson;
707 
708  //--Sherpa W->munu ??
709  // if(numOfParents==2&&(numOfDaug-NumOfquark-NumOfgluon)==2&&(NumOfEl==1||NumOfPos==1)&&NumOfElNeut==1) return WBoson;
710  if ((numOfDaug - NumOfquark - NumOfgluon) == 2 && (NumOfMuPl == 1 || NumOfMuMin == 1) && NumOfMuNeut == 1) return WBoson;
711 
712  //--Sherpa ZZ,ZW
713  if ((numOfDaug - NumOfquark - NumOfgluon) == 4 &&
714  (NumOfEl + NumOfPos + NumOfMuPl + NumOfMuMin + NumOfTau + NumOfElNeut + NumOfMuNeut + NumOfTauNeut == 4) &&
715  (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return DiBoson;
716 
717  //--Sherpa VVV -- Note, have to allow for prompt photon radiation or these get lost
718  if ((numOfDaug - NumOfquark - NumOfgluon - NumOfPhot) == 6 &&
719  (NumOfEl + NumOfPos + NumOfMuPl + NumOfMuMin + NumOfTau + NumOfElNeut + NumOfMuNeut + NumOfTauNeut == 6) &&
720  (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return MultiBoson;
721 }
722 
723  //--New Sherpa Z->mumu
724  if (partOriVert == mothOriVert) {
725  int NumOfMuLoop = 0;
726  int NumOfMuNeuLoop = 0;
727  int NumOfLepLoop = 0;
728  for (const auto & pout: partOriVert->particles_out()) {
729  for (const auto & pin: partOriVert->particles_in()) {
730  if (!pout) continue;
731  if (!pin) continue;
732  if (HepMC::is_same_particle(pout,pin)) {
733  if (MC::isMuon(pout)) NumOfMuLoop++;
734  if (std::abs(pout->pdg_id()) == MC::NU_MU) NumOfMuNeuLoop++;
735  if (MC::isSMLepton(pout)) NumOfLepLoop++;
736  }
737  }
738  }
739  if (NumOfMuLoop == 2 && NumOfMuNeuLoop == 0) return ZBoson;
740  if (NumOfMuLoop == 1 && NumOfMuNeuLoop == 1) return WBoson;
741  if ((NumOfMuLoop == 4 && NumOfMuNeuLoop == 0) || (NumOfMuLoop == 3 && NumOfMuNeuLoop == 1) || (NumOfMuLoop == 2 && NumOfMuNeuLoop == 2)) return DiBoson;
742  if (NumOfLepLoop == 4) return DiBoson;
743  }
744 
745  //-- McAtNLo
746 
747  if (MC::isHiggs(motherPDG)) return Higgs;
748 
749  if (MC::isMSSMHiggs(motherPDG)) return HiggsMSSM; // MSSM Higgs bosons
750 
751  if (MC::isHeavyBoson(motherPDG)) return HeavyBoson; // Heavy bosons( Z', Z'', W'+)
752 
753  if (abs(motherPDG) == 9900024) return WBosonLRSM; // Left-right symmetric model WBoson
754  if (abs(motherPDG) == 9900012) return NuREle; // Right-handed NU_E
755  if (abs(motherPDG) == 9900014) return NuRMu; // Right-handed NU_MU
756  if (abs(motherPDG) == 9900016) return NuRTau; // Right-handed NU_TAU
757  if (MC::isLeptoQuark(motherPDG) || NumOfLQ != 0) return LQ;
758  if (MC::isSUSY(motherPDG)) return SUSY;
759  if (MC::isBSM(motherPDG)) return OtherBSM;
760 
761  ParticleType pType = defTypeOfHadron(motherPDG);
762  if ((pType == BBbarMesonPart || pType == CCbarMesonPart) && mothOriVert != nullptr && MC::isHardScatteringVertex(mothOriVert)) isPrompt = true;
763 
764  return convHadronTypeToOrig(pType, motherPDG);
765 }
766 
768  const xAOD::TruthParticle* thePart,
769  int motherPDG,
770  MCTruthPartClassifier::Info* infoin) const
771 {
773  ATH_MSG_DEBUG("Executing DefOrigOfTau ");
774  const xAOD::TruthParticle* thePriPart = MC::findMatching(mcTruthTES, thePart);
775  if (!thePriPart) return NonDefined;
776  if (!MC::isTau(thePriPart)) return NonDefined;
777 
778  const xAOD::TruthVertex* partOriVert = thePriPart->hasProdVtx() ? thePriPart->prodVtx() : nullptr;
779 
780  //-- to define tau outcome status
781  if (MC::isPhysical(thePriPart) && info) info->particleOutCome = defOutComeOfTau(thePriPart);
782 
784  if (!info) { info = &tmpinfo; }
785  if (!partOriVert) return NonDefined;
786 
787  int numOfParents = partOriVert->nIncomingParticles();
788  if (numOfParents > 1) ATH_MSG_DEBUG("DefOrigOfTau:: tau has more than one mother ");
789 
790  const xAOD::TruthParticle* mother = MC::findMother(thePriPart);
791  if (info) info->setMotherProperties(mother);
792  if (!mother) {
793  return NonDefined;
794  }
795 
796  const xAOD::TruthVertex* mothOriVert = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
797 
798  const xAOD::TruthParticle* MotherParent(nullptr);
799 
800  if (MC::isW(motherPDG) && mothOriVert != nullptr) {
801  MotherParent = MC::findMother(mother);
802  int pPDG(0);
803 
804  if (MotherParent) {//MotherParent checked here...
805  pPDG = MotherParent->pdgId();
806  if (MC::isTop(pPDG)) {
807  mother = MotherParent; //...so mother cannot be nullptr
808  if (info) info->setMotherProperties(mother);
809  }
810  }
811  }
812 
813  motherPDG = mother->pdgId();
814  if (info) info->setMotherProperties(mother);
815  mothOriVert = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
816  partOriVert = mother->decayVtx();
817  if (!partOriVert) return NonDefined;
818 
819  numOfParents = partOriVert->nIncomingParticles();
820  auto DP = DecayProducts(partOriVert);
821  const int numOfDaug = DP.size();
822  const int NumOfPhot = DP.pd(MC::PHOTON);
823  const int NumOfEl = DP.pd(MC::ELECTRON);
824  const int NumOfPos = DP.pd(MC::POSITRON);
825  const int NumOfElNeut = DP.apd(MC::NU_E);
826  const int NumOfMuNeut = DP.apd(MC::NU_MU);
827  /* const int NumOfLQ = DP.apd(MC::LEPTOQUARK); */
828  const int NumOfquark = DP.apd({MC::DQUARK,MC::UQUARK,MC::SQUARK,MC::CQUARK,MC::BQUARK,MC::TQUARK});
829  const int NumOfgluon = DP.apd(MC::GLUON);
830  const int NumOfMuPl = DP.pd(-MC::MUON);
831  const int NumOfMuMin = DP.pd(MC::MUON);
832  const int NumOfTau = DP.apd(MC::TAU);
833  const int NumOfTauNeut = DP.apd(MC::NU_TAU);
834 
835  if (MC::isTop(motherPDG)) return top;
836  if (MC::isW(motherPDG) && mothOriVert != nullptr && mothOriVert->nIncomingParticles() != 0) {
837  const xAOD::TruthVertex* prodVert = mothOriVert;
838  const xAOD::TruthParticle* itrP;
839  do {
840  itrP = prodVert->incomingParticle(0);
841  prodVert = itrP->hasProdVtx() ? itrP->prodVtx() : nullptr;
842  } while (MC::isW(itrP) && prodVert != nullptr);
843 
844  if (prodVert && prodVert->nIncomingParticles() == 1 ) {
845  if (abs(itrP->pdgId()) == 9900012) return NuREle; // Right-handed NU_E
846  if (abs(itrP->pdgId()) == 9900014) return NuRMu; // Right-handed NU_MU
847  if (abs(itrP->pdgId()) == 9900016) return NuRTau; // Right-handed NU_TAU
848  }
849  return WBoson;
850  }
851  if (MC::isW(motherPDG)) { return WBoson;}
852  if (MC::isZ(motherPDG)) { return ZBoson;}
853  if (numOfParents == 1 && numOfDaug > 4 && (MC::isSMQuark(motherPDG) || MC::isGluon(motherPDG))) {
854  bool isZboson = false;
855  bool isWboson = false;
856  bool skipnext = false;
857  for (unsigned int ipOut = 0; ipOut + 1 < partOriVert->nOutgoingParticles(); ipOut++) {
858  if (skipnext) {
859  skipnext = false;
860  continue;
861  }
862  const xAOD::TruthParticle* theDaug = partOriVert->outgoingParticle(ipOut);
863  if (!theDaug) continue;
864  const xAOD::TruthParticle* theNextDaug = nullptr;
865  for (unsigned int ipOut1 = ipOut + 1; ipOut1 < partOriVert->nOutgoingParticles(); ipOut1++) {
866  theNextDaug = partOriVert->outgoingParticle(ipOut1);
867  if (theNextDaug != nullptr) break;
868  }
869  if (!theNextDaug) {
870  continue;
871  }
872  if (MC::isTau(theDaug) && MC::isTau(theNextDaug)) {
873  // Zboson
874  if (thePriPart == theDaug || thePriPart == theNextDaug) {
875  isZboson = true;
876  break;
877  }
878  skipnext = true;
879  } else if (MC::isTau(theDaug) && abs(theNextDaug->pdgId()) == MC::NU_TAU) {
880  // WBoson
881  if (thePriPart == theDaug || thePriPart == theNextDaug) {
882  isWboson = true;
883  break;
884  }
885  skipnext = true;
886  }
887  }
888  if (isWboson) return WBoson;
889  if (isZboson) return ZBoson;
890  }
891  if (numOfParents == 2 ) {
892  const int pdg1 = partOriVert->incomingParticle(0)->pdgId();
893  const int pdg2 = partOriVert->incomingParticle(1)->pdgId();
894  //--Sherpa Z->tautau
895  if ((numOfDaug - NumOfquark - NumOfgluon) == 2 && NumOfTau == 2 && (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return ZBoson;
896 
897  //--Sherpa W->taunu new
898  if ((numOfDaug - NumOfquark - NumOfgluon) == 2 && NumOfTau == 1 && NumOfTauNeut == 1) return WBoson;
899 
900  //--Sherpa ZZ,ZW
901  if ((numOfDaug - NumOfquark - NumOfgluon) == 4 &&
902  (NumOfEl + NumOfPos + NumOfMuPl + NumOfMuMin + NumOfTau + NumOfElNeut + NumOfMuNeut + NumOfTauNeut == 4) &&
903  (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return DiBoson;
904 
905  //--Sherpa VVV -- Note, have to allow for prompt photon radiation or these get lost
906  if ((numOfDaug - NumOfquark - NumOfgluon - NumOfPhot) == 6 &&
907  (NumOfEl + NumOfPos + NumOfMuPl + NumOfMuMin + NumOfTau + NumOfElNeut + NumOfMuNeut + NumOfTauNeut == 6) &&
908  (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return MultiBoson;
909  }
910  // New Sherpa Z->tautau
911  if (partOriVert == mothOriVert) {
912  int NumOfTauLoop = 0;
913  int NumOfTauNeuLoop = 0;
914  int NumOfLepLoop = 0;
915  for ( const auto *const pout: partOriVert->particles_out()) {
916  if (!pout) continue;
917  for (const auto *const pin: partOriVert->particles_in()) {
918  if (!pin) continue;
919  if (!HepMC::is_same_particle(pout,pin)) continue;
920  if (MC::isTau(pout)) NumOfTauLoop++;
921  if (std::abs(pout->pdgId()) == MC::NU_TAU) NumOfTauNeuLoop++;
922  if (MC::isSMLepton(pout)) NumOfLepLoop++;
923  }
924  }
925  if (NumOfTauLoop == 2 && NumOfTauNeuLoop == 0) return ZBoson;
926  if (NumOfTauLoop == 1 && NumOfTauNeuLoop == 1) return WBoson;
927  if ((NumOfTauLoop == 4 && NumOfTauNeuLoop == 0) || (NumOfTauLoop == 3 && NumOfTauNeuLoop == 1) || (NumOfTauLoop == 2 && NumOfTauNeuLoop == 2)) return DiBoson;
928  if (NumOfLepLoop == 4) return DiBoson;
929  }
930 
931  if (MC::isHiggs(motherPDG)) return Higgs;
932  if (MC::isMSSMHiggs(motherPDG)) return HiggsMSSM; // MSSM Higgs bosons
933  if (MC::isHeavyBoson(motherPDG)) return HeavyBoson; // Heavy bosons( Z', Z'', W'+)
934  if (abs(motherPDG) == 9900024) return WBosonLRSM; // Left-right symmetric model WBoson
935  if (abs(motherPDG) == 9900016) return NuRTau; // Right-handed NU_TAU
936  if (MC::isSUSY(motherPDG)) return SUSY;
937  if (MC::isBSM(motherPDG)) return OtherBSM;
938  if (abs(motherPDG) == MC::JPSI) return JPsi;
939 
940  ParticleType pType = defTypeOfHadron(motherPDG);
941  return convHadronTypeToOrig(pType, motherPDG);
942 }
943 
945  const xAOD::TruthParticle* thePart,
946  bool& isPrompt,
947  MCTruthPartClassifier::Info* infoin) const
948 {
949  if (!thePart) return NonDefined;
950  if (!mcTruthTES) return NonDefined;
952  ATH_MSG_DEBUG("Executing DefOrigOfPhoton ");
953 
955  if (!info) { info = &tmpinfo; }
956  if (info) {
957  info->resetMotherProperties();
958  info->photonMother = nullptr;
959  }
960 
961  const xAOD::TruthParticle* thePriPart = MC::findMatching(mcTruthTES, thePart);
962  if (!thePriPart) return NonDefined;
963  if (!MC::isPhoton(thePriPart)) return NonDefined;
964 
965  const xAOD::TruthVertex* partOriVert = thePriPart->hasProdVtx() ? thePriPart->prodVtx() : nullptr;
966 
967  //-- to define photon outcome status
968  if (info) info->particleOutCome = defOutComeOfPhoton(thePriPart);
969  if (!partOriVert) return NonDefined;
970 
971  int numOfParents = partOriVert->nIncomingParticles();
972  if (partOriVert->nIncomingParticles() > 1) ATH_MSG_DEBUG("DefOrigOfPhoton:: photon has more than one mother ");
973 
974 
975  const xAOD::TruthParticle* mother = MC::findMother(thePriPart);
976  if (info) info->setMotherProperties(mother);
977  if (!mother) return NonDefined;
978  int motherPDG = mother->pdgId();
979  const xAOD::TruthVertex* mothOriVert = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
980  if (info) info->setMotherProperties(mother);
981  partOriVert = mother->decayVtx();
982  numOfParents = partOriVert->nIncomingParticles();
983  int numOfDaug = partOriVert->nOutgoingParticles();
984  int NumOfNucFr(0);
985  int NumOfEl(0);
986  int NumOfPos(0);
987  int NumOfMu(0);
988  int NumOfTau(0);
989  int NumOfPht(0);
990  int NumOfLQ(0);
991  int DaugType(0);
992  long NumOfLep(0);
993  long NumOfNeut(0);
994  long NumOfPartons(0);
995  const xAOD::TruthParticle* Daug = nullptr;
996  for (const auto& pout: partOriVert->particles_out()) {
997  if (!pout) continue;
998  DaugType = pout->pdg_id();
999  if (numOfParents == 1 && (MC::isPhoton(motherPDG) || MC::isElectron(motherPDG) || abs(motherPDG) == MC::PIPLUS) &&
1000  (DaugType > 1000000000 || DaugType == 0 || DaugType == MC::PROTON || DaugType == MC::NEUTRON))
1001  NumOfNucFr++;
1002  if (DaugType == MC::PHOTON) NumOfPht++;
1003  else if (DaugType == MC::ELECTRON) NumOfEl++;
1004  else if (DaugType == MC::POSITRON) NumOfPos++;
1005  else if (MC::isMuon(DaugType)) NumOfMu++;
1006  else if (MC::isTau(DaugType)) NumOfTau++;
1007  else if (MC::isLeptoQuark(DaugType)) NumOfLQ++;
1008  else if (MC::isElectron(DaugType) || MC::isMuon(DaugType) || MC::isTau(DaugType)) NumOfLep++;
1009  else if (abs(DaugType) == MC::NU_E || abs(DaugType) == MC::NU_MU || abs(DaugType) == MC::NU_TAU) NumOfNeut++;
1010  if (abs(DaugType) < MC::ELECTRON || (abs(DaugType) > MC::NU_TAU && abs(DaugType) < 43 && !MC::isPhoton(DaugType))) NumOfPartons++; // FIXME Too loose? This definition picks up 4th generation quarks and leptons as well as all gauge bosons and leptoquarks.
1011  if (DaugType == motherPDG) {
1012  Daug = pout;
1013  }
1014  }
1015 
1016  bool foundISR = false;
1017  bool foundFSR = false;
1018  if (numOfParents == 1 && numOfDaug == 2 && Daug && info && info->Mother() && HepMC::is_same_generator_particle(Daug, info->Mother())) return BremPhot;
1019  if (numOfParents == 1 && numOfDaug == 2 && MC::isElectron(motherPDG) && NumOfPht == 2) return ElMagProc;
1020 
1021  // decay of W,Z and Higgs to lepton with FSR generated by Pythia
1022  if (numOfParents == 1 && numOfDaug == 2 && (MC::isElectron(motherPDG) || MC::isMuon(motherPDG) || MC::isTau(motherPDG)) &&
1023  !(Daug && info && info->Mother() && HepMC::is_same_generator_particle(Daug, info->Mother())) && mothOriVert != nullptr &&
1024  mothOriVert->nIncomingParticles() == 1) {
1025  int itr = 0;
1026  int PartPDG = 0;
1027  const xAOD::TruthVertex* prodVert = mothOriVert;
1028  const xAOD::TruthVertex* Vert = nullptr;
1029  do {
1030  Vert = prodVert;
1031  for (const auto & pin: Vert->particles_in()) {
1032  if (!pin) continue;
1033  PartPDG = abs(pin->pdgId());
1034  prodVert = pin->prodVtx();
1035  if (MC::isZ(PartPDG) || MC::isW(PartPDG) || MC::isHiggs(PartPDG)) foundFSR = true;
1036  }
1037  itr++;
1038  if (itr > 100) { // FIXME Improve loop detection here?
1039  ATH_MSG_WARNING("DefOrigOfPhoton:: infinite while");
1040  break;
1041  }
1042  } while (prodVert != nullptr && abs(motherPDG) == PartPDG);
1043 
1044  if (foundFSR) return FSRPhot;
1045  }
1046 
1047  // Nucl reaction
1048  // gamma+Nuclear=>gamma+Nucl.Fr+Nuclons+pions
1049  // e+Nuclear=>e+gamma+Nucl.Fr+Nuclons+pions
1050  // pi+Nuclear=>gamma+Nucl.Fr+Nuclons+pions
1051 
1052  if ((numOfParents == 1 && (MC::isPhoton(motherPDG) || MC::isElectron(motherPDG)) && numOfDaug > 2 && NumOfNucFr != 0) ||
1053  (numOfParents == 1 && abs(motherPDG) == MC::PIPLUS && numOfDaug > 10 && NumOfNucFr != 0) ||
1054  (numOfParents == 1 && MC::isPhoton(motherPDG) && numOfDaug > 10 && MC::isStable(mother)) ||
1055  (numOfParents == 1 && motherPDG > 1000000000))
1056  return NucReact;
1057 
1058  if (MC::isMuon(motherPDG) && NumOfMu == 0) return Mu;
1059  if (MC::isTau(motherPDG) && NumOfTau == 0) return TauLep;
1060 
1061  if (numOfParents == 1 && mother && mother->status() == 3) return (foundISR)? ISRPhot:UndrPhot;
1062 
1063  //-- to find initial and final state raiation and underline photons
1064  //-- SUSY
1065  if (numOfParents == 1 && (MC::isSMQuark(motherPDG) || MC::isGluon(motherPDG)) &&
1066  (numOfDaug != NumOfPht + NumOfPartons || !MC::Pythia8::isConditionA(mother))) {
1067  for (const auto& pout: partOriVert->particles_out()) {
1068  if (!pout) continue;
1069  if (motherPDG != pout->pdgId()) continue;
1070  const xAOD::TruthVertex* Vrtx = pout->decayVtx();
1071  if (!Vrtx) continue;
1072  if (Vrtx->nOutgoingParticles() != 1 && Vrtx->nIncomingParticles() == 1) continue;
1073  if (!Vrtx->outgoingParticle(0)) continue;
1074  if (Vrtx->outgoingParticle(0)->pdgId() == 91) foundISR = true; // Herwig "cluster"
1075  }
1076  return (foundISR)?ISRPhot:UndrPhot;
1077  }
1078 
1079  //-- to find final state radiation
1080  //-- Exotics
1081 
1082  // FSR from Photos
1083  //-- Exotics- CompHep
1084  if (numOfParents == 2 && ((MC::isElectron(motherPDG) && NumOfEl == 1 && NumOfPos == 1) || (MC::isMuon(motherPDG) && NumOfMu == 2) || (MC::isTau(motherPDG) && NumOfTau == 2))) {
1085  if (abs(partOriVert->incomingParticle(0)->pdgId()) == abs(partOriVert->incomingParticle(1)->pdgId())) return FSRPhot;
1086  }
1087 
1088  if (numOfParents == 2 && NumOfLep == 1 && NumOfNeut == 1 && (MC::isElectron(motherPDG) || abs(motherPDG) == MC::NU_E)) return FSRPhot;
1089 
1090  //-- Exotics - CompHep
1091  if (MC::isElectron(motherPDG) && numOfParents == 1 && numOfDaug == 2 && (NumOfEl == 1 || NumOfPos == 1) && NumOfPht == 1 &&
1092  !( Daug && info && info->Mother() && HepMC::is_same_generator_particle(Daug, info->Mother())) && !HepMC::is_simulation_particle(Daug) && !HepMC::is_simulation_particle(info->Mother()))
1093  return FSRPhot;
1094 
1095  // FSR from Photos
1096  if (MC::isZ(motherPDG) && ((NumOfEl + NumOfPos == 2 || NumOfEl + NumOfPos == 4) || (NumOfMu == 2 || NumOfMu == 4) || (NumOfTau == 2 || NumOfTau == 4)) && NumOfPht > 0) return FSRPhot;
1097 
1098  if (NumOfPht > 0 && (abs(motherPDG) == 9900024 || abs(motherPDG) == 9900012 || abs(motherPDG) == 9900014 || abs(motherPDG) == 9900016)) return FSRPhot; // Left-right symmetric model WBoson || Right-handed NU_E || Right-handed NU_MU || Right-handed NU_TAU
1099 
1100  if (numOfParents == 2 && NumOfLQ == 1) return FSRPhot;
1101 
1102  //--- other process
1103 
1104  if (MC::isZ(motherPDG)) return ZBoson;
1105  if (MC::isW(motherPDG)) {
1106 
1107  if (NumOfLep == 1 && NumOfNeut == 1 && numOfDaug == NumOfLep + NumOfNeut + NumOfPht) return FSRPhot;
1108 
1109  if (mothOriVert != nullptr && mothOriVert->nIncomingParticles() != 0) {
1110  const xAOD::TruthVertex* prodVert = mothOriVert;
1111  const xAOD::TruthParticle* itrP;
1112  do {
1113  itrP = prodVert->incomingParticle(0);
1114  prodVert = itrP->hasProdVtx() ? itrP->prodVtx() : nullptr;
1115  } while (MC::isW(itrP) && prodVert != nullptr);
1116 
1117  if (prodVert && prodVert->nIncomingParticles() == 1 ) {
1118  if ( MC::isTau(itrP)) return TauLep;
1119  if ( MC::isMuon(itrP)) return Mu;
1120  if ( abs(itrP->pdgId()) == 9900012) return NuREle; // Right-handed NU_E
1121  if ( abs(itrP->pdgId()) == 9900014) return NuRMu; // Right-handed NU_MU
1122  if ( abs(itrP->pdgId()) == 9900016) return NuRTau; // Right-handed NU_TAU
1123  }
1124  } else
1125  return WBoson;
1126  }
1127 
1128  // MadGraphPythia ZWW*->lllnulnu
1129  if (numOfParents == 1 && numOfDaug > 4 && (MC::isSMQuark(motherPDG) || MC::isGluon(motherPDG))) {
1130  bool isZboson = false;
1131  bool isWboson = false;
1132  bool skipnext = false;
1133  for (unsigned int ipOut = 0; ipOut + 1 < partOriVert->nOutgoingParticles(); ipOut++) {
1134  if (skipnext) {
1135  skipnext = false;
1136  continue;
1137  }
1138  const xAOD::TruthParticle* theDaug = partOriVert->outgoingParticle(ipOut);
1139  if (!theDaug) continue;
1140  const xAOD::TruthParticle* theNextDaug = nullptr;
1141  for (unsigned int ipOut1 = ipOut + 1; ipOut1 < partOriVert->nOutgoingParticles(); ipOut1++) {
1142  theNextDaug = partOriVert->outgoingParticle(ipOut1);
1143  if (theNextDaug != nullptr) break;
1144  }
1145  if (!theNextDaug) continue;
1146  if (MC::isTau(theDaug) && MC::isTau(theNextDaug)) {
1147  // Zboson
1148  if (thePriPart == theDaug || thePriPart == theNextDaug) {
1149  isZboson = true;
1150  break;
1151  }
1152  skipnext = true;
1153  } else if (MC::isTau(theDaug) && abs(theNextDaug->pdgId()) == MC::NU_TAU) {
1154  // WBoson
1155  if (thePriPart == theDaug || thePriPart == theNextDaug) {
1156  isWboson = true;
1157  break;
1158  }
1159  skipnext = true;
1160  }
1161  }
1162  if (isWboson) return WBoson;
1163  if (isZboson) return ZBoson;
1164  }
1165 
1166  //--Sherpa ZZ,ZW+FSR
1167  if (numOfParents == 4 && (numOfDaug - NumOfPht) == 4 && (NumOfLep + NumOfNeut == 4)) {
1168  if (MC::isSMLepton(partOriVert->incomingParticle(0))&&MC::isSMLepton(partOriVert->incomingParticle(1))
1169  && MC::isSMLepton(partOriVert->incomingParticle(2))&&MC::isSMLepton(partOriVert->incomingParticle(3)))
1170  return FSRPhot;
1171  }
1172 
1173  //--New Sherpa single photon
1174  if (partOriVert == mothOriVert && partOriVert != nullptr) {
1175  int NumOfPhtLoop = 0;
1176  for (const auto *const pout: partOriVert->particles_out()) {
1177  if (!pout) continue;
1178  for (const auto *const pin: partOriVert->particles_in()) {
1179  if (!pin) continue;
1180  if (HepMC::is_same_particle(pout,pin) && MC::isPhoton(pout)) NumOfPhtLoop++;
1181  if (NumOfPhtLoop == 1) return SinglePhot;
1182  }
1183  }
1184  }
1185 
1186  if (MC::isHiggs(motherPDG)) return Higgs;
1187  if (abs(motherPDG) == MC::PI0) return PiZero;
1188  if (MC::isMSSMHiggs(motherPDG)) return HiggsMSSM; // MSSM Higgs bosons
1189  if (MC::isHeavyBoson(motherPDG) || std::abs(motherPDG) == 5100039 ) return HeavyBoson; // Heavy Bosons (Z' Z'' W'+) + KK excited graviton
1190 
1191  if (MC::isSUSY(motherPDG)) return SUSY;
1192  if (MC::isBSM(motherPDG)) return OtherBSM;
1193 
1194  // Pythia8 gamma+jet samples
1195  if (MC::Pythia8::isConditionA(mother) && MC::isStable(thePriPart) && NumOfPht == 1 && numOfDaug == (NumOfPht + NumOfPartons)) return PromptPhot;
1196 
1197  ParticleType pType = defTypeOfHadron(motherPDG);
1198  if ((pType == BBbarMesonPart || pType == CCbarMesonPart) && mothOriVert != nullptr && MC::isHardScatteringVertex(mothOriVert)) isPrompt = true;
1199  return convHadronTypeToOrig(pType, motherPDG);
1200 }
1201 
1204  const xAOD::TruthParticle* thePart,
1205  bool& isPrompt,
1206  MCTruthPartClassifier::Info* infoin) const
1207 {
1209  ATH_MSG_DEBUG("Executing DefOrigOfNeutrino ");
1210 
1211  const int nuFlav = abs(thePart->pdgId());
1212  const xAOD::TruthParticle* thePriPart = MC::findMatching(mcTruthTES, thePart);
1213  if (!thePriPart) return NonDefined;
1214  if (abs(thePriPart->pdgId()) != nuFlav) return NonDefined;
1215 
1216  const xAOD::TruthVertex* partOriVert = thePriPart->hasProdVtx() ? thePriPart->prodVtx() : nullptr;
1217 
1218  //-- to define neutrino outcome status
1219  if (info) info->particleOutCome = NonInteract;
1220 
1222  if (!info) { info = &tmpinfo; }
1223  if (!partOriVert) return NonDefined;
1224 
1225  int numOfParents = -1;
1226  numOfParents = partOriVert->nIncomingParticles();
1227  if (numOfParents > 1) ATH_MSG_DEBUG("DefOrigOfNeutrino:: neutrino has more than one mother ");
1228 
1229  const xAOD::TruthParticle* mother = MC::findMother(thePriPart);
1230  if (info) info->mother = mother;
1231  if (!mother) return NonDefined;
1232  int motherPDG = mother->pdgId();
1233  const xAOD::TruthVertex* mothOriVert = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
1234 
1235  // to resolve Sherpa loop
1236  bool samePart = false;
1237  if (mothOriVert && HepMC::is_same_vertex(mothOriVert,partOriVert)) samePart = true;
1238  //
1239  if ((abs(motherPDG) == nuFlav || MC::isTau(motherPDG) || MC::isW(motherPDG)) && mothOriVert != nullptr &&
1240  !samePart) {
1241  int pPDG(0);
1242  const xAOD::TruthParticle* MotherParent(nullptr);
1243  do {
1244  pPDG = 0;
1245  MotherParent = MC::findMother(mother);
1246  // to prevent Sherpa loop
1247  const xAOD::TruthVertex* mother_prdVtx(nullptr);
1248  const xAOD::TruthVertex* mother_endVtx(nullptr);
1249  if (mother) {
1250  mother_prdVtx = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
1251  mother_endVtx = mother->decayVtx();
1252  }
1253  const xAOD::TruthVertex* parent_prdVtx(nullptr);
1254  const xAOD::TruthVertex* parent_endVtx(nullptr);
1255  if (MotherParent) {
1256  parent_prdVtx = MotherParent->hasProdVtx() ? MotherParent->prodVtx() : nullptr;
1257  parent_endVtx = MotherParent->decayVtx();
1258  }
1259  if (mother_endVtx == parent_prdVtx && mother_prdVtx == parent_endVtx) {
1260  MotherParent = mother;
1261  break;
1262  }
1263  //
1264  if (MotherParent) {
1265  pPDG = MotherParent->pdgId();
1266  }
1267  // to prevent Sherpa loop
1268  if (mother == MotherParent) {
1269  break;
1270  }
1271  if (abs(pPDG) == nuFlav || MC::isTau(pPDG) || MC::isW(pPDG) ) {
1272  mother = MotherParent;
1273  if (info) info->setMotherProperties(mother);
1274  }
1275 
1276  } while ((std::abs(pPDG) == nuFlav || MC::isTau(pPDG) || MC::isW(pPDG)));
1277 
1278  if (std::abs(pPDG) == nuFlav || MC::isTau(pPDG) || MC::isW(pPDG) || MC::isZ(pPDG) || MC::isHiggs(pPDG) ||
1279  MC::isMSSMHiggs(pPDG) || MC::isHeavyBoson(pPDG) || MC::isTop(pPDG) || // MSSM Higgs bosons, Heavy bosons( Z', Z'', W'+)
1280  std::abs(pPDG) == 9900024 || std::abs(pPDG) == 9900012 || std::abs(pPDG) == 9900014 || std::abs(pPDG) == 9900016 || // Left-right symmetric model WBoson || Right-handed NU_E || Right-handed NU_MU || Right-handed NU_TAU
1281  (std::abs(pPDG) < 2000040 && std::abs(pPDG) > 1000001)) { // FIXME This is nearly MC::isSUSY(pPDG), but slightly stricter.
1282  mother = MotherParent;
1283  if (info) info->setMotherProperties(mother);
1284  }
1285  }
1286  //if mother is still nullptr, we have a problem
1287  if (!mother) return NonDefined;
1288 
1289  motherPDG = mother->pdgId();
1290  partOriVert = mother->decayVtx();
1291  mothOriVert = mother->hasProdVtx() ? mother->prodVtx() : nullptr;
1292  numOfParents = partOriVert->nIncomingParticles();
1293  int numOfDaug = partOriVert->nOutgoingParticles();
1294 
1295  if (info) info->setMotherProperties(mother);
1296 
1297  int NumOfPhot(0);
1298  int NumOfquark(0);
1299  int NumOfgluon(0);
1300  int NumOfLQ(0);
1301  int NumOfElNeut(0);
1302  int NumOfMuNeut(0);
1303  int NumOfTauNeut(0);
1304  int NumOfEl(0);
1305  int NumOfMu(0);
1306  int NumOfTau(0);
1307  samePart = false;
1308 
1309  for (const auto& theDaug: partOriVert->particles_out()) {
1310  if (!theDaug) continue;
1311  int DaugType = theDaug->pdgId();
1312  if (MC::isSMQuark(DaugType)) NumOfquark++;
1313  else if (MC::isGluon(DaugType)) NumOfgluon++;
1314  else if (abs(DaugType) == MC::NU_E) NumOfElNeut++;
1315  else if (std::abs(DaugType) == MC::NU_MU) NumOfMuNeut++;
1316  else if (std::abs(DaugType) == MC::NU_TAU) NumOfTauNeut++;
1317  else if (MC::isPhoton(DaugType)) NumOfPhot++;
1318  else if (MC::isElectron(DaugType)) NumOfEl++;
1319  else if (MC::isMuon(DaugType)) NumOfMu++;
1320  else if (MC::isTau(DaugType)) NumOfTau++;
1321  else if (MC::isLeptoQuark(DaugType)) NumOfLQ++;
1322  if (std::abs(DaugType) == std::abs(motherPDG) && theDaug && info && info->Mother() && HepMC::is_same_generator_particle(theDaug,info->Mother())) samePart = true;
1323  }
1324 
1325  // Quark weak decay
1326  if (MC::isQuark(motherPDG) && numOfParents == 1 && numOfDaug == 3 && NumOfquark == 1 && (NumOfEl == 1 || NumOfMu == 1 || NumOfTau == 1)) return QuarkWeakDec;
1327  if (MC::isTop(motherPDG)) return top;
1328 
1329  if (MC::isW(motherPDG) && mothOriVert != nullptr && mothOriVert->nIncomingParticles() != 0) {
1330  const xAOD::TruthVertex* prodVert = mothOriVert;
1331  const xAOD::TruthParticle* ptrPart;
1332  do {
1333  ptrPart = prodVert->incomingParticle(0);
1334  prodVert = ptrPart->hasProdVtx() ? ptrPart->prodVtx() : nullptr;
1335  } while (MC::isW(ptrPart) && prodVert != nullptr);
1336 
1337  if (prodVert && prodVert->nIncomingParticles() == 1) {
1338  if (abs(ptrPart->pdgId()) == 9900012) return NuREle; // Right-handed NU_E
1339  if (abs(ptrPart->pdgId()) == 9900014) return NuRMu; // Right-handed NU_MU
1340  if (abs(ptrPart->pdgId()) == 9900016) return NuRTau; // Right-handed NU_TAU
1341  }
1342  return WBoson;
1343  }
1344  if (MC::isW(motherPDG)) return WBoson;
1345  if (MC::isZ(motherPDG)) return ZBoson;
1346 
1347  //-- Exotics
1348 
1349  // MadGraphPythia ZWW*->lllnulnu or ZWW*->nunulnulnu (don't even know if the latter is generated)
1350  if (numOfParents == 1 && numOfDaug > 4 && (MC::isSMQuark(motherPDG) || MC::isGluon(motherPDG))) {
1351 
1352  const xAOD::TruthParticle* thePartToCheck = thePriPart;
1353  const xAOD::TruthParticle* theMother = thePriPart->hasProdVtx() ? thePriPart->prodVtx()->incomingParticle(0) : nullptr;
1354 
1355  if (MC::isElectron(theMother) && MC::isDecayed(theMother)) thePartToCheck = theMother;
1356  bool isZboson = false;
1357  bool isWboson = false;
1358  bool skipnext = false;
1359 
1360  for (unsigned int ipOut = 0; ipOut + 1 < partOriVert->nOutgoingParticles(); ++ipOut) {
1361  const xAOD::TruthParticle* theDaug = partOriVert->outgoingParticle(ipOut);
1362  if (!theDaug) continue;
1363  const xAOD::TruthParticle* theNextDaug = nullptr;
1364  for (unsigned int ipOut1 = ipOut + 1; ipOut1 < partOriVert->nOutgoingParticles(); ipOut1++) {
1365  theNextDaug = partOriVert->outgoingParticle(ipOut1);
1366  if (theNextDaug != nullptr) break;
1367  }
1368  if (!theNextDaug) continue;
1369 
1370  if (skipnext) {
1371  skipnext = false;
1372  continue;
1373  }
1374 
1375  const int apdgID1 = abs(theDaug->pdgId());
1376  const int apdgID2 = abs(theNextDaug->pdgId());
1377  if (apdgID1 == apdgID2 && MC::isSMNeutrino(apdgID1)) {
1378  // Zboson
1379  if (thePartToCheck == theDaug || thePartToCheck == theNextDaug) {
1380  isZboson = true;
1381  break;
1382  }
1383  skipnext = true;
1384  } else if ((apdgID1 == MC::ELECTRON && apdgID2 == MC::NU_E) ||
1385  (apdgID1 == MC::NU_E && apdgID2 == MC::ELECTRON) ||
1386  (apdgID1 == MC::MUON && apdgID2 == MC::NU_MU) ||
1387  (apdgID1 == MC::NU_MU && apdgID2 == MC::MUON) ||
1388  (apdgID1 == MC::TAU && apdgID2 == MC::NU_TAU) ||
1389  (apdgID1 == MC::NU_TAU && apdgID2 == MC::TAU)
1390  ) {
1391  // WBoson
1392  if (thePartToCheck == theDaug || thePartToCheck == theNextDaug) {
1393  isWboson = true;
1394  break;
1395  }
1396  skipnext = true;
1397  }
1398  }
1399  if (isWboson) return WBoson;
1400  if (isZboson) return ZBoson;
1401  }
1402 
1403  if (numOfParents == 2) {
1404  int pdg1 = partOriVert->incomingParticle(0)->pdgId();
1405  int pdg2 = partOriVert->incomingParticle(1)->pdgId();
1406  //--Sherpa Z->nunu
1407  if ( (numOfDaug - NumOfquark - NumOfgluon) == 2 && (NumOfElNeut == 2 || NumOfMuNeut == 2 || NumOfTauNeut == 2)) return ZBoson;
1408 
1409  //--Sherpa W->enu ??
1410  if ((numOfDaug - NumOfquark - NumOfgluon) == 2 && ((NumOfEl == 1 && NumOfElNeut == 1) || (NumOfMu == 1 && NumOfMuNeut == 1) || (NumOfTau == 1 && NumOfTauNeut == 1))) return WBoson;
1411 
1412  //--Sherpa ZZ,ZW
1413  if ( (numOfDaug - NumOfquark - NumOfgluon) == 4 && (NumOfEl + NumOfMu + NumOfTau + NumOfElNeut + NumOfMuNeut + NumOfTauNeut == 4) &&
1414  (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return DiBoson;
1415 
1416  //--Sherpa VVV -- Note, have to allow for prompt photon radiation or these get lost
1417  if ((numOfDaug - NumOfquark - NumOfgluon - NumOfPhot) == 6 && (NumOfEl + NumOfMu + NumOfTau + NumOfElNeut + NumOfMuNeut + NumOfTauNeut == 6) &&
1418  (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return MultiBoson;
1419  }
1420 
1421  // New Sherpa Z->nunu
1422  if (partOriVert == mothOriVert && partOriVert != nullptr) {
1423  int NumOfLepLoop = 0;
1424  int NumOfNeuLoop = 0;
1425  for (const auto *const pout: partOriVert->particles_out()) {
1426  if (!pout) continue;
1427  for (const auto *const pin: partOriVert->particles_in()) {
1428  if (!pin) continue;
1429  if (HepMC::is_same_particle(pin,pout)) continue;
1430  const int apdgid = abs(pout->pdgId());
1431  if (MC::isSMLepton(apdgid)) {
1432  if (MC::isSMNeutrino(apdgid)) { NumOfNeuLoop++; }
1433  else { NumOfLepLoop++; }
1434  }
1435  }
1436  }
1437  if (NumOfNeuLoop == 2 && NumOfLepLoop == 0) return ZBoson;
1438  if (NumOfNeuLoop == 1 && NumOfLepLoop == 1) return WBoson;
1439  if (NumOfNeuLoop + NumOfLepLoop == 4) return DiBoson;
1440  }
1441 
1442  //-- McAtNLo
1443 
1444  if (MC::isHiggs(motherPDG)) return Higgs;
1445  if (MC::isMSSMHiggs(motherPDG)) return HiggsMSSM; // MSSM Higgs bosons
1446  if (MC::isHeavyBoson(motherPDG)) return HeavyBoson; // Heavy bosons( Z', Z'', W'+)
1447 
1448  if (MC::isTau(motherPDG)) {
1449  ParticleOrigin tauOrig = defOrigOfTau(mcTruthTES, mother, motherPDG, info);
1450  ParticleType tautype = defTypeOfTau(tauOrig);
1451  return (tautype == IsoTau)?tauOrig:TauLep;
1452  }
1453 
1454  if (abs(motherPDG) == 9900024) return WBosonLRSM; // Left-right symmetric model WBoson
1455  if (abs(motherPDG) == 9900012) return NuREle; // Right-handed NU_E
1456  if (abs(motherPDG) == 9900014) return NuRMu; // Right-handed NU_MU
1457  if (abs(motherPDG) == 9900016) return NuRTau; // Right-handed NU_TAU
1458  if (MC::isLeptoQuark(motherPDG) || NumOfLQ != 0) return LQ;
1459  if (MC::isSUSY(motherPDG)) return SUSY;
1460  if (MC::isBSM(motherPDG)) return OtherBSM;
1461 
1462  ParticleType pType = defTypeOfHadron(motherPDG);
1463  if ((pType == BBbarMesonPart || pType == CCbarMesonPart) && mothOriVert != nullptr && MC::isHardScatteringVertex(mothOriVert)) isPrompt = true;
1464 
1465  return convHadronTypeToOrig(pType, motherPDG);
1466 }
grepfile.info
info
Definition: grepfile.py:38
MCTruthClassifier::defOrigOfNeutrino
MCTruthPartClassifier::ParticleOrigin defOrigOfNeutrino(const xAOD::TruthParticleContainer *m_xTruthParticleContainer, const xAOD::TruthParticle *, bool &isPrompt, MCTruthPartClassifier::Info *info) const
Definition: MCTruthClassifierGen.cxx:1203
xAOD::TruthVertex_v1::nOutgoingParticles
size_t nOutgoingParticles() const
Get the number of outgoing particles.
NuRMu
@ NuRMu
Definition: TruthClasses.h:73
temp
Definition: JetEventDict.h:21
IsoPhoton
@ IsoPhoton
Definition: TruthClasses.h:23
MCTruthPartClassifier::defOutComeOfMuon
ParticleOutCome defOutComeOfMuon(T thePart)
Definition: TruthClassifiers.h:213
MCTruthPartClassifier::defTypeOfPhoton
ParticleType defTypeOfPhoton(ParticleOrigin PhotOrig)
Definition: TruthClassifiers.h:115
calibdata.pout
def pout(output, newline=True)
Definition: calibdata.py:130
GenEvent.h
TauLep
@ TauLep
Definition: TruthClasses.h:63
HepMC::is_same_vertex
bool is_same_vertex(const T1 &p1, const T2 &p2)
Method to establish if two particles in the GenEvent actually represent the same vertex.
Definition: MagicNumbers.h:370
OtherBSMParticle
@ OtherBSMParticle
Definition: TruthClasses.h:32
isHeavyBoson
bool isHeavyBoson(const T &p)
APID: Additional "Heavy"/"prime" versions of W and Z bosons (Used in MCTruthClassifier)
Definition: AtlasPID.h:359
NuREle
@ NuREle
Definition: TruthClasses.h:72
Mu
@ Mu
Definition: TruthClasses.h:62
xAOD::TruthVertex_v1::particles_in
std::vector< const TruthParticle * > particles_in() const
Get the incoming particles.
Definition: TruthVertex_v1.cxx:59
WBosonLRSM
@ WBosonLRSM
Definition: TruthClasses.h:71
SingleTau
@ SingleTau
Definition: TruthClasses.h:57
isBSM
bool isBSM(const T &p)
APID: graviton and all Higgs extensions are BSM.
Definition: AtlasPID.h:783
PionDecay
@ PionDecay
Definition: TruthClasses.h:90
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:67
MCTruthPartClassifier::defOutComeOfTau
ParticleOutCome defOutComeOfTau(T thePart)
Definition: TruthClassifiers.h:240
GenVertex.h
ElMagProc
@ ElMagProc
Definition: TruthClasses.h:61
xAODTruthParticleLinkVector::find
ElementLink< xAOD::TruthParticleContainer > find(const HepMcParticleLink &hepMCLink) const
Definition: xAODTruthParticleLink.h:28
isNeutrino
bool isNeutrino(const T &p)
APID: the fourth generation neutrinos are neutrinos.
Definition: AtlasPID.h:201
IsoElectron
@ IsoElectron
Definition: TruthClasses.h:11
SingleElec
@ SingleElec
Definition: TruthClasses.h:54
PhotonConv
@ PhotonConv
Definition: TruthClasses.h:59
ISRPhot
@ ISRPhot
Definition: TruthClasses.h:95
PromptPhot
@ PromptPhot
Definition: TruthClasses.h:93
MCTruthPartClassifier::defTypeOfElectron
ParticleType defTypeOfElectron(ParticleOrigin EleOrig, bool isPrompt)
Definition: TruthClassifiers.h:61
MCTruthClassifier.h
NuRTau
@ NuRTau
Definition: TruthClasses.h:74
MCTruthClassifier::defOrigOfElectron
MCTruthPartClassifier::ParticleOrigin defOrigOfElectron(const xAOD::TruthParticleContainer *xTruthParticleContainer, const xAOD::TruthParticle *, bool &isPrompt, MCTruthPartClassifier::Info *info) const
Definition: MCTruthClassifierGen.cxx:226
UnknownMuon
@ UnknownMuon
Definition: TruthClasses.h:14
ZBoson
@ ZBoson
Definition: TruthClasses.h:67
UnknownTau
@ UnknownTau
Definition: TruthClasses.h:18
Hadron
@ Hadron
Definition: TruthClasses.h:26
xAOD::TruthParticle_v1::pdg_id
int pdg_id() const
PDG ID code.
Definition: TruthParticle_v1.h:52
HepMC::is_same_particle
bool is_same_particle(const T1 &p1, const T2 &p2)
Method to establish if two particles in the GenEvent actually represent the same particle.
Definition: MagicNumbers.h:367
GenParticle.h
WBoson
@ WBoson
Definition: TruthClasses.h:66
isSMQuark
bool isSMQuark(const T &p)
Definition: AtlasPID.h:162
isSMLepton
bool isSMLepton(const T &p)
Definition: AtlasPID.h:183
isGluon
bool isGluon(const T &p)
Definition: AtlasPID.h:346
FSRPhot
@ FSRPhot
Definition: TruthClasses.h:96
isHiggs
bool isHiggs(const T &p)
APID: HIGGS boson is only one particle.
Definition: AtlasPID.h:363
NuclFrag
@ NuclFrag
Definition: TruthClasses.h:28
BBbarMesonPart
@ BBbarMesonPart
Definition: TruthClasses.h:33
MCTruthPartClassifier::defTypeOfMuon
ParticleType defTypeOfMuon(ParticleOrigin MuOrig, bool isPrompt)
Definition: TruthClassifiers.h:80
NonDefined
@ NonDefined
Definition: TruthClasses.h:52
HiggsMSSM
@ HiggsMSSM
Definition: TruthClasses.h:69
MC::isPhysical
bool isPhysical(const T &p)
Identify if the particle is physical, i.e. is stable or decayed.
Definition: HepMCHelpers.h:51
CCbarMesonPart
@ CCbarMesonPart
Definition: TruthClasses.h:35
isQuark
bool isQuark(const T &p)
PDG rule 2: Quarks and leptons are numbered consecutively starting from 1 and 11 respectively; to dot...
Definition: AtlasPID.h:158
Neutrino
@ Neutrino
Definition: TruthClasses.h:27
HepMC::is_same_generator_particle
bool is_same_generator_particle(const T1 &p1, const T2 &p2)
Method to establish if two particles in the GenEvent actually represent the same generated particle.
Definition: MagicNumbers.h:364
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:355
MC::findMatching
T findMatching(C TruthContainer, T p)
Function to find a particle in container.
Definition: HepMCHelpers.h:127
UnknownElectron
@ UnknownElectron
Definition: TruthClasses.h:10
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
ParticleOrigin
ParticleOrigin
Definition: TruthClasses.h:51
PiZero
@ PiZero
Definition: TruthClasses.h:98
xAOD::EgammaHelpers::isElectron
bool isElectron(const xAOD::Egamma *eg)
is the object an electron (not Fwd)
Definition: EgammaxAODHelpers.cxx:12
MultiBoson
@ MultiBoson
Definition: TruthClasses.h:101
isZ
bool isZ(const T &p)
Definition: AtlasPID.h:352
xAOD::TruthParticle_v1
Class describing a truth particle in the MC record.
Definition: TruthParticle_v1.h:37
xAOD::TruthParticle_v1::hasProdVtx
bool hasProdVtx() const
Check for a production vertex on this particle.
Definition: TruthParticle_v1.cxx:74
DiBoson
@ DiBoson
Definition: TruthClasses.h:99
IsoTau
@ IsoTau
Definition: TruthClasses.h:19
LQ
@ LQ
Definition: TruthClasses.h:75
JPsi
@ JPsi
Definition: TruthClasses.h:84
MC::findMother
T findMother(T thePart)
Function to get a mother of particle. MCTruthClassifier legacy.
Definition: HepMCHelpers.h:94
xAOD::TruthVertex_v1::incomingParticle
const TruthParticle_v1 * incomingParticle(size_t index) const
Get one of the incoming particles.
Definition: TruthVertex_v1.cxx:69
IsoMuon
@ IsoMuon
Definition: TruthClasses.h:15
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
MCTruthPartClassifier::defTypeOfHadron
ParticleType defTypeOfHadron(int pdg)
Definition: TruthClassifiers.h:46
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
isSUSY
bool isSUSY(const T &p)
PDG rule 11d Fundamental supersymmetric particles are identified by adding a nonzero n to the particl...
Definition: AtlasPID.h:428
isTau
bool isTau(const T &p)
Definition: AtlasPID.h:197
MC::Pythia8::isConditionA
bool isConditionA(const T &p)
To be understood.
Definition: HepMCHelpers.h:23
xAOD::TruthParticle_v1::decayVtx
const TruthVertex_v1 * decayVtx() const
The decay vertex of this particle.
xAOD::TruthParticle_v1::prodVtx
const TruthVertex_v1 * prodVtx() const
The production vertex of this particle.
Definition: TruthParticle_v1.cxx:80
GetAllXsec.entry
list entry
Definition: GetAllXsec.py:132
HepMC::ConstGenParticlePtr
const GenParticle * ConstGenParticlePtr
Definition: GenParticle.h:38
xAOD::TruthVertex_v1
Class describing a truth vertex in the MC record.
Definition: TruthVertex_v1.h:37
SinglePhot
@ SinglePhot
Definition: TruthClasses.h:56
MCTruthPartClassifier::defOutComeOfElectron
ParticleOutCome defOutComeOfElectron(T thePart)
Definition: TruthClassifiers.h:189
MagicNumbers.h
ReadHandle.h
Handle class for reading from StoreGate.
SUSY
@ SUSY
Definition: TruthClasses.h:77
HeavyBoson
@ HeavyBoson
Definition: TruthClasses.h:70
isHadron
bool isHadron(const T &p)
Definition: AtlasPID.h:324
xAOD::TruthVertex_v1::particles_out
std::vector< const TruthParticle * > particles_out() const
Get the outgoing particles.
Definition: TruthVertex_v1.cxx:64
UndrPhot
@ UndrPhot
Definition: TruthClasses.h:94
SG::VarHandleBase::key
virtual const std::string & key() const override final
Return the StoreGate ID for the referenced object.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleBase.cxx:64
NucReact
@ NucReact
Definition: TruthClasses.h:97
NonInteract
@ NonInteract
Definition: TruthClasses.h:110
SG::ReadHandle::ptr
const_pointer_type ptr()
Dereference the pointer.
QuarkWeakDec
@ QuarkWeakDec
Definition: TruthClasses.h:65
MCTruthClassifier::defOrigOfMuon
MCTruthPartClassifier::ParticleOrigin defOrigOfMuon(const xAOD::TruthParticleContainer *m_xTruthParticleContainer, const xAOD::TruthParticle *, bool &isPrompt, MCTruthPartClassifier::Info *info) const
Definition: MCTruthClassifierGen.cxx:531
isW
bool isW(const T &p)
Definition: AtlasPID.h:355
xAOD::TruthParticle_v1::status
int status() const
Status code.
MC::isStable
bool isStable(const T &p)
Identify if the particle is stable, i.e. has not decayed.
Definition: HepMCHelpers.h:45
SUSYParticle
@ SUSYParticle
Definition: TruthClasses.h:31
isTop
bool isTop(const T &p)
Definition: AtlasPID.h:175
MCTruthPartClassifier::isPrompt
int isPrompt(const unsigned int classify, bool allow_prompt_tau_decays=true)
Definition: TruthClassifiers.h:180
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
MCTruthPartClassifier
Definition: TruthClassifiers.h:12
xAOD::TruthVertex_v1::nIncomingParticles
size_t nIncomingParticles() const
Get the number of incoming particles.
Definition: TruthVertex_v1.cxx:47
MC::isDecayed
bool isDecayed(const T &p)
Identify if the particle decayed.
Definition: HepMCHelpers.h:42
xAOD::EgammaHelpers::isPhoton
bool isPhoton(const xAOD::Egamma *eg)
is the object a photon
Definition: EgammaxAODHelpers.cxx:21
MC::isBeam
bool isBeam(const T &p)
Identify if the particle is beam particle.
Definition: HepMCHelpers.h:39
MCTruthPartClassifier::convHadronTypeToOrig
ParticleOrigin convHadronTypeToOrig(ParticleType pType, int motherPDG)
Definition: TruthClassifiers.h:15
NonPrimary
@ NonPrimary
Definition: TruthClasses.h:29
Higgs
@ Higgs
Definition: TruthClasses.h:68
MCTruthClassifier::defOrigOfTau
MCTruthPartClassifier::ParticleOrigin defOrigOfTau(const xAOD::TruthParticleContainer *m_xTruthParticleContainer, const xAOD::TruthParticle *, int motherPDG, MCTruthPartClassifier::Info *info) const
Definition: MCTruthClassifierGen.cxx:767
isSMNeutrino
bool isSMNeutrino(const T &p)
Definition: AtlasPID.h:204
DecayProducts.h
SingleMuon
@ SingleMuon
Definition: TruthClasses.h:55
MCTruthPartClassifier::defTypeOfTau
ParticleType defTypeOfTau(ParticleOrigin TauOrig)
Definition: TruthClassifiers.h:98
OtherBSM
@ OtherBSM
Definition: TruthClasses.h:78
DalitzDec
@ DalitzDec
Definition: TruthClasses.h:60
top
@ top
Definition: TruthClasses.h:64
MCTruthPartClassifier::defOutComeOfPhoton
ParticleOutCome defOutComeOfPhoton(T thePart)
Definition: TruthClassifiers.h:268
KaonDecay
@ KaonDecay
Definition: TruthClasses.h:91
xAOD::TruthParticle_v1::pdgId
int pdgId() const
PDG ID code.
MCTruthClassifier::particleTruthClassifier
virtual std::pair< MCTruthPartClassifier::ParticleType, MCTruthPartClassifier::ParticleOrigin > particleTruthClassifier(const xAOD::TruthParticle *, MCTruthPartClassifier::Info *info=nullptr) const override final
Definition: MCTruthClassifierGen.cxx:75
ParticleType
ParticleType
Definition: TruthClasses.h:8
MCTruthPartClassifier::Info
Definition: IMCTruthClassifier.h:49
Neutrino
Definition: Neutrino.h:33
xAOD::TruthVertex_v1::outgoingParticle
const TruthParticle_v1 * outgoingParticle(size_t index) const
Get one of the outgoing particles.
Definition: TruthVertex_v1.cxx:119
UnknownPhoton
@ UnknownPhoton
Definition: TruthClasses.h:22
isMSSMHiggs
bool isMSSMHiggs(const T &p)
APID: Additional Higgs bosons for MSSM (Used in MCTruthClassifier)
Definition: AtlasPID.h:367
BremPhot
@ BremPhot
Definition: TruthClasses.h:92
GenParticle
@ GenParticle
Definition: TruthClasses.h:30
MC::isHardScatteringVertex
bool isHardScatteringVertex(T pVert)
Function to classify the vertex as hard scattering vertex.
Definition: HepMCHelpers.h:168
MCTruthClassifier::defOrigOfPhoton
MCTruthPartClassifier::ParticleOrigin defOrigOfPhoton(const xAOD::TruthParticleContainer *m_xTruthParticleContainer, const xAOD::TruthParticle *, bool &isPrompt, MCTruthPartClassifier::Info *info) const
Definition: MCTruthClassifierGen.cxx:944
isMuon
bool isMuon(const T &p)
Definition: AtlasPID.h:194
isLeptoQuark
bool isLeptoQuark(const T &p)
PDG rule 11c: “One-of-a-kind” exotic particles are assigned numbers in the range 41–80.
Definition: AtlasPID.h:381