HepMCCandAlgos/src/GenParticlesHelper.cc

0001 #include "PhysicsTools/HepMCCandAlgos/interface/GenParticlesHelper.h"
0002 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0003
0004 using namespace reco;
0005
0006 namespace GenParticlesHelper {
0007
0008   void findParticles(const reco::GenParticleCollection& sourceParticles,
0009                      reco::GenParticleRefVector& particleRefs,
0010                      int pdgId,
0011                      int status) {
0012     unsigned index = 0;
0013     for (IG ig = sourceParticles.begin(); ig != sourceParticles.end(); ++ig, ++index) {
0014       const GenParticle& gen = *ig;
0015
0016       // status has been specified, and this one does not have the correct
0017       // status
0018       if (status && gen.status() != status)
0019         continue;
0020
0021       if (std::abs(gen.pdgId()) == pdgId) {
0022         GenParticleRef genref(&sourceParticles, index);
0023         particleRefs.push_back(genref);
0024       }
0025     }
0026   }
0027
0028   void findDescendents(const reco::GenParticleRef& base,
0029                        reco::GenParticleRefVector& descendents,
0030                        int status,
0031                        int pdgId) {
0032     const GenParticleRefVector& daughterRefs = base->daughterRefVector();
0033
0034     for (IGR idr = daughterRefs.begin(); idr != daughterRefs.end(); ++idr) {
0035       if ((*idr)->status() == status && (!pdgId || std::abs((*idr)->pdgId()) == pdgId)) {
0036         // cout<<"adding "<<(*idr)<<endl;
0037         descendents.push_back(*idr);
0038       } else
0039         findDescendents(*idr, descendents, status, pdgId);
0040     }
0041   }
0042
0043   void findSisters(const reco::GenParticleRef& baseSister, GenParticleRefVector& sisterRefs) {
0044     assert(baseSister->numberOfMothers() > 0);
0045
0046     // get first mother
0047     const GenParticleRefVector& mothers = baseSister->motherRefVector();
0048
0049     // get sisters
0050     const GenParticleRefVector allRefs = mothers[0]->daughterRefVector();
0051
0052     typedef GenParticleRefVector::const_iterator IT;
0053     for (IT id = allRefs.begin(); id != allRefs.end(); ++id) {
0054       if (*id == baseSister) {
0055         continue;  // this is myself
0056       } else
0057         sisterRefs.push_back(*id);
0058     }
0059   }
0060
0061   bool isDirect(const reco::GenParticleRef& particle) {
0062     assert((particle->status() != 0) && (particle->status() < 4));
0063     if (particle->status() == 3)
0064       return true;
0065     else {
0066       assert(particle->numberOfMothers() > 0);
0067
0068       // get first mother
0069       const GenParticleRefVector& mothers = particle->motherRefVector();
0070       if (mothers[0]->status() == 3)
0071         return true;
0072       else
0073         return false;
0074     }
0075   }
0076
0077   bool hasAncestor(const reco::GenParticle* particle, int pdgId, int status) {
0078     if (particle->pdgId() == pdgId && particle->status() == status)
0079       return true;
0080
0081     const GenParticleRefVector& mothers = particle->motherRefVector();
0082
0083     for (IGR im = mothers.begin(); im != mothers.end(); ++im) {
0084       const GenParticle& part = **im;
0085       if (hasAncestor(&part, pdgId, status))
0086         return true;
0087     }
0088
0089     return false;
0090   }
0091
0092 }  // namespace GenParticlesHelper
0093
0094 namespace edm {
0095   std::ostream& operator<<(std::ostream& out, const reco::GenParticleRef& genRef) {
0096     if (!out)
0097       return out;
0098
0099     out << "key:" << genRef.key() << " pt:" << genRef->pt();
0100
0101     return out;
0102   }
0103
0104 }  // namespace edm