[Rivet-svn] r1643 - in trunk: . include/Rivet include/Rivet/Projections src src/Analyses src/Projections

[blackhole at projects.hepforge.org]

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;