20 #include "GaudiKernel/SystemOfUnits.h"
22 #include "GaudiKernel/ThreadLocalContext.h"
43 const std::string&
t,
const std::string&
n,
const IInterface*
p ) :
54 "truth jet minumum pt");
57 "truth jet maximum abs(eta)");
60 "jet constituent minimum pt");
63 "jet constituent photon minimum pt");
66 "list of abs(pdgID) to keep");
68 declareProperty(
"IsolRadius",
m_isolR,
69 "isolation radius for leptons and photons");
72 "isolation particle minimum pt");
74 declareProperty(
"MaxCount",
76 "maximum number of events to print");
95 ATH_CHECK( m_hardParticleKey.initialize() );
97 ATH_CHECK( m_truthParticleName.initialize (m_streamName) );
98 ATH_CHECK( m_truthVertexName.initialize (m_streamName) );
104 if( m_jetPhotonPtCut < m_jetConstPtCut ) m_jetPhotonPtCut =m_jetConstPtCut;
106 return StatusCode::SUCCESS;
113 if( m_errCount > 0 ){
114 ATH_MSG_WARNING(
"TruthHard/TruthEvent pdgId mismatches " <<m_errCount);
116 ATH_MSG_INFO(
"No TruthHard/TruthEvent pdgId mismatches");
119 return StatusCode::SUCCESS;
129 const EventContext& ctx = Gaudi::Hive::currentContext();
132 bool doPrint = m_evtCount < m_maxCount;
136 (m_truthParticleName, ctx);
138 (m_truthVertexName, ctx);
142 ATH_MSG_ERROR(
"Cannot retrieve TruthParticleContainer " << m_hardParticleKey);
143 return StatusCode::FAILURE;
151 std::vector<const xAOD::TruthParticle*> hardPart;
152 std::vector<TLorentzVector> pLepGam;
154 for (
const auto* pItr : *inHardParts) {
156 hardPart.push_back(pItr);
158 int ida = pItr->absPdgId();
160 pLepGam.push_back( pItr->p4() );
167 long long int evtNum{};
172 return StatusCode::FAILURE;
174 evtNum =
evt->eventNumber();
176 printxAODTruth(evtNum, inTruthParts.cptr());
183 std::vector<bool> partMask;
184 std::vector<bool> vertMask;
185 unsigned int inPartNum = inTruthParts->size();
186 unsigned int inVertNum = inTruthVerts->size();
187 partMask.assign(inPartNum,
false);
188 vertMask.assign(inVertNum,
false);
190 ATH_MSG_DEBUG(
"Initial particles/vertices = " <<inPartNum <<
" " <<inVertNum);
199 std::vector<const xAOD::TruthParticle*> kids;
201 for (
const auto* pItr : *inTruthParts) {
205 for(
unsigned int i=0;
i<hardPart.size(); ++
i){
208 if( pItr->pdgId() != (hardPart[
i])->pdgId() ){
212 <<
" " <<(hardPart[
i])->pdgId());
220 partMask[ pItr->index() ] =
true;
224 int ida = pItr->absPdgId();
226 for(
unsigned int i=0;
i<m_keepIds.size(); ++
i){
227 if( ida == m_keepIds[
i] ){
234 partMask[pItr->index()] =
true;
235 int nkids = getDescendants( pItr, kids );
236 for(
int i=0;
i<nkids; ++
i){
239 partMask[ (kids[
i])->
index() ] =
true;
241 if(
v ) vertMask[
v->index() ] =
true;
246 if( !pLepGam.empty() && pItr->pt() > m_isolPtCut ){
247 TLorentzVector pp4 = pItr->p4();
248 for(
unsigned int lep=0; lep<pLepGam.size(); ++lep){
249 double r = pp4.DeltaR( pLepGam[lep] );
252 partMask[ pItr->index() ] =
true;
263 if( !m_truthJetsKey.empty() ){
267 ATH_MSG_ERROR(
"Cannot retrieve JetContainer " << m_truthJetsKey);
268 return StatusCode::FAILURE;
271 std::vector<int> uidJetConst;
273 for(
const auto* ajet : *inJets){
274 if( ajet->pt() < m_jetPtCut )
continue;
275 if( std::abs(ajet->eta()) > m_jetEtaCut )
continue;
280 for( ; aItr != aItrE; ++aItr){
298 for (
const auto* pItr : *inTruthParts) {
301 for(
unsigned int i=0;
i<uidJetConst.size(); ++
i){
302 if(
uid == uidJetConst[
i] ){
309 partMask[ pItr->index() ] =
true;
317 inTruthParts.keep (partMask);
318 inTruthVerts.
keep (vertMask);
322 for(
unsigned int i=0;
i<partMask.size(); ++
i){
323 if( partMask[
i] ) ++outPartNum;
326 for(
unsigned int i=0;
i<vertMask.size(); ++
i){
327 if( vertMask[
i] ) ++outVertNum;
330 ATH_MSG_DEBUG(
"Final particles/vertices = " <<outPartNum <<
" " <<outVertNum);
333 std::cout <<
"======================================================================================" <<std::endl;
334 std::cout <<
"HardTruthThinning complete for event " <<evtNum <<std::endl;
335 std::cout <<
"Saved " <<outPartNum <<
" particles" <<std::endl;
336 std::cout <<
"Particle unique IDs = ";
337 for(
unsigned int i=0;
i<partMask.size(); ++
i){
340 std::cout <<std::endl;
342 std::cout <<
"Saved " <<outVertNum <<
" vertices" <<std::endl;
343 std::cout <<
"Vertex unique IDs = ";
344 for(
unsigned int i=0;
i<vertMask.size(); ++
i){
347 std::cout <<std::endl;
348 std::cout <<
"======================================================================================" <<std::endl;
351 return StatusCode::SUCCESS;
367 std::vector<const xAOD::TruthParticle*>& descendants ) {
369 if( ! (
p->hasDecayVtx()) )
return 0;
371 if( !dvtx )
return 0;
374 const std::vector< ElementLink< xAOD::TruthParticleContainer > >& outPart =
376 for(
unsigned int k=0;
k<outPart.size(); ++
k){
377 if( ! (outPart[
k]).
isValid() )
continue;
379 descendants.push_back(kid);
383 int nstop = descendants.size();
385 while( nstop > nstart ){
386 for(
int i=nstart;
i<nstop; ++
i){
393 const std::vector< ElementLink< xAOD::TruthParticleContainer > >&
395 for(
unsigned int k=0;
k<outPart2.size(); ++
k){
396 if( ! (outPart2[
k]).isValid() )
continue;
400 for(
unsigned int kk=0;
kk<descendants.size(); ++
kk){
401 if(kpp==descendants[
kk]) isIn =
true;
403 if( !isIn ) descendants.push_back(kpp);
407 nstop = descendants.size();
420 std::vector<int> uidPars;
421 std::vector<int> uidKids;
423 std::cout <<
"======================================================================================\n" ;
424 std::cout <<
"xAODTruth Event " <<evnum <<
"\n";
425 std::cout <<
" Unique ID PDG Id Status px(GeV) py(GeV) pz(GeV) E(GeV) Parent: Decay\n" ;
426 std::cout <<
" -----------------------------------------------------------------------------------\n" ;
428 for (
const auto* tpItr : *truths) {
431 int id = tpItr->pdgId();
432 int stat = tpItr->status();
433 float px = tpItr->px()/
GeV;
434 float py = tpItr->py()/
GeV;
435 float pz = tpItr->pz()/
GeV;
436 float e = tpItr->e()/
GeV;
440 if( tpItr->hasProdVtx() ){
443 const std::vector< ElementLink< xAOD::TruthParticleContainer > >&
pars =
445 for(
unsigned int k=0;
k<
pars.size(); ++
k){
452 if( tpItr->hasDecayVtx() ){
455 const std::vector< ElementLink< xAOD::TruthParticleContainer > >& kids =
457 for(
unsigned int k=0;
k<kids.size(); ++
k){
458 if( ! (kids[
k]).isValid() )
continue;
465 std::cout <<
std::format(
"{:>10}{:>12}{:>8}{:>10.2f}{:>10.2f}{:>10.2f}{:>10.2f} P: ",
467 for(
unsigned int k=0;
k<uidPars.size(); ++
k){
468 std::cout <<uidPars[
k] <<
" ";
471 for(
unsigned int k=0;
k<uidKids.size(); ++
k){
472 std::cout <<uidKids[
k] <<
" ";
476 std::cout <<
"======================================================================================" <<std::endl;