[Rivet-svn] r1993 - in trunk: . include/Rivet include/Rivet/Projections pyext src src/Projections

[blackhole at projects.hepforge.org]

Author: buckley
Date: Wed Nov  4 01:41:01 2009
New Revision: 1993

Log:
Adding first-event beam consistency checking

Modified:
   trunk/ChangeLog
   trunk/include/Rivet/Projections/Beam.hh
   trunk/include/Rivet/Run.hh
   trunk/pyext/rivet.i
   trunk/src/Projections/Beam.cc
   trunk/src/Run.cc

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Wed Nov  4 00:46:05 2009	(r1992)
+++ trunk/ChangeLog	Wed Nov  4 01:41:01 2009	(r1993)
@@ -1,3 +1,8 @@
+2009-11-04  Andy Buckley  <andy at insectnation.org>
+
+	* Adding consistence checking on beam ID and sqrt(s) vs. those
+	from first event.
+
 2009-11-03  Andy Buckley  <andy at insectnation.org>
 
 	* Adding more assertion checks to linear algebra testing.

Modified: trunk/include/Rivet/Projections/Beam.hh
==============================================================================
--- trunk/include/Rivet/Projections/Beam.hh	Wed Nov  4 00:46:05 2009	(r1992)
+++ trunk/include/Rivet/Projections/Beam.hh	Wed Nov  4 01:41:01 2009	(r1993)
@@ -65,6 +65,7 @@
 
   ///////////////////////////////////////////////////////
 
+
   /// @name Stand-alone functions
   //@{
 
@@ -74,6 +75,9 @@
   /// Function to get beam particle IDs from an event
   BeamPair beamIds(const Event& e);
 
+  /// Function to get beam centre of mass energy from an event
+  double sqrtS(const Event& e);
+
   //@}
 
 }

Modified: trunk/include/Rivet/Run.hh
==============================================================================
--- trunk/include/Rivet/Run.hh	Wed Nov  4 00:46:05 2009	(r1992)
+++ trunk/include/Rivet/Run.hh	Wed Nov  4 01:41:01 2009	(r1993)
@@ -27,25 +27,80 @@
 
   public:
 
-    /// Get the name of this run.
-    Run& setCrossSection(const double& xs);
-    Run& setListAnalyses(const bool& );
-    
+    /// @name Set run properties
+    //@{
+
+    /// Get the cross-section for this run.
+    Run& setCrossSection(const double xs);
+
+    /// Declare whether to list available analyses
+    Run& setListAnalyses(const bool dolist);
+
+    //@}
+
+
+    /// @name Get run conditions
+    //@{
+
+    /// Get beam IDs for this run, determined from first event
+    const BeamPair& beams() const;
+
+    /// Get energy for this run, determined from first event
+    double sqrtS() const;
+
+    //@}
+
+
+    /// @name File processing stages
+    //@{
+
+    /// Set up HepMC file readers
     bool prepareFile(const std::string& evtfile);
+
+    /// Handle next event
     bool processEvent(bool firstEvent);
+
+    /// Close up HepMC I/O
     bool finalizeFile();
+
+    //@}
+
     
   private:
 
     /// AnalysisHandler object
     AnalysisHandler& _ah;
     
+    /// @name Run variables obtained from events or command line
+    //@{
+
+    /// Cross-section from command line
     double _xs;
+
+    /// Centre of mass energy, determined from first event
+    double _sqrts;
+
+    /// Beam IDs, determined from first event
+    BeamPair _beams;
+
+    //@}
+
+
+    /// Flag to show list of analyses
     bool _listAnalyses;
-    
+
+
+    /// @name HepMC I/O members
+    //@{
+
+    /// HepMC's own reader from streams   
     HepMC::IO_GenEvent* m_io;
+
+    /// STL istream, used by IO_GenEvent if input is not a file
     std::istream* m_istr;
 
+    //@}
+
   };
 
 

