rivet is hosted by Hepforge, IPPP Durham
Jet.cc
Go to the documentation of this file.
00001 #include "Rivet/Jet.hh"
00002 #include "Rivet/Tools/Cuts.hh"
00003 #include "Rivet/Tools/ParticleName.hh"
00004 #include "Rivet/Tools/Logging.hh"
00005 #include "Rivet/Tools/ParticleIdUtils.hh"
00006 
00007 namespace Rivet {
00008 
00009 
00010   Jet& Jet::clear() {
00011     _momentum = FourMomentum();
00012     _pseudojet.reset(0,0,0,0);
00013     _particles.clear();
00014     return *this;
00015   }
00016 
00017 
00018   Jet& Jet::setState(const FourMomentum& mom, const Particles& particles, const Particles& tags) {
00019     clear();
00020     _momentum = mom;
00021     _pseudojet = fastjet::PseudoJet(mom.px(), mom.py(), mom.pz(), mom.E());
00022     _particles = particles;
00023     _tags = tags;
00024     return *this;
00025   }
00026 
00027 
00028   Jet& Jet::setState(const fastjet::PseudoJet& pj, const Particles& particles, const Particles& tags) {
00029     clear();
00030     _pseudojet = pj;
00031     _momentum = FourMomentum(pj.e(), pj.px(), pj.py(), pj.pz());
00032     _particles = particles;
00033     _tags = tags;
00034     // if (_particles.empty()) {
00035     //   foreach (const fastjet::PseudoJet pjc, _pseudojet.constituents()) {
00036     //     // If there is no attached user info, we can't create a meaningful particle, so skip
00037     //     if (!pjc.has_user_info<RivetFJInfo>()) continue;
00038     //     const RivetFJInfo& fjinfo = pjc.user_info<RivetFJInfo>();
00039     //     // Don't add ghosts to the particles list
00040     //     if (fjinfo.isGhost) continue;
00041     //     // Otherwise construct a Particle from the PseudoJet, preferably from an associated GenParticle
00042     //     ?if (fjinfo.genParticle != NULL) {
00043     //       _particles.push_back(Particle(fjinfo.genParticle));
00044     //     } else {
00045     //       if (fjinfo.pid == 0) continue; // skip if there is a null PID entry in the FJ info
00046     //       const FourMomentum pjcmom(pjc.e(), pjc.px(), pjc.py(), pjc.pz());
00047     //       _particles.push_back(Particle(fjinfo.pid, pjcmom));
00048     //     }
00049     //   }
00050     // }
00051     return *this;
00052   }
00053 
00054 
00055   Jet& Jet::setParticles(const Particles& particles) {
00056     _particles = particles;
00057     return *this;
00058   }
00059 
00060 
00061   bool Jet::containsParticle(const Particle& particle) const {
00062     const int barcode = particle.genParticle()->barcode();
00063     foreach (const Particle& p, particles()) {
00064       if (p.genParticle()->barcode() == barcode) return true;
00065     }
00066     return false;
00067   }
00068 
00069 
00070   bool Jet::containsParticleId(PdgId pid) const {
00071     foreach (const Particle& p, particles()) {
00072       if (p.pid() == pid) return true;
00073     }
00074     return false;
00075   }
00076 
00077 
00078   bool Jet::containsParticleId(const vector<PdgId>& pids) const {
00079     foreach (const Particle& p, particles()) {
00080       foreach (PdgId pid, pids) {
00081         if (p.pid() == pid) return true;
00082       }
00083     }
00084     return false;
00085   }
00086 
00087 
00088   /// @todo Jet::containsMatch(Matcher m) { ... if m(pid) return true; ... }
00089 
00090 
00091   Jet& Jet::transformBy(const LorentzTransform& lt) {
00092     _momentum = lt.transform(_momentum);
00093     for (Particle& p : _particles) p.transformBy(lt);
00094     for (Particle& t : _tags) t.transformBy(lt);
00095     _pseudojet.reset(_momentum.px(), _momentum.py(), _momentum.pz(), _momentum.E()); //< lose ClusterSeq etc.
00096     return *this;
00097   }
00098 
00099 
00100   double Jet::neutralEnergy() const {
00101     double e_neutral = 0.0;
00102     foreach (const Particle& p, particles()) {
00103       const PdgId pid = p.pid();
00104       if (PID::threeCharge(pid) == 0) {
00105         e_neutral += p.E();
00106       }
00107     }
00108     return e_neutral;
00109   }
00110 
00111 
00112   double Jet::hadronicEnergy() const {
00113     double e_hadr = 0.0;
00114     foreach (const Particle& p, particles()) {
00115       const PdgId pid = p.pid();
00116       if (PID::isHadron(pid)) {
00117         e_hadr += p.E();
00118       }
00119     }
00120     return e_hadr;
00121   }
00122 
00123 
00124   bool Jet::containsCharm(bool include_decay_products) const {
00125     foreach (const Particle& p, particles()) {
00126       const PdgId pid = p.pid();
00127       if (abs(pid) == PID::CQUARK) return true;
00128       if (PID::isHadron(pid) && PID::hasCharm(pid)) return true;
00129       if (include_decay_products) {
00130         const HepMC::GenVertex* gv = p.genParticle()->production_vertex();
00131         if (gv) {
00132           foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) {
00133             const PdgId pid2 = pi->pdg_id();
00134             if (PID::isHadron(pid2) && PID::hasCharm(pid2)) return true;
00135           }
00136         }
00137       }
00138     }
00139     return false;
00140   }
00141 
00142 
00143   bool Jet::containsBottom(bool include_decay_products) const {
00144     foreach (const Particle& p, particles()) {
00145       const PdgId pid = p.pid();
00146       if (abs(pid) == PID::BQUARK) return true;
00147       if (PID::isHadron(pid) && PID::hasBottom(pid)) return true;
00148       if (include_decay_products) {
00149         const HepMC::GenVertex* gv = p.genParticle()->production_vertex();
00150         if (gv) {
00151           foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) {
00152             const PdgId pid2 = pi->pdg_id();
00153             if (PID::isHadron(pid2) && PID::hasBottom(pid2)) return true;
00154           }
00155         }
00156       }
00157     }
00158     return false;
00159   }
00160 
00161 
00162   Particles Jet::tags(const Cut& c) const {
00163     return filter_select(tags(), c);
00164   }
00165 
00166   Particles Jet::bTags(const Cut& c) const {
00167     Particles rtn;
00168     for (const Particle& tp : tags()) {
00169       if (hasBottom(tp) && c->accept(tp)) rtn.push_back(tp);
00170     }
00171     return rtn;
00172   }
00173 
00174   Particles Jet::cTags(const Cut& c) const {
00175     Particles rtn;
00176     for (const Particle& tp : tags()) {
00177       /// @todo Is making b and c tags exclusive the right thing to do?
00178       if (hasCharm(tp) && !hasBottom(tp) && c->accept(tp)) rtn.push_back(tp);
00179     }
00180     return rtn;
00181   }
00182 
00183   Particles Jet::tauTags(const Cut& c) const {
00184     Particles rtn;
00185     for (const Particle& tp : tags()) {
00186       if (isTau(tp) && c->accept(tp)) rtn.push_back(tp);
00187     }
00188     return rtn;
00189   }
00190 
00191 
00192 }