22 #include "GaudiKernel/SystemOfUnits.h"
24 #include "GaudiKernel/ThreadLocalContext.h"
34 const IInterface*
p ) :
36 m_eventsKey(
"TruthEvents"),
41 declareProperty (
"EventsKey",
43 "TruthEvent container name");
45 declareProperty (
"WritePartons",
49 declareProperty (
"WriteHadrons",
53 declareProperty (
"WriteBHadrons",
57 declareProperty (
"WriteCHadrons",
61 declareProperty (
"WriteGeant",
63 "Keep Geant particles?");
65 declareProperty (
"GeantPhotonPtThresh",
67 "Write Geant photons with Pt above this threshold. "
68 "Set to < 0 to not write any.");
70 declareProperty (
"WriteTauHad",
72 "Keep hadronic taus?");
74 declareProperty (
"PartonPtThresh",
76 "Write partons with Pt above this threshold.");
78 declareProperty (
"WriteBSM",
80 "Keep BSM particles?");
82 declareProperty (
"WriteBosons",
86 declareProperty (
"PhotonPtThresh",
88 "Photon pt cut with WriteBosons");
90 declareProperty (
"WriteBSMProducts",
92 "Keep BSM particle decay products?");
94 declareProperty (
"WriteBosonProducts",
96 "Keep boson decay products?");
98 declareProperty (
"WriteTopAndDecays",
100 "Keep top partons and immediate decay products?");
102 declareProperty(
"WriteEverything",
104 "Keep absolutely everything (overrides all other flags)");
106 declareProperty(
"WriteAllLeptons",
108 "Keep absolutely all leptons");
110 declareProperty(
"WriteLeptonsNotFromHadrons",
112 "Keep leptons not from hadron decays");
114 declareProperty(
"WriteAllStable",
116 "Keep all stable particles");
118 declareProperty(
"WriteNotPhysical",
120 "Save also non-physical particles");
122 declareProperty(
"WriteFirstN",
124 "Keep first N particles in record");
126 declareProperty(
"PreserveDescendants",
128 "Preserve descendants of retained particles");
130 declareProperty(
"PreserveGeneratorDescendants",
132 "Preserve descendants of retained particles excluding Geant particles");
134 declareProperty(
"PreserveAncestors",
136 "Preserve ancestors of retained particles");
138 declareProperty (
"PreserveParentsSiblingsChildren",
140 "Preserve the parents, siblings and children of retained particles");
142 declareProperty (
"PreserveHadronizationVertices",
144 "Preserve hadronization vertices with parents/children.");
146 declareProperty (
"WritettHFHadrons",
148 "Keep tt+HF hadrons?");
150 declareProperty (
"PDGIDsToKeep",
152 "List of PDG IDs to always keep");
154 declareProperty (
"LongLivedPDGIDs",
156 "List of PDG IDs of long lived particles so that one can store their children");
166 if (m_preserveDescendants && m_preserveGeneratorDescendants) {
167 ATH_MSG_FATAL(
"You are asking to keep both all descendants, and only those from the event generator. Please check your job options.");
168 return StatusCode::FAILURE;
170 if (m_preserveImmediate && (m_preserveDescendants || m_preserveGeneratorDescendants || m_preserveAncestors)) {
171 ATH_MSG_FATAL(
"You are asking to preserve only parents/children/siblings of retained states, and also more distant relatives. Please check your job options.");
172 return StatusCode::FAILURE;
178 ATH_CHECK( m_particlesKey.initialize (m_streamName) );
179 ATH_CHECK( m_verticesKey.initialize (m_streamName) );
180 return StatusCode::SUCCESS;
186 ATH_MSG_INFO(
"Total of " << m_totpart <<
" truth particles");
187 ATH_MSG_INFO(
"Removed "<< m_removedpart <<
" truth particles");
188 return StatusCode::SUCCESS;
194 const EventContext& ctx = Gaudi::Hive::currentContext();
198 (m_particlesKey, ctx);
199 m_totpart += importedTruthParticles->size();
202 (m_verticesKey, ctx);
205 if (evtStore()->
retrieve(importedTruthEvents,m_eventsKey).isFailure()) {
206 ATH_MSG_ERROR(
"No TruthEventContainer with name " << m_eventsKey <<
" found in StoreGate!");
207 return StatusCode::FAILURE;
217 std::vector<bool> particleMask, vertexMask;
218 int nTruthParticles = importedTruthParticles->size();
219 int nTruthVertices = importedTruthVertices->size();
220 particleMask.assign(nTruthParticles,
true);
221 vertexMask.assign(nTruthVertices,
false);
224 for (
int particleCounter = 0; particleCounter < nTruthParticles; ++particleCounter) {
227 if ( (particleCounter > m_writeFirstN) && (!
isAccepted(particle)) ){
229 particleMask[particleCounter] =
false;
241 std::unordered_set<int> encounteredUniqueIDs;
242 if (m_preserveDescendants || m_preserveGeneratorDescendants || m_preserveAncestors) {
243 for (
int i=0;
i<nTruthParticles; ++
i) {
244 bool toKeep = particleMask[
i];
245 if (!toKeep)
continue;
247 encounteredUniqueIDs.clear();
248 if (m_preserveDescendants) decayHelper.
descendants(particle,particleMask,vertexMask,encounteredUniqueIDs,
true);
249 encounteredUniqueIDs.clear();
250 if (m_preserveGeneratorDescendants) decayHelper.
descendants(particle,particleMask,vertexMask,encounteredUniqueIDs,
false);
251 encounteredUniqueIDs.clear();
252 if (m_preserveAncestors) decayHelper.
ancestors(particle,particleMask,vertexMask,encounteredUniqueIDs);
253 encounteredUniqueIDs.clear();
258 if (m_preserveImmediate) {
261 std::vector<bool> particleMaskCopy = particleMask;
262 for (
int i=0;
i<nTruthParticles; ++
i) {
263 bool toKeep = particleMask[
i];
264 if (!toKeep)
continue;
269 particleMask = particleMaskCopy;
273 importedTruthParticles.
keep (particleMask);
274 importedTruthVertices.
keep (vertexMask);
276 return StatusCode::SUCCESS;
284 int pdg_id = std::abs(
p->pdgId());
288 if(!m_longLivedPdgIds.empty() && parentIsLongLived(
p)) ok =
true;
292 if (!(
MC::isPhoton(pdg_id) && m_geantPhotonPtThresh >= 0 &&
p->pt() > m_geantPhotonPtThresh) )
297 if (m_writeEverything) ok =
true;
304 (m_partonPtThresh<0 || p->
pt()>m_partonPtThresh) )
321 int motherPDGID = 999999999;
328 if (mother) motherPDGID = mother->
pdgId();
330 if( m_writePartons && !
MC::isHadron( motherPDGID ) ) ok =
true;
331 if( m_writeHadrons &&
MC::isHadron( motherPDGID ) ) ok =
true;
337 if(isLeptonFromTau(
p)) ok =
true;
339 if(isFsrFromLepton(
p)) ok =
true;
343 std::unordered_set<int> uniqueID_trace;
355 if (m_writettHFHadrons && isttHFHadron(
p)) ok =
true;
358 if(m_writeTopAndDecays &&
MC::isTop(
p)) ok =
true;
368 ok = !(matchHadronIncTau(
p) || matchQuarkIncTau(
p) || isOrphanIncTau(
p));
371 if ((m_writeBSMProducts || m_writeBosonProducts || m_writeTopAndDecays) &&
p->hasProdVtx()){
374 for(
unsigned int itr=0; itr<nIncoming; ++itr) {
376 if (!incomingParticle)
continue;
377 if ((m_writeBSMProducts &&
MC::isBSM(incomingParticle)) ||
378 (m_writeBosonProducts && (
MC::isZ(incomingParticle)||
MC::isW(incomingParticle)||
MC::isHiggs(incomingParticle) || (
MC::isPhoton(incomingParticle) && incomingParticle->
pt()>m_photonPtCut) )) ||
379 (m_writeTopAndDecays &&
MC::isTop(incomingParticle)) ){
388 for (
const auto id : m_pdgIdsToKeep){
389 if (pdg_id==
id) ok=
true;
398 std::vector<int> hadrons = std::vector<int>{111};
399 std::vector<int> taus = std::vector<int>{15,-15};
400 return matchGenParticle(
part, hadrons, taus,
true);
406 std::vector<int> quarks = std::vector<int>{1,6};
407 std::vector<int> taus = std::vector<int>{15,-15};
408 return matchGenParticle(
part, quarks, taus,
true);
416 int pdgId =
part->pdgId();
418 if (!(
part->hasProdVtx()))
return true;
422 unsigned int itrParent = 0;
425 if (prodVtx==
part->decayVtx())
return false;
428 if(nIncoming == 0)
return true;
430 unsigned int nParents = 0;
431 while(nParents==0 && itrParent<nIncoming) {
435 if(incomingParticle->
pdgId() == pdgId) {
436 if(!isOrphanIncTau(incomingParticle)) nParents++;
439 if(!isOrphanIncTau(incomingParticle)) nParents++;
446 return (nParents==0);
450 std::vector<int> &intermediatePdgIds,
bool targetsAreRange)
const {
454 int pdgId =
part->pdgId();
459 std::vector<int>::const_iterator itrPdgId, itrPdgIdEnd;
461 if (!(
part->hasProdVtx()))
return false;
468 if (prodVtx==
part->decayVtx())
return false;
470 unsigned int itrParent = 0;
471 while(!
found && itrParent<nIncoming) {
474 if(!targetsAreRange) {
477 itrPdgId = targetPdgIds.begin();
478 itrPdgIdEnd = targetPdgIds.end();
479 for(;itrPdgId != itrPdgIdEnd; ++itrPdgId) {
480 if(incomingParticle->
pdgId() == (*itrPdgId))
return true;
484 int absPdgId = abs(incomingParticle->
pdgId());
487 if(targetPdgIds.size() == 1) {
488 if(absPdgId >= targetPdgIds.at(0))
return true;
490 else if(targetPdgIds.size() >= 2) {
491 if(absPdgId >= targetPdgIds.at(0) &&
492 absPdgId <= targetPdgIds.at(1))
return true;
497 if(incomingParticle->
pdgId() == pdgId) {
498 found = matchGenParticle(incomingParticle, targetPdgIds, intermediatePdgIds, targetsAreRange);
502 itrPdgId = intermediatePdgIds.begin();
503 itrPdgIdEnd = intermediatePdgIds.end();
504 bool foundIntermediate =
false;
505 while(!foundIntermediate && itrPdgId != itrPdgIdEnd) {
506 if(incomingParticle->
pdgId() == (*itrPdgId)) foundIntermediate =
true;
509 if(foundIntermediate) {
510 found = matchGenParticle(incomingParticle, targetPdgIds, intermediatePdgIds, targetsAreRange);
521 int pdg =
part->pdgId();
526 if(!prod)
return false;
529 if (prod==
part->decayVtx())
return false;
533 for(
unsigned int itr = 0; itr<nIncoming; ++itr){
536 ATH_MSG_DEBUG(
"Particle with pdgId = " << pdg <<
", matched to tau");
554 std::unordered_set<int>& uniqueID_trace)
const {
556 int pdg =
part->pdgId();
559 if(!prod)
return false;
562 if (prod==
part->decayVtx())
return false;
565 std::unordered_set<int>::const_iterator foundVtx = uniqueID_trace.find(
HepMC::uniqueID(prod) );
566 if( foundVtx != uniqueID_trace.end() ) {
567 ATH_MSG_DEBUG(
"Found a loop (a la Sherpa sample). Backing out." );
574 for (
unsigned int pitr = 0; pitr<nIncoming; ++pitr){
576 if (!itrParent)
continue;
582 bool has_fsr =
false;
585 for (
unsigned int citr = 0; citr<nChildren; ++citr) {
587 if (!itrChild)
continue;
594 if (has_fsr)
return false;
595 ATH_MSG_DEBUG(
"Particle with pdgId = " << pdg <<
", matched to tau");
600 if(
isFromTau(itrParent, uniqueID_trace)) {
612 int ttHFClassification=6;
616 ttHFClassification = TopHadronOriginFlagAcc(*
part);
622 if (ttHFClassification < 6 )
630 int pdg =
part->pdgId();
634 if(!prod)
return false;
636 if (prod==
part->decayVtx())
return false;
639 for(
unsigned int pitr=0; pitr<nIncoming; ++pitr){
650 if(isFsrFromLepton(itrParent))
return true;
659 for(
size_t parent_itr = 0; parent_itr <
part->nParents(); parent_itr++){
660 if(!
part->parent(parent_itr))
continue;
662 const int parent_abs_pdgid = abs(
parent->pdgId());
663 if(
std::find(m_longLivedPdgIds.begin(), m_longLivedPdgIds.end(), parent_abs_pdgid) != m_longLivedPdgIds.end() ){