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

[blackhole at projects.hepforge.org]

Author: buckley
Date: Tue Oct  6 16:39:28 2009
New Revision: 1879

Log:
Improving assert check in FourMomentum::mass() function in response to numerical precision fail, and taking the opportunity of merging some stats functions into MathUtils.hh since I'm anyway forcing a rebuild...

Modified:
   trunk/include/Rivet/Math/MathUtils.hh
   trunk/include/Rivet/Math/Vector4.hh
   trunk/src/Analyses/UA5_1988_S1867512.cc
   trunk/src/Run.cc

Modified: trunk/include/Rivet/Math/MathUtils.hh
==============================================================================
--- trunk/include/Rivet/Math/MathUtils.hh	Tue Oct  6 14:20:58 2009	(r1878)
+++ trunk/include/Rivet/Math/MathUtils.hh	Tue Oct  6 16:39:28 2009	(r1879)
@@ -7,6 +7,9 @@
 namespace Rivet {
 
 
+  /// @name Number comparisons etc.
+  //@{
+
   /// Compare a floating point number to zero with a degree 
   /// of fuzziness expressed by the absolute @a tolerance parameter.
   inline bool isZero(double val, double tolerance=1E-8) {
@@ -86,6 +89,54 @@
     return a*a;
   }
 
+  //@}
+
+
+
+  /// @name Statistics functions
+  //@{
+
+  /// Calculate the mean of a sample
+  inline double mean(const vector<int>& sample) {
+    double mean = 0.0;
+    foreach (const int& i, sample) {
+      mean += i;
+    }
+    return mean/sample.size();
+  }
+  
+  
+  /// Calculate the covariance (variance) between two samples  
+  inline double covariance(const vector<int>& sample1, const vector<int>& sample2) {
+    double mean1 = mean(sample1);
+    double mean2 = mean(sample2);
+    int N = sample1.size();
+    double cov = 0.0;
+    for (int i = 0; i < N; i++) {
+      double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
+      cov += cov_i;
+    }
+    if (N > 1) return cov/(N-1);
+    else return 0.0;
+  }
+  
+  
+  /// Calculate the correlation strength between two samples
+  inline double correlation(const vector<int>& sample1, const vector<int>& sample2) {
+    const double cov = covariance(sample1, sample2);
+    const double var1 = covariance(sample1, sample1);
+    const double var2 = covariance(sample2, sample2);
+    const double correlation = cov/sqrt(var1*var2);
+    const double corr_strength = correlation*sqrt(var2/var1);
+    return corr_strength;
+  }
+
+  //@}
+
+
+  /// @name Angle range mappings
+  //@{
+
   /// Reduce any number to the range [-2PI, 2PI] by repeated addition or
   /// subtraction of 2PI as required. Used to normalise angular measures.
   inline double _mapAngleM2PITo2Pi(double angle) {
@@ -123,6 +174,12 @@
     return rtn;
   }
 
+  //@}
+
+
+  /// @name Phase space measure helpers
+  //@{
+
   /// Calculate the difference between two angles in radians, 
   /// returning in the range [0, PI].
   inline double deltaPhi(double phi1, double phi2) {
@@ -152,4 +209,6 @@
 
 }
 
+//@}
+
 #endif

Modified: trunk/include/Rivet/Math/Vector4.hh
==============================================================================
--- trunk/include/Rivet/Math/Vector4.hh	Tue Oct  6 14:20:58 2009	(r1878)
+++ trunk/include/Rivet/Math/Vector4.hh	Tue Oct  6 16:39:28 2009	(r1879)
@@ -364,7 +364,7 @@
 
     /// Get mass \f$ m = \sqrt{E^2 - p^2} \f$ (the Lorentz self-invariant).
     double mass() const { 
-      assert(mass2() >= 0);
+      assert(Rivet::isZero(mass2()) || mass2() > 0);
       return sqrt(mass2()); 
     }
 
@@ -374,7 +374,7 @@
     }
 
     /// Calculate squared transverse momentum \f$ p_T^2 \f$.
-    double pT2() const { 
+    double pT2() const {
       return vector3().polarRadius2();
     }
 

Modified: trunk/src/Analyses/UA5_1988_S1867512.cc
==============================================================================
--- trunk/src/Analyses/UA5_1988_S1867512.cc	Tue Oct  6 14:20:58 2009	(r1878)
+++ trunk/src/Analyses/UA5_1988_S1867512.cc	Tue Oct  6 16:39:28 2009	(r1879)
@@ -7,49 +7,6 @@
 #include "Rivet/Projections/TriggerUA5.hh"
 
 namespace Rivet {
-
-  /// @todo Move these into the MathUtils header.
-
-  // A simple function to calculate the mean of a sample
-  double mean(const vector<int>& sample) {
-    double mean = 0.0;
-    foreach (const int& i, sample) {
-      mean += i;
-    }
-    return mean/sample.size();
-  }
-  
-  
-  // Calculate the covariance (variance) between two samples  
-  double covariance(const vector<int>& sample1, const vector<int>& sample2) {
-    double mean1 = mean(sample1);
-    double mean2 = mean(sample2);
-    int N = sample1.size();
-    double cov = 0.0;
-    for (int i = 0; i < N; i++) {
-      double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
-      cov += cov_i;
-    }
-    if (N > 1) return cov/(N-1);
-    else return 0.0;
-  }
-  
-  
-  // Calculate the correlation strength between two samples
-  double correlation(const vector<int>& sample1, const vector<int>& sample2) {
-    const double cov = covariance(sample1, sample2);
-    const double var1 = covariance(sample1, sample1);
-    const double var2 = covariance(sample2, sample2);
-    const double correlation = cov/sqrt(var1*var2);
-    const double corr_strength = correlation*sqrt(var2/var1);
-    return corr_strength;
-  }
-
-}
-
-
-
-namespace Rivet {
   
 
   class UA5_1988_S1867512 : public Analysis {

Modified: trunk/src/Run.cc
==============================================================================
--- trunk/src/Run.cc	Tue Oct  6 14:20:58 2009	(r1878)
+++ trunk/src/Run.cc	Tue Oct  6 16:39:28 2009	(r1879)
@@ -55,6 +55,8 @@
     
     GenEvent* evt = new GenEvent();
     while (io->fill_next_event(evt)) {
+
+      // Set up system based on properties of first event
       if (_numEvents == 0) {
         int num_anas_requested = _ah.analysisNames().size();
         _ah.removeIncompatibleAnalyses(beamIds(*evt));
@@ -93,7 +95,8 @@
           return false;
         }
       }
-
+      /// @todo If NOT first event, check that beams aren't changed
+      
       // Analyze event and delete HepMC event object      
       _ah.analyze(*evt);
       delete evt;