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 
12 #include "MCTruthClassifier/MCTruthClassifier.h"
13 #include "AsgDataHandles/ReadHandle.h"
14 #include "TruthUtils/MagicNumbers.h"
15 #ifndef XAOD_ANALYSIS
16 #include "AtlasHepMC/GenEvent.h"
17 #include "AtlasHepMC/GenVertex.h"
18 #include "AtlasHepMC/GenParticle.h"
19 #endif
20 #include "TruthUtils/DecayProducts.h"
21 using namespace MCTruthPartClassifier;
22 using std::abs;
23 
24 #ifndef XAOD_ANALYSIS
25 std::pair<ParticleType, ParticleOrigin>
26 MCTruthClassifier::particleTruthClassifier(const HepMcParticleLink& theLink, MCTruthPartClassifier::Info* info /*= nullptr*/) const {
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>
42 MCTruthClassifier::particleTruthClassifier(HepMC::ConstGenParticlePtr theGenPart, MCTruthPartClassifier::Info* info /*= nullptr*/) const {
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)) {
63  ATH_MSG_DEBUG(
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>
75 MCTruthClassifier::particleTruthClassifier(const xAOD::TruthParticle* thePart, MCTruthPartClassifier::Info* infoin /*= nullptr*/) const {
76  MCTruthPartClassifier::Info* info = infoin;
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();
86  MCTruthPartClassifier::Info tmpinfo;
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 
226 ParticleOrigin MCTruthClassifier::defOrigOfElectron(const xAOD::TruthParticleContainer* mcTruthTES,
227  const xAOD::TruthParticle* thePart,
228  bool& isPrompt,
229  MCTruthPartClassifier::Info* infoin) const
230 {
231  MCTruthPartClassifier::Info* info = infoin;
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);
242  MCTruthPartClassifier::Info tmpinfo;
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 ((abs(motherPDG) == 13 || abs(motherPDG) == 15 || abs(motherPDG) == 24) && mothOriVert != nullptr && !samePart) {
271  long 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 (abs(pPDG) == 13 || abs(pPDG) == 15 || abs(pPDG) == 24) mother = MotherParent;
298 
299  } while ((abs(pPDG) == 13 || abs(pPDG) == 15 || abs(pPDG) == 24));
300 
301  if (abs(pPDG) == 13 || abs(pPDG) == 15 || abs(pPDG) == 24 || abs(pPDG) == 23 || abs(pPDG) == 25 ||
302  abs(pPDG) == 35 || abs(pPDG) == 36 || abs(pPDG) == 37 || abs(pPDG) == 32 || abs(pPDG) == 33 ||
303  abs(pPDG) == 34 || abs(pPDG) == 6 || abs(pPDG) == 9900024 || abs(pPDG) == 9900012 || abs(pPDG) == 9900014 ||
304  abs(pPDG) == 9900016 || (abs(pPDG) < 2000040 && abs(pPDG) > 1000001))
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 (abs(DaugType) < 7) NumOfquark++;
335  if (abs(DaugType) == 21) NumOfgluon++;
336  if (abs(DaugType) == 12) NumOfElNeut++;
337  if (abs(DaugType) == 14) NumOfMuNeut++;
338  if (DaugType == 22) NumOfPhot++;
339  if (DaugType == 11) NumOfEl++;
340  if (DaugType == -11) NumOfPos++;
341  if (DaugType == 13) NumOfMuMin++;
342  if (DaugType == -13) NumOfMuPl++;
343  if (abs(DaugType) == 15) NumOfTau++;
344  if (abs(DaugType) == 16) NumOfTauNeut++;
345  if (abs(DaugType) == 42) 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  (motherPDG == 22 || abs(motherPDG) == 11 || abs(motherPDG) == 13 || abs(motherPDG) == 211) &&
349  (DaugType > 1000000000 || DaugType == 0 || DaugType == 2212 || DaugType == 2112 || abs(DaugType) == 211 ||
350  abs(DaugType) == 111))
351  NumOfNucFr++;
352  }
353 
354  if (motherPDG == 22 && 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 ((motherPDG == 22 && numOfDaug == 2 && NumOfEl == 1 && NumOfPos == 1) || (motherPDG == 22 && numOfDaug == 1 && (NumOfEl == 1 || NumOfPos == 1))) return PhotonConv;
363 
364  // e,gamma,pi+Nuclear->NuclearFragments+nuclons+e
365  if ((numOfParents == 1 && (abs(motherPDG) == 22 || abs(motherPDG) == 11 || abs(motherPDG) == 15)) && numOfDaug > 1 && NumOfNucFr != 0) return ElMagProc;
366 
367  if (numOfParents == 1 && abs(motherPDG) == 211 && numOfDaug > 2 && NumOfNucFr != 0) return ElMagProc;
368 
369  // nuclear photo fission
370  if (motherPDG == 22 && numOfDaug > 4 && NumOfNucFr != 0) return ElMagProc;
371 
372  // unknown process el(pos)->el+pos??
373  if (abs(motherPDG) == 11 && numOfDaug == 2 && NumOfEl == 1 && NumOfPos == 1) return ElMagProc;
374 
375  // unknown process el->el+el??
376  if (motherPDG == 11 && numOfDaug == 2 && NumOfEl == 2 && NumOfPos == 0) return ElMagProc;
377 
378  // unknown process pos->pos+pos??
379  if (motherPDG == -11 && numOfDaug == 2 && NumOfEl == 0 && NumOfPos == 2) return ElMagProc;
380 
381  // unknown process pos/el->pos/el??
382  if (abs(motherPDG) == 11 && !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) && abs(motherPDG) != 11 && samePart) return ElMagProc;
387 
388  if ((motherPDG == 111 && numOfDaug == 3 && NumOfPhot == 1 && NumOfEl == 1 && NumOfPos == 1) ||
389  (motherPDG == 111 && numOfDaug == 4 && NumOfPhot == 0 && NumOfEl == 2 && NumOfPos == 2))
390  return DalitzDec;
391 
392  // Quark weak decay
393  if (abs(motherPDG) < 7 && numOfParents == 1 && numOfDaug == 3 && NumOfquark == 1 && NumOfElNeut == 1) return QuarkWeakDec;
394 
395  if (abs(motherPDG) == 13 && NumOfNucFr != 0) return ElMagProc;
396 
397  if (abs(motherPDG) == 6) 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;
410  if (abs(ptrPart->pdgId()) == 9900014) return NuRMu;
411  if (abs(ptrPart->pdgId()) == 9900016) return NuRTau;
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 && (abs(motherPDG) < 7 || motherPDG == 21)) {
420 
421  const xAOD::TruthParticle* thePartToCheck = thePriPart;
422  const xAOD::TruthParticle* theMother = thePriPart->hasProdVtx() ? thePriPart->prodVtx()->incomingParticle(0) : nullptr;
423  if (theMother != nullptr && abs(theMother->pdgId()) == 11 && 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 (abs(theDaug->pdgId()) == 11 && abs(theNextDaug->pdgId()) == 11) {
444  // Zboson
445  if (thePartToCheck == theDaug || thePartToCheck == theNextDaug) {
446  isZboson = true;
447  break;
448  }
449  skipnext = true;
450  } else if (abs(theDaug->pdgId()) == 11 && abs(theNextDaug->pdgId()) == 12) {
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  int pdg1 = partOriVert->incomingParticle(0)->pdgId();
464  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()) == 12) 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 (abs(motherPDG) == 25) return Higgs;
506 
507  if (abs(motherPDG) == 35 || abs(motherPDG) == 36 || abs(motherPDG) == 37) return HiggsMSSM;
508 
509  if (abs(motherPDG) == 32 || abs(motherPDG) == 33 || abs(motherPDG) == 34) return HeavyBoson;
510 
511  if (abs(motherPDG) == 13) return Mu;
512  if (abs(motherPDG) == 15) {
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;
519  if (abs(motherPDG) == 9900012) return NuREle;
520  if (abs(motherPDG) == 9900014) return NuRMu;
521  if (abs(motherPDG) == 9900016) return NuRTau;
522  if (abs(motherPDG) == 42 || 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 
531 ParticleOrigin MCTruthClassifier::defOrigOfMuon(const xAOD::TruthParticleContainer* mcTruthTES,
532  const xAOD::TruthParticle* thePart,
533  bool& isPrompt,
534  MCTruthPartClassifier::Info* infoin) const
535 {
536  MCTruthPartClassifier::Info* info = infoin;
537  ATH_MSG_DEBUG("Executing DefOrigOfMuon ");
538 
539  const xAOD::TruthParticle* thePriPart = MC::findMatching(mcTruthTES, thePart);
540  if (!thePriPart) return NonDefined;
541  if (abs(thePriPart->pdgId()) != 13) 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 
548  MCTruthPartClassifier::Info tmpinfo;
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  long 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 (abs(pPDG) == 13 || abs(pPDG) == 15 || abs(pPDG) == 24) {
598  mother = MotherParent;
599  if (info) info->setMotherProperties(mother);
600  }
601  }
602  } while ((abs(pPDG) == 13 || abs(pPDG) == 15 || abs(pPDG) == 24));
603 
604  if (abs(pPDG) == 15 || abs(pPDG) == 24 || abs(pPDG) == 23 || abs(pPDG) == 25 || abs(pPDG) == 35 ||
605  abs(pPDG) == 36 || abs(pPDG) == 37 || abs(pPDG) == 32 || abs(pPDG) == 33 || abs(pPDG) == 34 || abs(pPDG) == 6 ||
606  abs(pPDG) == 9900024 || abs(pPDG) == 9900012 || abs(pPDG) == 9900014 || abs(pPDG) == 9900016 ||
607  (abs(pPDG) < 2000040 && abs(pPDG) > 1000001)) {
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  int NumOfPhot = DP.pd(22);
621  int NumOfEl = DP.pd(11);
622  int NumOfPos = DP.pd(-11);
623  int NumOfElNeut = DP.apd(12);
624  int NumOfMuNeut = DP.apd(14);
625  int NumOfLQ = DP.apd(42);
626  int NumOfquark = DP.apd({1,2,3,4,5,6});
627  int NumOfgluon = DP.apd(21);
628  int NumOfMuPl = DP.pd(-13);
629  int NumOfMuMin = DP.pd(13);
630  int NumOfTau = DP.apd(15);
631  int NumOfTauNeut = DP.apd(16);
632  if (std::abs(motherPDG) == 211 && numOfDaug == 2 && NumOfMuNeut == 1) return PionDecay;
633  if (std::abs(motherPDG) == 321 && 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 (std::abs(motherPDG) == 6) return top;
641  // Quark weak decay
642  if (abs(motherPDG) < 7 && 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;
654  if (abs(itrP->pdgId()) == 9900014) return NuRMu;
655  if (abs(itrP->pdgId()) == 9900016) return NuRTau;
656  }
657  return WBoson;
658  }
659  if (MC::isW(motherPDG)) return WBoson;
660  if (MC::isZ(motherPDG)) return ZBoson;
661  if (motherPDG == 22 && numOfDaug == 2 && NumOfMuMin == 1 && NumOfMuPl == 1) return PhotonConv;
662  //-- Exotics
663 
664  // MadGraphPythia ZWW*->lllnulnu
665 
666  if (numOfParents == 1 && numOfDaug > 4 && (abs(motherPDG) < 7 || motherPDG == 21)) {
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 (abs(theDaug->pdgId()) == 13 && abs(theNextDaug->pdgId()) == 13) {
684  // Zboson
685  if (thePriPart == theDaug || thePriPart == theNextDaug) {
686  isZboson = true;
687  break;
688  }
689  skipnext = true;
690  } else if (abs(theDaug->pdgId()) == 13 && abs(theNextDaug->pdgId()) == 14) {
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  int pdg1 = partOriVert->incomingParticle(0)->pdgId();
704  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 (std::abs(pout->pdg_id()) == 13) NumOfMuLoop++;
734  if (std::abs(pout->pdg_id()) == 14) 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 (abs(motherPDG) == 25) return Higgs;
748 
749  if (abs(motherPDG) == 35 || abs(motherPDG) == 36 || abs(motherPDG) == 37) return HiggsMSSM;
750 
751  if (abs(motherPDG) == 32 || abs(motherPDG) == 33 || abs(motherPDG) == 34) return HeavyBoson;
752 
753  if (abs(motherPDG) == 9900024) return WBosonLRSM;
754  if (abs(motherPDG) == 9900012) return NuREle;
755  if (abs(motherPDG) == 9900014) return NuRMu;
756  if (abs(motherPDG) == 9900016) return NuRTau;
757  if (abs(motherPDG) == 42 || 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 
767 ParticleOrigin MCTruthClassifier::defOrigOfTau(const xAOD::TruthParticleContainer* mcTruthTES,
768  const xAOD::TruthParticle* thePart,
769  int motherPDG,
770  MCTruthPartClassifier::Info* infoin) const
771 {
772  MCTruthPartClassifier::Info* info = infoin;
773  ATH_MSG_DEBUG("Executing DefOrigOfTau ");
774  const xAOD::TruthParticle* thePriPart = MC::findMatching(mcTruthTES, thePart);
775  if (!thePriPart) return NonDefined;
776  if (abs(thePriPart->pdgId()) != 15) 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 
783  MCTruthPartClassifier::Info tmpinfo;
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  long pPDG(0);
803 
804  if (MotherParent) {//MotherParent checked here...
805  pPDG = MotherParent->pdgId();
806  if (abs(pPDG) == 6) {
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  int numOfDaug = DP.size();
822  int NumOfPhot = DP.pd(22);
823  int NumOfEl = DP.pd(11);
824  int NumOfPos = DP.pd(-11);
825  int NumOfElNeut = DP.apd(12);
826  int NumOfMuNeut = DP.apd(14);
827  /* int NumOfLQ = DP.apd(42); */
828  int NumOfquark = DP.apd({1,2,3,4,5,6});
829  int NumOfgluon = DP.apd(21);
830  int NumOfMuPl = DP.pd(-13);
831  int NumOfMuMin = DP.pd(13);
832  int NumOfTau = DP.apd(15);
833  int NumOfTauNeut = DP.pd(16);
834 
835  if (abs(motherPDG) == 6) 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;
846  if (abs(itrP->pdgId()) == 9900014) return NuRMu;
847  if (abs(itrP->pdgId()) == 9900016) return NuRTau;
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 && (abs(motherPDG) < 7 || motherPDG == 21)) {
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 (abs(theDaug->pdgId()) == 15 && abs(theNextDaug->pdgId()) == 15) {
873  // Zboson
874  if (thePriPart == theDaug || thePriPart == theNextDaug) {
875  isZboson = true;
876  break;
877  }
878  skipnext = true;
879  } else if (abs(theDaug->pdgId()) == 15 && abs(theNextDaug->pdgId()) == 16) {
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  int pdg1 = partOriVert->incomingParticle(0)->pdgId();
893  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 (std::abs(pout->pdgId()) == 15) NumOfTauLoop++;
921  if (std::abs(pout->pdgId()) == 16) 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 (abs(motherPDG) == 25) return Higgs;
932  if (abs(motherPDG) == 35 || abs(motherPDG) == 36 || abs(motherPDG) == 37) return HiggsMSSM;
933  if (abs(motherPDG) == 32 || abs(motherPDG) == 33 || abs(motherPDG) == 34) return HeavyBoson;
934  if (abs(motherPDG) == 9900024) return WBosonLRSM;
935  if (abs(motherPDG) == 9900016) return NuRTau;
936  if (MC::isSUSY(motherPDG)) return SUSY;
937  if (MC::isBSM(motherPDG)) return OtherBSM;
938  if (abs(motherPDG) == 443) return JPsi;
939 
940  ParticleType pType = defTypeOfHadron(motherPDG);
941  return convHadronTypeToOrig(pType, motherPDG);
942 }
943 
944 ParticleOrigin MCTruthClassifier::defOrigOfPhoton(const xAOD::TruthParticleContainer* mcTruthTES,
945  const xAOD::TruthParticle* thePart,
946  bool& isPrompt,
947  MCTruthPartClassifier::Info* infoin) const
948 {
949  if (!thePart) return NonDefined;
950  if (!mcTruthTES) return NonDefined;
951  MCTruthPartClassifier::Info* info = infoin;
952  ATH_MSG_DEBUG("Executing DefOrigOfPhoton ");
953 
954  MCTruthPartClassifier::Info tmpinfo;
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 (abs(thePriPart->pdgId()) != 22) 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  long 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 && (motherPDG == 22 || abs(motherPDG) == 11 || abs(motherPDG) == 211) &&
1000  (DaugType > 1000000000 || DaugType == 0 || DaugType == 2212 || DaugType == 2112))
1001  NumOfNucFr++;
1002  if (DaugType == 22) NumOfPht++;
1003  if (DaugType == 11) NumOfEl++;
1004  if (DaugType == -11) NumOfPos++;
1005  if (abs(DaugType) == 13) NumOfMu++;
1006  if (abs(DaugType) == 15) NumOfTau++;
1007  if (abs(DaugType) == 42) NumOfLQ++;
1008  if (abs(DaugType) == 11 || abs(DaugType) == 13 || abs(DaugType) == 15) NumOfLep++;
1009  if (abs(DaugType) == 12 || abs(DaugType) == 14 || abs(DaugType) == 16) NumOfNeut++;
1010  if (abs(DaugType) < 11 || (abs(DaugType) > 16 && abs(DaugType) < 43 && abs(DaugType) != 22)) NumOfPartons++;
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 && abs(motherPDG) == 11 && 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 && (abs(motherPDG) == 11 || abs(motherPDG) == 13 || abs(motherPDG) == 15) &&
1023  !(Daug && info && info->Mother() && HepMC::is_same_generator_particle(Daug, info->Mother())) && mothOriVert != nullptr &&
1024  mothOriVert->nIncomingParticles() == 1) {
1025  int itr = 0;
1026  long 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 (PartPDG == 23 || PartPDG == 24 || PartPDG == 25) foundFSR = true;
1036  }
1037  itr++;
1038  if (itr > 100) {
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 && (abs(motherPDG) == 22 || abs(motherPDG) == 11) && numOfDaug > 2 && NumOfNucFr != 0) ||
1053  (numOfParents == 1 && abs(motherPDG) == 211 && numOfDaug > 10 && NumOfNucFr != 0) ||
1054  (numOfParents == 1 && motherPDG == 22 && 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 && (abs(motherPDG) < 7 || motherPDG == 21) &&
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;
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 && ((abs(motherPDG) == 11 && NumOfEl == 1 && NumOfPos == 1) || (abs(motherPDG) == 13 && NumOfMu == 2) || (abs(motherPDG) == 15 && 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 && (abs(motherPDG) == 11 || abs(motherPDG) == 12)) return FSRPhot;
1089 
1090  //-- Exotics - CompHep
1091  if (abs(motherPDG) == 11 && 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;
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;
1121  if ( abs(itrP->pdgId()) == 9900014) return NuRMu;
1122  if ( abs(itrP->pdgId()) == 9900016) return NuRTau;
1123  }
1124  } else
1125  return WBoson;
1126  }
1127 
1128  // MadGraphPythia ZWW*->lllnulnu
1129  if (numOfParents == 1 && numOfDaug > 4 && (abs(motherPDG) < 7 || motherPDG == 21)) {
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 (abs(theDaug->pdgId()) == 15 && abs(theNextDaug->pdgId()) == 15) {
1147  // Zboson
1148  if (thePriPart == theDaug || thePriPart == theNextDaug) {
1149  isZboson = true;
1150  break;
1151  }
1152  skipnext = true;
1153  } else if (abs(theDaug->pdgId()) == 15 && abs(theNextDaug->pdgId()) == 16) {
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 (abs(motherPDG) == 25) return Higgs;
1187  if (abs(motherPDG) == 111) return PiZero;
1188  if (abs(motherPDG) == 35 || abs(motherPDG) == 36 || abs(motherPDG) == 37) return HiggsMSSM;
1189  if (abs(motherPDG) == 32 || abs(motherPDG) == 33 || abs(motherPDG) == 34 || abs(motherPDG) == 5100039 ) return HeavyBoson;
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 
1202 ParticleOrigin
1203 MCTruthClassifier::defOrigOfNeutrino(const xAOD::TruthParticleContainer* mcTruthTES,
1204  const xAOD::TruthParticle* thePart,
1205  bool& isPrompt,
1206  MCTruthPartClassifier::Info* infoin) const
1207 {
1208  MCTruthPartClassifier::Info* info = infoin;
1209  ATH_MSG_DEBUG("Executing DefOrigOfNeutrino ");
1210 
1211  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 
1221  MCTruthPartClassifier::Info tmpinfo;
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 || abs(motherPDG) == 15 || MC::isW(motherPDG)) && mothOriVert != nullptr &&
1240  !samePart) {
1241  long 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 || abs(pPDG) == 15 || abs(pPDG) == 24 ) {
1272  mother = MotherParent;
1273  if (info) info->setMotherProperties(mother);
1274  }
1275 
1276  } while ((std::abs(pPDG) == nuFlav || std::abs(pPDG) == 15 || std::abs(pPDG) == 24));
1277 
1278  if (std::abs(pPDG) == nuFlav || std::abs(pPDG) == 15 || std::abs(pPDG) == 24 || std::abs(pPDG) == 23 || std::abs(pPDG) == 25 ||
1279  std::abs(pPDG) == 35 || std::abs(pPDG) == 36 || std::abs(pPDG) == 37 || std::abs(pPDG) == 32 || std::abs(pPDG) == 33 ||
1280  std::abs(pPDG) == 34 || std::abs(pPDG) == 6 || std::abs(pPDG) == 9900024 || std::abs(pPDG) == 9900012 || std::abs(pPDG) == 9900014 ||
1281  std::abs(pPDG) == 9900016 || (std::abs(pPDG) < 2000040 && std::abs(pPDG) > 1000001)) {
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 (abs(DaugType) < 7) NumOfquark++;
1313  if (abs(DaugType) == 21) NumOfgluon++;
1314  if (abs(DaugType) == 12) NumOfElNeut++;
1315  if (std::abs(DaugType) == 14) NumOfMuNeut++;
1316  if (std::abs(DaugType) == 16) NumOfTauNeut++;
1317  if (DaugType == 22) NumOfPhot++;
1318  if (std::abs(DaugType) == 11) NumOfEl++;
1319  if (std::abs(DaugType) == 13) NumOfMu++;
1320  if (std::abs(DaugType) == 15) NumOfTau++;
1321  if (std::abs(DaugType) == 42) 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 (std::abs(motherPDG) == 6) 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;
1339  if (abs(ptrPart->pdgId()) == 9900014) return NuRMu;
1340  if (abs(ptrPart->pdgId()) == 9900016) return NuRTau;
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 && (abs(motherPDG) < 7 || motherPDG == 21)) {
1351 
1352  const xAOD::TruthParticle* thePartToCheck = thePriPart;
1353  const xAOD::TruthParticle* theMother = thePriPart->hasProdVtx() ? thePriPart->prodVtx()->incomingParticle(0) : nullptr;
1354 
1355  if (abs(theMother->pdgId()) == 11 && 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  int apdgID1 = abs(theDaug->pdgId());
1376  int apdgID2 = abs(theNextDaug->pdgId());
1377  if ((apdgID1 == 12 && apdgID2 == 12) || (apdgID1 == 14 && apdgID2 == 14) || (apdgID1 == 16 && apdgID2 == 16)) {
1378  // Zboson
1379  if (thePartToCheck == theDaug || thePartToCheck == theNextDaug) {
1380  isZboson = true;
1381  break;
1382  }
1383  skipnext = true;
1384  } else if ((apdgID1 == 11 && apdgID2 == 12) || (apdgID1 == 14 && apdgID2 == 14) ||
1385  (apdgID1 == 16 && apdgID2 == 16)) {
1386  // WBoson
1387  if (thePartToCheck == theDaug || thePartToCheck == theNextDaug) {
1388  isWboson = true;
1389  break;
1390  }
1391  skipnext = true;
1392  }
1393  }
1394  if (isWboson) return WBoson;
1395  if (isZboson) return ZBoson;
1396  }
1397 
1398  if (numOfParents == 2) {
1399  int pdg1 = partOriVert->incomingParticle(0)->pdgId();
1400  int pdg2 = partOriVert->incomingParticle(1)->pdgId();
1401  //--Sherpa Z->nunu
1402  if ( (numOfDaug - NumOfquark - NumOfgluon) == 2 && (NumOfElNeut == 2 || NumOfMuNeut == 2 || NumOfTauNeut == 2)) return ZBoson;
1403 
1404  //--Sherpa W->enu ??
1405  if ((numOfDaug - NumOfquark - NumOfgluon) == 2 && ((NumOfEl == 1 && NumOfElNeut == 1) || (NumOfMu == 1 && NumOfMuNeut == 1) || (NumOfTau == 1 && NumOfTauNeut == 1))) return WBoson;
1406 
1407  //--Sherpa ZZ,ZW
1408  if ( (numOfDaug - NumOfquark - NumOfgluon) == 4 && (NumOfEl + NumOfMu + NumOfTau + NumOfElNeut + NumOfMuNeut + NumOfTauNeut == 4) &&
1409  (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return DiBoson;
1410 
1411  //--Sherpa VVV -- Note, have to allow for prompt photon radiation or these get lost
1412  if ((numOfDaug - NumOfquark - NumOfgluon - NumOfPhot) == 6 && (NumOfEl + NumOfMu + NumOfTau + NumOfElNeut + NumOfMuNeut + NumOfTauNeut == 6) &&
1413  (MC::isQuark(pdg1)||MC::isGluon(pdg1)) && (MC::isQuark(pdg2)||MC::isGluon(pdg2))) return MultiBoson;
1414  }
1415 
1416  // New Sherpa Z->nunu
1417  if (partOriVert == mothOriVert && partOriVert != nullptr) {
1418  int NumOfLepLoop = 0;
1419  int NumOfNeuLoop = 0;
1420  for (const auto *const pout: partOriVert->particles_out()) {
1421  if (!pout) continue;
1422  for (const auto *const pin: partOriVert->particles_in()) {
1423  if (!pin) continue;
1424  if (HepMC::is_same_particle(pin,pout)) continue;
1425  int apdgid = abs(pout->pdgId());
1426  if (apdgid == 12 || apdgid == 14 || apdgid == 16) NumOfNeuLoop++;
1427  if (apdgid == 11 || apdgid == 13 || apdgid == 15) NumOfLepLoop++;
1428  }
1429  }
1430  if (NumOfNeuLoop == 2 && NumOfLepLoop == 0) return ZBoson;
1431  if (NumOfNeuLoop == 1 && NumOfLepLoop == 1) return WBoson;
1432  if (NumOfNeuLoop + NumOfLepLoop == 4) return DiBoson;
1433  }
1434 
1435  //-- McAtNLo
1436 
1437  if (abs(motherPDG) == 25) return Higgs;
1438  if (abs(motherPDG) == 35 || abs(motherPDG) == 36 || abs(motherPDG) == 37) return HiggsMSSM;
1439  if (abs(motherPDG) == 32 || abs(motherPDG) == 33 || abs(motherPDG) == 34) return HeavyBoson;
1440 
1441  if (abs(motherPDG) == 15) {
1442  ParticleOrigin tauOrig = defOrigOfTau(mcTruthTES, mother, motherPDG, info);
1443  ParticleType tautype = defTypeOfTau(tauOrig);
1444  return (tautype == IsoTau)?tauOrig:TauLep;
1445  }
1446 
1447  if (abs(motherPDG) == 9900024) return WBosonLRSM;
1448  if (abs(motherPDG) == 9900012) return NuREle;
1449  if (abs(motherPDG) == 9900014) return NuRMu;
1450  if (abs(motherPDG) == 9900016) return NuRTau;
1451  if (abs(motherPDG) == 42 || NumOfLQ != 0) return LQ;
1452  if (MC::isSUSY(motherPDG)) return SUSY;
1453  if (MC::isBSM(motherPDG)) return OtherBSM;
1454 
1455  ParticleType pType = defTypeOfHadron(motherPDG);
1456  if ((pType == BBbarMesonPart || pType == CCbarMesonPart) && mothOriVert != nullptr && MC::isHardScatteringVertex(mothOriVert)) isPrompt = true;
1457 
1458  return convHadronTypeToOrig(pType, motherPDG);
1459 }
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:357
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:59
WBosonLRSM
@ WBosonLRSM
Definition: TruthClasses.h:71
SingleTau
@ SingleTau
Definition: TruthClasses.h:57
isBSM
bool isBSM(const T &p)
Definition: AtlasPID.h:224
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:161
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:354
GenParticle.h
WBoson
@ WBoson
Definition: TruthClasses.h:66
isSMLepton
bool isSMLepton(const T &p)
Definition: AtlasPID.h:143
isGluon
bool isGluon(const T &p)
Definition: AtlasPID.h:167
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
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:122
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:351
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:342
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:173
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: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:210
isTau
bool isTau(const T &p)
Definition: AtlasPID.h:157
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:218
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:176
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
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
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
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:154
Generated on Thu Nov 7 2024 21:20:34 for ATLAS Offline Software by   doxygen 1.8.18