20 #include "GaudiKernel/SystemOfUnits.h"
22 #include "GaudiKernel/ThreadLocalContext.h"
32 const IInterface*
p ) :
38 declareProperty (
"WritePartons",
42 declareProperty (
"WriteHadrons",
46 declareProperty (
"WriteBHadrons",
50 declareProperty (
"WriteCHadrons",
54 declareProperty (
"WriteGeant",
56 "Keep Geant particles?");
58 declareProperty (
"GeantPhotonPtThresh",
60 "Write Geant photons with Pt above this threshold. "
61 "Set to < 0 to not write any.");
63 declareProperty (
"WriteTauHad",
65 "Keep hadronic taus?");
67 declareProperty (
"PartonPtThresh",
69 "Write partons with Pt above this threshold.");
71 declareProperty (
"WriteBSM",
73 "Keep BSM particles?");
75 declareProperty (
"WriteBosons",
79 declareProperty (
"PhotonPtThresh",
81 "Photon pt cut with WriteBosons");
83 declareProperty (
"WriteBSMProducts",
85 "Keep BSM particle decay products?");
87 declareProperty (
"WriteBosonProducts",
89 "Keep boson decay products?");
91 declareProperty (
"WriteTopAndDecays",
93 "Keep top partons and immediate decay products?");
95 declareProperty(
"WriteEverything",
97 "Keep absolutely everything (overrides all other flags)");
99 declareProperty(
"WriteAllLeptons",
101 "Keep absolutely all leptons");
103 declareProperty(
"WriteLeptonsNotFromHadrons",
105 "Keep leptons not from hadron decays");
107 declareProperty(
"WriteAllStable",
109 "Keep all stable particles");
111 declareProperty(
"WriteNotPhysical",
113 "Save also non-physical particles");
115 declareProperty(
"WriteFirstN",
117 "Keep first N particles in record");
119 declareProperty(
"PreserveDescendants",
121 "Preserve descendants of retained particles");
123 declareProperty(
"PreserveGeneratorDescendants",
125 "Preserve descendants of retained particles excluding Geant particles");
127 declareProperty(
"PreserveAncestors",
129 "Preserve ancestors of retained particles");
131 declareProperty (
"PreserveParentsSiblingsChildren",
133 "Preserve the parents, siblings and children of retained particles");
135 declareProperty (
"PreserveHadronizationVertices",
137 "Preserve hadronization vertices with parents/children.");
139 declareProperty (
"WritettHFHadrons",
141 "Keep tt+HF hadrons?");
143 declareProperty (
"PDGIDsToKeep",
145 "List of PDG IDs to always keep");
147 declareProperty (
"LongLivedPDGIDs",
149 "List of PDG IDs of long lived particles so that one can store their children");
159 if (m_preserveDescendants && m_preserveGeneratorDescendants) {
160 ATH_MSG_FATAL(
"You are asking to keep both all descendants, and only those from the event generator. Please check your job options.");
161 return StatusCode::FAILURE;
163 if (m_preserveImmediate && (m_preserveDescendants || m_preserveGeneratorDescendants || m_preserveAncestors)) {
164 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.");
165 return StatusCode::FAILURE;
171 ATH_CHECK( m_particlesKey.initialize (m_streamName) );
172 ATH_CHECK( m_verticesKey.initialize (m_streamName) );
173 return StatusCode::SUCCESS;
179 ATH_MSG_INFO(
"Total of " << m_totpart <<
" truth particles");
180 ATH_MSG_INFO(
"Removed "<< m_removedpart <<
" truth particles");
181 return StatusCode::SUCCESS;
187 const EventContext& ctx = Gaudi::Hive::currentContext();
191 (m_particlesKey, ctx);
192 m_totpart += importedTruthParticles->size();
195 (m_verticesKey, ctx);
204 std::vector<bool> particleMask, vertexMask;
205 int nTruthParticles = importedTruthParticles->size();
206 int nTruthVertices = importedTruthVertices->size();
207 particleMask.assign(nTruthParticles,
true);
208 vertexMask.assign(nTruthVertices,
false);
211 for (
int particleCounter = 0; particleCounter < nTruthParticles; ++particleCounter) {
216 particleMask[particleCounter] =
false;
228 std::unordered_set<int> encounteredUniqueIDs;
229 if (m_preserveDescendants || m_preserveGeneratorDescendants || m_preserveAncestors) {
230 for (
int i=0;
i<nTruthParticles; ++
i) {
231 bool toKeep = particleMask[
i];
232 if (!toKeep)
continue;
234 encounteredUniqueIDs.clear();
235 if (m_preserveDescendants) decayHelper.
descendants(
particle,particleMask,vertexMask,encounteredUniqueIDs,
true);
236 encounteredUniqueIDs.clear();
237 if (m_preserveGeneratorDescendants) decayHelper.
descendants(
particle,particleMask,vertexMask,encounteredUniqueIDs,
false);
238 encounteredUniqueIDs.clear();
239 if (m_preserveAncestors) decayHelper.
ancestors(
particle,particleMask,vertexMask,encounteredUniqueIDs);
240 encounteredUniqueIDs.clear();
245 if (m_preserveImmediate) {
248 std::vector<bool> particleMaskCopy = particleMask;
249 for (
int i=0;
i<nTruthParticles; ++
i) {
250 bool toKeep = particleMask[
i];
251 if (!toKeep)
continue;
256 particleMask = particleMaskCopy;
260 importedTruthParticles.
keep (particleMask);
261 importedTruthVertices.
keep (vertexMask);
263 return StatusCode::SUCCESS;
271 int pdg_id =
p->absPdgId();
275 if(!m_longLivedPdgIds.empty() && parentIsLongLived(
p)) ok =
true;
279 if (!(
MC::isPhoton(pdg_id) && m_geantPhotonPtThresh >= 0 &&
p->pt() > m_geantPhotonPtThresh) )
284 if (m_writeEverything) ok =
true;
291 (m_partonPtThresh<0 || p->
pt()>m_partonPtThresh) )
308 int motherPDGID = 999999999;
315 if (mother) motherPDGID = mother->
pdgId();
317 if( m_writePartons && !
MC::isHadron( motherPDGID ) ) ok =
true;
318 if( m_writeHadrons &&
MC::isHadron( motherPDGID ) ) ok =
true;
324 if(isLeptonFromTau(
p)) ok =
true;
326 if(isFsrFromLepton(
p)) ok =
true;
330 std::unordered_set<int> uniqueID_trace;
342 if (m_writettHFHadrons && isttHFHadron(
p)) ok =
true;
345 if(m_writeTopAndDecays &&
MC::isTop(
p)) ok =
true;
355 ok = !(matchHadronIncTau(
p) || matchQuarkIncTau(
p) || isOrphanIncTau(
p));
358 if ((m_writeBSMProducts || m_writeBosonProducts || m_writeTopAndDecays) &&
p->hasProdVtx()){
361 for(
unsigned int itr=0; itr<nIncoming; ++itr) {
363 if (!incomingParticle)
continue;
364 if ((m_writeBSMProducts &&
MC::isBSM(incomingParticle)) ||
365 (m_writeBosonProducts && (
MC::isZ(incomingParticle)||
MC::isW(incomingParticle)||
MC::isHiggs(incomingParticle) || (
MC::isPhoton(incomingParticle) && incomingParticle->
pt()>m_photonPtCut) )) ||
366 (m_writeTopAndDecays &&
MC::isTop(incomingParticle)) ){
375 for (
const auto id : m_pdgIdsToKeep){
376 if (pdg_id==
id) ok=
true;
385 std::vector<int>
hadrons = std::vector<int>{111};
386 std::vector<int> taus = std::vector<int>{15,-15};
387 return matchGenParticle(
part,
hadrons, taus,
true);
393 std::vector<int> quarks = std::vector<int>{1,6};
394 std::vector<int> taus = std::vector<int>{15,-15};
395 return matchGenParticle(
part, quarks, taus,
true);
405 if (!
part->hasProdVtx())
return true;
409 unsigned int itrParent = 0;
412 if (prodVtx==
part->decayVtx())
return false;
415 if(nIncoming == 0)
return true;
417 unsigned int nParents = 0;
418 while(nParents==0 && itrParent<nIncoming) {
423 if(!isOrphanIncTau(incomingParticle)) nParents++;
426 if(!isOrphanIncTau(incomingParticle)) nParents++;
433 return (nParents==0);
437 std::vector<int> &intermediatePdgIds,
bool targetsAreRange)
const {
446 std::vector<int>::const_iterator itrPdgId, itrPdgIdEnd;
448 if (!
part->hasProdVtx())
return false;
455 if (prodVtx==
part->decayVtx())
return false;
457 unsigned int itrParent = 0;
458 while(!
found && itrParent<nIncoming) {
461 if(!targetsAreRange) {
464 itrPdgId = targetPdgIds.begin();
465 itrPdgIdEnd = targetPdgIds.end();
466 for(;itrPdgId != itrPdgIdEnd; ++itrPdgId) {
467 if(incomingParticle->
pdgId() == (*itrPdgId))
return true;
471 int absPdgId = incomingParticle->
absPdgId();
474 if(targetPdgIds.size() == 1) {
475 if(absPdgId >= targetPdgIds.at(0))
return true;
477 else if(targetPdgIds.size() >= 2) {
478 if(absPdgId >= targetPdgIds.at(0) &&
479 absPdgId <= targetPdgIds.at(1))
return true;
485 found = matchGenParticle(incomingParticle, targetPdgIds, intermediatePdgIds, targetsAreRange);
489 itrPdgId = intermediatePdgIds.begin();
490 itrPdgIdEnd = intermediatePdgIds.end();
491 bool foundIntermediate =
false;
492 while(!foundIntermediate && itrPdgId != itrPdgIdEnd) {
493 if(incomingParticle->
pdgId() == (*itrPdgId)) foundIntermediate =
true;
496 if(foundIntermediate) {
497 found = matchGenParticle(incomingParticle, targetPdgIds, intermediatePdgIds, targetsAreRange);
508 int pdg =
part->pdgId();
513 if(!prod)
return false;
516 if (prod==
part->decayVtx())
return false;
520 for(
unsigned int itr = 0; itr<nIncoming; ++itr){
523 ATH_MSG_DEBUG(
"Particle with pdgId = " << pdg <<
", matched to tau");
541 std::unordered_set<int>& uniqueID_trace)
const {
543 int pdg =
part->pdgId();
546 if(!prod)
return false;
549 if (prod==
part->decayVtx())
return false;
552 std::unordered_set<int>::const_iterator foundVtx = uniqueID_trace.find(
HepMC::uniqueID(prod) );
553 if( foundVtx != uniqueID_trace.end() ) {
554 ATH_MSG_DEBUG(
"Found a loop (a la Sherpa sample). Backing out." );
561 for (
unsigned int pitr = 0; pitr<nIncoming; ++pitr){
563 if (!itrParent)
continue;
569 bool has_fsr =
false;
572 for (
unsigned int citr = 0; citr<nChildren; ++citr) {
574 if (!itrChild)
continue;
581 if (has_fsr)
return false;
582 ATH_MSG_DEBUG(
"Particle with pdgId = " << pdg <<
", matched to tau");
587 if(
isFromTau(itrParent, uniqueID_trace)) {
599 int ttHFClassification=6;
603 ttHFClassification = TopHadronOriginFlagAcc(*
part);
609 if (ttHFClassification < 6 )
617 int pdg =
part->pdgId();
621 if(!prod)
return false;
623 if (prod==
part->decayVtx())
return false;
626 for(
unsigned int pitr=0; pitr<nIncoming; ++pitr){
637 if(isFsrFromLepton(itrParent))
return true;
646 for(
size_t parent_itr = 0; parent_itr <
part->nParents(); parent_itr++){
647 if(!
part->parent(parent_itr))
continue;
649 const int parent_abs_pdgid = abs(
parent->pdgId());
650 if(
std::find(m_longLivedPdgIds.begin(), m_longLivedPdgIds.end(), parent_abs_pdgid) != m_longLivedPdgIds.end() ){