Modified: trunk/pyext/rivet.i
==============================================================================
--- trunk/pyext/rivet.i	Wed Nov  4 00:46:05 2009	(r1992)
+++ trunk/pyext/rivet.i	Wed Nov  4 01:41:01 2009	(r1993)
@@ -67,6 +67,8 @@
     return beamIds(Event(e));
   }
 
+  double sqrtS(const Event& e);
+
 
   class Analysis {
   public:

Modified: trunk/src/Projections/Beam.cc
==============================================================================
--- trunk/src/Projections/Beam.cc	Wed Nov  4 00:46:05 2009	(r1992)
+++ trunk/src/Projections/Beam.cc	Wed Nov  4 01:41:01 2009	(r1993)
@@ -52,5 +52,11 @@
     return beamproj.beamIDs();
   }
 
+  double sqrtS(const Event& e) {
+    Beam beamproj;
+    beamproj.project(e);
+    return beamproj.sqrtS();
+  }
+
 
 }

Modified: trunk/src/Run.cc
==============================================================================
--- trunk/src/Run.cc	Wed Nov  4 00:46:05 2009	(r1992)
+++ trunk/src/Run.cc	Wed Nov  4 01:41:01 2009	(r1993)
@@ -17,14 +17,14 @@
   }
   
   
-  Run& Run::setCrossSection(const double& xs) {
+  Run& Run::setCrossSection(const double xs) {
     _xs = xs;
     return *this;
   }
   
   
-  Run& Run::setListAnalyses(const bool& l) {
-    _listAnalyses = l;
+  Run& Run::setListAnalyses(const bool dolist) {
+    _listAnalyses = dolist;
     return *this;
   }
   
@@ -60,6 +60,28 @@
       return false;
     }
     
+    // Get beam details from first event, and ensure they match for all following events
+    if (evt->particles_size() != 0) {
+      const BeamPair beams = beamIds(*evt);
+      const double sqrts = Rivet::sqrtS(*evt);
+      Log::getLog("Rivet.Run") << Log::DEBUG << "Beams: "
+                               << beams << " @ " << sqrts/GeV << " GeV" << endl;
+      if (firstEvent) {
+        _beams = beams;
+        _sqrts = sqrts;
+        Log::getLog("Rivet.Run") << Log::INFO << "First event beams: "
+                                 << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl;
+      } else {
+        if (_beams != _beams || !fuzzyEquals(sqrts, sqrtS())) {
+          Log::getLog("Rivet.Run") << Log::ERROR << "Event beams mismatch: "
+                                   << beams << " @ " << sqrts/GeV << " GeV" << " vs. first beams "
+                                   << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl;
+          delete evt;
+          return false;
+        }
+      }
+    }
+
     // Set up system based on properties of first event
     if (firstEvent) {
       // If empty
@@ -68,8 +90,9 @@
         delete evt;
         return false;
       }
-      size_t num_anas_requested = _ah.analysisNames().size();
-      _ah.removeIncompatibleAnalyses(beamIds(*evt));
+
+      const size_t num_anas_requested = _ah.analysisNames().size();
+      _ah.removeIncompatibleAnalyses(beams());
       if (num_anas_requested > 0 && _ah.analysisNames().size() == 0) {
         Log::getLog("Rivet.Run") << Log::ERROR
             << "All analyses were incompatible with the first event's beams\n"
@@ -80,7 +103,7 @@
       
       if (_listAnalyses) {
         foreach (const std::string& ana, _ah.analysisNames()) {
-          cout<<ana<<endl;
+          cout << ana << endl;
         }
       }
     }
@@ -118,6 +141,7 @@
     return true;
   }
 
+
   bool Run::finalizeFile() {
     // Final HepMC object clean-up
     delete m_io;
@@ -126,5 +150,15 @@
     return true;
   }
 
+
+  const BeamPair& Run::beams() const {
+    return _beams;
+  }
+
+
+  double Run::sqrtS() const {
+    return _sqrts;
+  }
+
   
 }