|
[Rivet-svn] r1643 - in trunk: . include/Rivet include/Rivet/Projections src src/Analyses src/Projectionsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgWed Jun 24 17:27:45 BST 2009
Author: buckley Date: Wed Jun 24 17:23:36 2009 New Revision: 1643 Log: Fixing Jet::containsBottom() and Jet::containsCharm() to dig into the event record a bit, using the new HepMC accessor functions to make HepMC looping a little more pleasant than being repeatedly poked in the eye. Modified: trunk/ChangeLog trunk/include/Rivet/Jet.hh trunk/include/Rivet/Projections/FastJets.hh trunk/include/Rivet/Projections/TotalVisibleMomentum.hh trunk/include/Rivet/Rivet.hh trunk/src/Analyses/ExampleAnalysis.cc trunk/src/Jet.cc trunk/src/Projections/FastJets.cc trunk/src/Projections/FinalState.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Wed Jun 24 17:21:37 2009 (r1642) +++ trunk/ChangeLog Wed Jun 24 17:23:36 2009 (r1643) @@ -1,5 +1,11 @@ 2009-06-24 Andy Buckley <andy at insectnation.org> + * Fixing Jet's containsBottom and containsCharm methods, since B + hadrons are not necessarily to be found in the final + state. Discovered at the same time that HepMC::GenParticle defines + a massively unhelpful copy constructor that actually loses the + tree information; it would be better to hide it entirely! + * Adding RivetHepMC.hh, which defines container-type accessors to HepMC particles and vertices, making it possible to use Boost foreach and hence avoiding the usual huge boilerplate for-loops. Modified: trunk/include/Rivet/Jet.hh ============================================================================== --- trunk/include/Rivet/Jet.hh Wed Jun 24 17:21:37 2009 (r1642) +++ trunk/include/Rivet/Jet.hh Wed Jun 24 17:23:36 2009 (r1643) @@ -81,6 +81,9 @@ /// Check whether this jet contains a certain particle type. bool containsParticleId(PdgId pid) const; + /// Check whether this jet contains at least one of certain particle types. + bool containsParticleId(vector<PdgId> pids) const; + /// Check whether this jet contains a charm-flavoured hadron. bool containsCharm() const; Modified: trunk/include/Rivet/Projections/FastJets.hh ============================================================================== --- trunk/include/Rivet/Projections/FastJets.hh Wed Jun 24 17:21:37 2009 (r1642) +++ trunk/include/Rivet/Projections/FastJets.hh Wed Jun 24 17:23:36 2009 (r1643) @@ -70,21 +70,10 @@ public: /// Reset the projection. Jet def etc are unchanged. - void reset() { - _yscales.clear(); - _particles.clear(); - /// @todo Clear _cseq - //_cseq = fastjet::ClusterSequence(); - } + void reset(); /// Number of jets above the \f$ p_\perp \f$ cut. - size_t numJets(double ptmin = 0.0) const { - if (_cseq.get() != 0) { - return _cseq->inclusive_jets(ptmin).size(); - } else { - return 0; - } - } + size_t numJets(double ptmin = 0.0) const; /// Number of jets. size_t size() const { @@ -92,33 +81,19 @@ } /// Get the jets (unordered). - Jets jets(double ptmin = 0.0) const { - return _pseudojetsToJets(pseudoJets(ptmin)); - } + Jets jets(double ptmin = 0.0) const; /// Get the jets, ordered by \f$ p_T \f$. - Jets jetsByPt(double ptmin = 0.0) const { - return _pseudojetsToJets(pseudoJetsByPt(ptmin)); - } + Jets jetsByPt(double ptmin = 0.0) const; /// Get the jets, ordered by \f$ E \f$. - Jets jetsByE(double ptmin = 0.0) const { - return _pseudojetsToJets(pseudoJetsByE(ptmin)); - } + Jets jetsByE(double ptmin = 0.0) const; /// Get the jets, ordered by rapidity. - Jets jetsByRapidity(double ptmin = 0.0) const { - return _pseudojetsToJets(pseudoJetsByRapidity(ptmin)); - } + Jets jetsByRapidity(double ptmin = 0.0) const; /// Get the pseudo jets (unordered). - PseudoJets pseudoJets(double ptmin = 0.0) const { - if (_cseq.get() != 0) { - return _cseq->inclusive_jets(ptmin); - } else { - return PseudoJets(); - } - } + PseudoJets pseudoJets(double ptmin = 0.0) const; /// Get the pseudo jets, ordered by \f$ p_T \f$. PseudoJets pseudoJetsByPt(double ptmin = 0.0) const { Modified: trunk/include/Rivet/Projections/TotalVisibleMomentum.hh ============================================================================== --- trunk/include/Rivet/Projections/TotalVisibleMomentum.hh Wed Jun 24 17:21:37 2009 (r1642) +++ trunk/include/Rivet/Projections/TotalVisibleMomentum.hh Wed Jun 24 17:23:36 2009 (r1643) @@ -17,7 +17,7 @@ public: - /// Constructor. + /// Constructor. Make sure you supply an appropriately vetoed FS! TotalVisibleMomentum(const FinalState& fsp) { setName("TotalVisibleMomentum"); Modified: trunk/include/Rivet/Rivet.hh ============================================================================== --- trunk/include/Rivet/Rivet.hh Wed Jun 24 17:21:37 2009 (r1642) +++ trunk/include/Rivet/Rivet.hh Wed Jun 24 17:23:36 2009 (r1643) @@ -20,12 +20,6 @@ #include <cassert> #include <fstream> -#include "HepMC/GenEvent.h" -#include "HepMC/GenVertex.h" -#include "HepMC/GenParticle.h" - - -/// This is the main namespace in which all Rivet classes are defined. namespace Rivet { // Convenient imports of common STL classes and functions. @@ -52,10 +46,6 @@ using std::setw; using std::endl; - using HepMC::GenEvent; - using HepMC::GenParticle; - using HepMC::GenVertex; - /// A sensible default maximum value of rapidity for Rivet analyses to use. static const double MaxRapidity = 100000.0; @@ -71,6 +61,9 @@ // Pull some Boost defns into the Rivet namespace #include "Rivet/RivetBoost.hh" +// HepMC headers and helper functions +#include "Rivet/RivetHepMC.hh" + // Now import some Rivet classes #include "Rivet/Exceptions.hh" #include "Rivet/Math/MathUtils.hh" Modified: trunk/src/Analyses/ExampleAnalysis.cc ============================================================================== --- trunk/src/Analyses/ExampleAnalysis.cc Wed Jun 24 17:21:37 2009 (r1642) +++ trunk/src/Analyses/ExampleAnalysis.cc Wed Jun 24 17:23:36 2009 (r1643) @@ -4,6 +4,7 @@ #include "Rivet/Analyses/ExampleAnalysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" +#include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/Multiplicity.hh" #include "Rivet/Projections/Thrust.hh" #include "Rivet/Projections/Sphericity.hh" @@ -17,6 +18,7 @@ const ChargedFinalState cfs(-4, 4, 2*GeV); addProjection(cnfs, "FS"); addProjection(cfs, "CFS"); + addProjection(FastJets(cnfs, FastJets::KT, 0.7), "Jets"); addProjection(Multiplicity(cfs), "CMult"); addProjection(Multiplicity(cnfs), "CNMult"); addProjection(Thrust(cfs), "Thrust"); @@ -60,10 +62,10 @@ // Analyse and print some info const Multiplicity& cm = applyProjection<Multiplicity>(event, "CMult"); const Multiplicity& cnm = applyProjection<Multiplicity>(event, "CNMult"); - getLog() << Log::DEBUG << "Total multiplicity = " << cnm.totalMultiplicity() << endl; - getLog() << Log::DEBUG << "Total charged multiplicity = " << cm.totalMultiplicity() << endl; - getLog() << Log::DEBUG << "Hadron multiplicity = " << cnm.hadronMultiplicity() << endl; - getLog() << Log::DEBUG << "Hadron charged multiplicity = " << cm.hadronMultiplicity() << endl; + getLog() << Log::DEBUG << "Total multiplicity = " << cnm.totalMultiplicity() << endl; + getLog() << Log::DEBUG << "Total charged multiplicity = " << cm.totalMultiplicity() << endl; + getLog() << Log::DEBUG << "Hadron multiplicity = " << cnm.hadronMultiplicity() << endl; + getLog() << Log::DEBUG << "Hadron charged multiplicity = " << cm.hadronMultiplicity() << endl; const Thrust& t = applyProjection<Thrust>(event, "Thrust"); getLog() << Log::DEBUG << "Thrust = " << t.thrust() << endl; @@ -72,6 +74,13 @@ getLog() << Log::DEBUG << "Sphericity = " << s.sphericity() << endl; getLog() << Log::DEBUG << "Aplanarity = " << s.aplanarity() << endl; + size_t num_b_jets = 0; + const Jets jets = applyProjection<FastJets>(event, "Jets").jets(); + foreach (const Jet& j, jets) { + if (j.containsBottom()) ++num_b_jets; + } + getLog() << Log::DEBUG << "#B-jets = " << num_b_jets << endl; + // Fill histograms const double weight = event.weight(); _histTot->fill(cnm.totalMultiplicity(), weight); Modified: trunk/src/Jet.cc ============================================================================== --- trunk/src/Jet.cc Wed Jun 24 17:21:37 2009 (r1642) +++ trunk/src/Jet.cc Wed Jun 24 17:23:36 2009 (r1643) @@ -1,4 +1,5 @@ #include "Rivet/Jet.hh" +#include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/ParticleIDMethods.hh" namespace Rivet { @@ -35,33 +36,57 @@ bool Jet::containsParticle(const Particle& particle) const { const int barcode = particle.genParticle().barcode(); - foreach (const Particle& p, _fullParticles) { - const GenParticle& part = p.genParticle(); - if (part.barcode() == barcode) return true; + foreach (const Particle& p, particles()) { + if (p.genParticle().barcode() == barcode) return true; } return false; } bool Jet::containsParticleId(PdgId pid) const { - foreach (const Particle& p, _fullParticles) { + foreach (const Particle& p, particles()) { if (p.pdgId() == pid) return true; } return false; } + bool Jet::containsParticleId(vector<PdgId> pids) const { + foreach (const Particle& p, particles()) { + foreach (PdgId pid, pids) { + if (p.pdgId() == pid) return true; + } + } + return false; + } + + + /// @todo Jet::containsMatch(Matcher m) { ... if m(pid) return true; ... } + + bool Jet::containsCharm() const { - foreach (const Particle& p, _fullParticles) { - if (PID::hasCharm(p.pdgId())) return true; + foreach (const Particle& p, particles()) { + HepMC::GenVertex* gv = p.genParticle().production_vertex(); + if (gv) { + foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) { + if (PID::hasCharm(pi->pdg_id())) return true; + } + } } return false; } bool Jet::containsBottom() const { - foreach (const Particle& p, _fullParticles) { - if (PID::hasBottom(p.pdgId())) return true; + foreach (const Particle& p, particles()) { + HepMC::GenVertex* gv = p.genParticle().production_vertex(); + if (gv) { + foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) { + const PdgId pid = pi->pdg_id(); + //Log::getLog("Rivet") << Log::INFO << "PID = " << pid << endl; + if (PID::hasBottom(pid)) return true; + } + } } return false; } Modified: trunk/src/Projections/FastJets.cc ============================================================================== --- trunk/src/Projections/FastJets.cc Wed Jun 24 17:21:37 2009 (r1642) +++ trunk/src/Projections/FastJets.cc Wed Jun 24 17:23:36 2009 (r1643) @@ -46,7 +46,8 @@ _plugin.reset(new fastjet::CDFMidPointPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold)); } else if (alg == D0ILCONE) { // There's a numerical instability in the D0 IL cone which makes it really hate a pTmin of zero! - const double MIN_ET = pTmin + 1E-10; + const double MIN_ET = pTmin + (isZero(pTmin) ? 1E-10 : 0.0); + assert(MIN_ET > 1E-11); _plugin.reset(new fastjet::D0RunIIConePlugin(rparameter, MIN_ET)); } else if (alg == JADE) { _plugin.reset(new fastjet::JadePlugin()); @@ -121,7 +122,6 @@ } - Jets FastJets::_pseudojetsToJets(const PseudoJets& pjets) const { Jets rtn; foreach (const fastjet::PseudoJet& pj, pjets) { @@ -130,14 +130,8 @@ const PseudoJets parts = clusterSeq()->constituents(pj); foreach (const fastjet::PseudoJet& p, parts) { map<int, Particle>::const_iterator found = _particles.find(p.user_index()); - if (found != _particles.end()) { - // New way keeping full particle info - j.addParticle(found->second); - } else { - // Old way storing just the momentum - const FourMomentum particle(p.E(), p.px(), p.py(), p.pz()); - j.addParticle(particle); - } + assert(found != _particles.end()); + j.addParticle(found->second); } rtn.push_back(j); } @@ -145,6 +139,51 @@ } + void FastJets::reset() { + _yscales.clear(); + _particles.clear(); + /// @todo _cseq = fastjet::ClusterSequence(); + } + + + size_t FastJets::numJets(double ptmin) const { + if (_cseq.get() != 0) { + return _cseq->inclusive_jets(ptmin).size(); + } else { + return 0; + } + } + + + Jets FastJets::jets(double ptmin) const { + Jets rtn = _pseudojetsToJets(pseudoJets(ptmin)); + return rtn; + } + + + Jets FastJets::jetsByPt(double ptmin) const { + return _pseudojetsToJets(pseudoJetsByPt(ptmin)); + } + + + Jets FastJets::jetsByE(double ptmin) const { + return _pseudojetsToJets(pseudoJetsByE(ptmin)); + } + + + Jets FastJets::jetsByRapidity(double ptmin) const { + return _pseudojetsToJets(pseudoJetsByRapidity(ptmin)); + } + + + PseudoJets FastJets::pseudoJets(double ptmin) const { + if (_cseq.get() != 0) { + return _cseq->inclusive_jets(ptmin); + } else { + return PseudoJets(); + } + } + vector<double> FastJets::ySubJet(const fastjet::PseudoJet& jet) const { assert(clusterSeq()); Modified: trunk/src/Projections/FinalState.cc ============================================================================== --- trunk/src/Projections/FinalState.cc Wed Jun 24 17:21:37 2009 (r1642) +++ trunk/src/Projections/FinalState.cc Wed Jun 24 17:23:36 2009 (r1643) @@ -62,12 +62,10 @@ if (_etaRanges.empty() && _ptmin == 0) { getLog() << Log::TRACE << "Open FS processing: should only see this once per event (" << e.genEvent().event_number() << ")" << endl; - for (GenEvent::particle_const_iterator p = e.genEvent().particles_begin(); p != e.genEvent().particles_end(); ++p) { - if ((*p)->status() == 1) { - // pair<GenParticle*,GenParticle*> bps = e.genEvent().beam_particles(); - // const bool isbeam = (*p == bps.first || *p == bps.second); - // getLog() << Log::TRACE << *p << std::boolalpha << (isbeam ? " (beam)" : "") << endl; - _theParticles.push_back(Particle(**p)); + foreach (const GenParticle* p, Rivet::particles(e.genEvent())) { + if (p->status() == 1) { + //cout << "FS GV = " << p->production_vertex() << endl; + _theParticles.push_back(Particle(*p)); } } return;
More information about the Rivet-svn mailing list |