[Rivet-svn] r2111 - in trunk: include/Rivet/Math include/Rivet/Projections src/Analyses src/Projections

[blackhole at projects.hepforge.org]

Author: buckley
Date: Mon Nov 30 12:10:02 2009
New Revision: 2111

Log:
Remove compile warning in BELLE analysis and deprecate InitialQuarks

Modified:
   trunk/include/Rivet/Math/MathUtils.hh
   trunk/include/Rivet/Projections/InitialQuarks.hh
   trunk/src/Analyses/BELLE_2006_S6265367.cc
   trunk/src/Projections/InitialQuarks.cc

Modified: trunk/include/Rivet/Math/MathUtils.hh
==============================================================================
--- trunk/include/Rivet/Math/MathUtils.hh	Sun Nov 29 17:07:36 2009	(r2110)
+++ trunk/include/Rivet/Math/MathUtils.hh	Mon Nov 30 12:10:02 2009	(r2111)
@@ -84,6 +84,24 @@
     }
   }
 
+
+  /// Determine if @a value is in the range @a low to @a high, with boundary
+  /// types defined by @a lowbound and @a highbound.
+  /// @todo Optimise to one-line at compile time?
+  inline bool inRange(long value, long low, long high,
+                      RangeBoundary lowbound=CLOSED, RangeBoundary highbound=CLOSED) {
+    if (lowbound == OPEN && highbound == OPEN) {
+      return (value > low && value < high);
+    } else if (lowbound == OPEN && highbound == CLOSED) {
+      return (value > low && value <= high);
+    } else if (lowbound == CLOSED && highbound == OPEN) {
+      return (value >= low && value < high);
+    } else { // if (lowbound == CLOSED && highbound == CLOSED) {
+      return (value >= low && value <= high);
+    }
+  }
+
+
   /// Named number-type squaring operation.
   template <typename Num>
   inline Num sqr(Num a) {

Modified: trunk/include/Rivet/Projections/InitialQuarks.hh
==============================================================================
--- trunk/include/Rivet/Projections/InitialQuarks.hh	Sun Nov 29 17:07:36 2009	(r2110)
+++ trunk/include/Rivet/Projections/InitialQuarks.hh	Mon Nov 30 12:10:02 2009	(r2111)
@@ -9,7 +9,8 @@
 namespace Rivet {
 
 
-  /// Project out all final-state particles in an event.
+  /// Project out quarks from the hard process in e+e- -> Z0 events
+  /// @deprecated This is a very dangerous and specific projection! Use e.g. PID::hasBottom and friends instead
   class InitialQuarks : public Projection {
  
   public:

Modified: trunk/src/Analyses/BELLE_2006_S6265367.cc
==============================================================================
--- trunk/src/Analyses/BELLE_2006_S6265367.cc	Sun Nov 29 17:07:36 2009	(r2110)
+++ trunk/src/Analyses/BELLE_2006_S6265367.cc	Mon Nov 30 12:10:02 2009	(r2111)
@@ -4,9 +4,7 @@
 #include "Rivet/RivetAIDA.hh"
 #include "Rivet/Tools/ParticleIdUtils.hh"
 #include "Rivet/Projections/Beam.hh"
-// #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/UnstableFinalState.hh"
-// #include "Rivet/Projections/InitialQuarks.hh"
 
 namespace Rivet {
 
@@ -25,7 +23,7 @@
 
     void analyze(const Event& e) {
       //
-      // TODO: apply BELLE hadron selection cuts
+      /// @todo Apply BELLE hadron selection cuts
       //
 
       const double weight = e.weight();
@@ -33,15 +31,13 @@
       // Loop through unstable FS particles and look for charmed mesons/baryons
       const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(e, "UFS");
 
+      /// @todo Implement sqrtS() for asymm. beams in beam projection
       const Beam beamproj = applyProjection<Beam>(e, "Beams");
       const ParticlePair& beams = beamproj.beams();
-      FourMomentum mom_tot = beamproj.beams().first.momentum() + 
-                             beamproj.beams().second.momentum();
+      FourMomentum mom_tot = beams.first.momentum() + beams.second.momentum();
       LorentzTransform cms_boost(-mom_tot.boostVector());
+      const double s = sqr(beamproj.sqrtS());
 
-      const double s = beamproj.sqrtS()*beamproj.sqrtS();
-
-      // TODO: implement sqrtS() for asymm. beams in beam projection
       const bool onresonance = true;
       // const bool onresonance = fuzzyEquals(beamproj.sqrtS(), 10.58, 1E-4);
 

Modified: trunk/src/Projections/InitialQuarks.cc
==============================================================================
--- trunk/src/Projections/InitialQuarks.cc	Sun Nov 29 17:07:36 2009	(r2110)
+++ trunk/src/Projections/InitialQuarks.cc	Mon Nov 30 12:10:02 2009	(r2111)
@@ -3,8 +3,10 @@
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Cmp.hh"
 
+
 #define IS_PARTON_PDGID(id) ( abs(id) <= 100 && abs(id) != 22 && (abs(id) < 11 || abs(id) > 18) )
 
+
 namespace Rivet {
 
 
@@ -14,58 +16,56 @@
 
 
   void InitialQuarks::project(const Event& e) {
-    Log& log = getLog();
     _theParticles.clear();
 
-    for (GenEvent::particle_const_iterator p = e.genEvent().particles_begin();
-         p != e.genEvent().particles_end(); ++p) {
-      const GenVertex* pv = (*p)->production_vertex();
-      const GenVertex* dv = (*p)->end_vertex();
-      bool passed = abs((*p)->pdg_id()) >= 1 && abs((*p)->pdg_id()) <= 5;
+    /// @todo This is all fragile and application-specific: remove!
+
+    foreach (const GenParticle* p, Rivet::particles(e.genEvent())) {
+      const GenVertex* pv = p->production_vertex();
+      const GenVertex* dv = p->end_vertex();
+      const PdgId pid = abs(p->pdg_id());
+      bool passed = (inRange(pid, 1, 5));
       if (passed) {
-        if (pv!=NULL) {
-          for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin() ;
-              pp != pv->particles_in_const_end() ; ++pp) {
+        if (pv != 0) {
+          foreach (const GenParticle* pp, particles_in(pv)) {
             // Only accept if parent is electron or Z0
-            if ( abs((*pp)->pdg_id()) != 11 && abs((*pp)->pdg_id()) != 23 )
-              passed = false;
+            const PdgId pid = abs(pp->pdg_id());
+            passed = (pid == ELECTRON || abs(pp->pdg_id()) == ZBOSON);
           }
-        }
-        else {
+        } else {
           passed = false;
         }
       }
 
-      if (log.isActive(Log::TRACE)) {
-        const int st = (*p)->status();
-        const double pT = (*p)->momentum().perp();
-        const double eta = (*p)->momentum().eta();
-        log << Log::TRACE << std::boolalpha
-            << "ID = " << (*p)->pdg_id() << ", status = " << st << ", pT = " << pT
-            << ", eta = " << eta << ": result = " << passed << endl;
-        if (pv != NULL) {
-          for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin() ;
-               pp != pv->particles_in_const_end() ; ++pp) {
-            log << Log::TRACE << std::boolalpha
-                << "     parent ID = " << (*pp)->pdg_id() << endl;
+      if (getLog().isActive(Log::TRACE)) {
+        const int st = p->status();
+        const double pT = p->momentum().perp();
+        const double eta = p->momentum().eta();
+        getLog() << Log::TRACE << std::boolalpha
+                 << "ID = " << p->pdg_id() << ", status = " << st << ", pT = " << pT
+                 << ", eta = " << eta << ": result = " << passed << endl;
+        if (pv != 0) {
+          foreach (const GenParticle* pp, particles_in(pv)) {
+            getLog() << Log::TRACE << std::boolalpha
+                     << "     parent ID = " << pp->pdg_id() << endl;
           }
         }
-        if (dv != NULL) {
-          for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin() ;
-               pp != dv->particles_out_const_end() ; ++pp) {
-            log << Log::TRACE << std::boolalpha
-                << "     child ID  = " << (*pp)->pdg_id() << endl;
+        if (dv != 0) {
+          foreach (const GenParticle* pp, particles_out(dv)) {
+            getLog() << Log::TRACE << std::boolalpha
+                     << "     child ID  = " << pp->pdg_id() << endl;
           }
         }
       }
-      if (passed) _theParticles.push_back(Particle(**p));
+      if (passed) _theParticles.push_back(Particle(*p));
     }
-    log << Log::DEBUG << "Number of initial quarks = "
-        << _theParticles.size() << endl;
-    if (not _theParticles.empty())
+    getLog() << Log::DEBUG << "Number of initial quarks = "
+             << _theParticles.size() << endl;
+    if (! _theParticles.empty())
       for (size_t i=0 ; i < _theParticles.size() ; i++)
-        log << Log::DEBUG << "Initial quark[" << i << "] = "
-            << _theParticles[i].pdgId() << std::endl;
+        getLog() << Log::DEBUG << "Initial quark[" << i << "] = "
+                 << _theParticles[i].pdgId() << std::endl;
   }
 
+
 }