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