[yoda-svn] r320 - in trunk: . include/YODA src

[blackhole at projects.hepforge.org]

Author: buckley
Date: Mon Aug 22 12:50:04 2011
New Revision: 320

Log:
First untested attempt to add binomial error combination in 1D histo division

Modified:
   trunk/TODO
   trunk/include/YODA/AnalysisObject.h
   trunk/include/YODA/Bin1D.h
   trunk/include/YODA/Exceptions.h
   trunk/include/YODA/Histo1D.h
   trunk/src/Histo1D.cc

Modified: trunk/TODO
==============================================================================
--- trunk/TODO	Mon Aug 22 12:23:29 2011	(r319)
+++ trunk/TODO	Mon Aug 22 12:50:04 2011	(r320)
@@ -4,6 +4,7 @@
 NOW
 
 * Bin division: add quadratic / binomial error treatment enum option (AB)
+  AB: done... needs to be checked. Add a test that y +- 1sigma in [0, 1]
 
 * Rebinning in 1D: merges of n adjacent bins (AB)
 

Modified: trunk/include/YODA/AnalysisObject.h
==============================================================================
--- trunk/include/YODA/AnalysisObject.h	Mon Aug 22 12:23:29 2011	(r319)
+++ trunk/include/YODA/AnalysisObject.h	Mon Aug 22 12:50:04 2011	(r320)
@@ -13,6 +13,18 @@
 
 namespace YODA {
 
+
+  /// @name Enums
+  ///
+  /// Definitions of enums used by many/all of the classes
+  //@{
+
+  enum ErrorCombination { QUAD, BINOMIAL };
+
+  //@}
+
+
+
   /// AnalysisObject is the base class for histograms and scatters
   class AnalysisObject {
 

Modified: trunk/include/YODA/Bin1D.h
==============================================================================
--- trunk/include/YODA/Bin1D.h	Mon Aug 22 12:23:29 2011	(r319)
+++ trunk/include/YODA/Bin1D.h	Mon Aug 22 12:50:04 2011	(r320)
@@ -187,6 +187,11 @@
       return _dbn.numEntries();
     }
 
+    /// The number of entries
+    unsigned long effNumEntries() const {
+      return _dbn.effNumEntries();
+    }
+
     /// The sum of weights
     double sumW() const {
       return _dbn.sumW();

Modified: trunk/include/YODA/Exceptions.h
==============================================================================
--- trunk/include/YODA/Exceptions.h	Mon Aug 22 12:23:29 2011	(r319)
+++ trunk/include/YODA/Exceptions.h	Mon Aug 22 12:50:04 2011	(r320)
@@ -23,6 +23,13 @@
   };
 
 
+  /// Error for general binning problems.
+  class BinningError : public Exception {
+  public:
+    BinningError(const std::string& what) : Exception(what) {}
+  };
+
+
   /// Error for e.g. use of invalid bin ranges.
   class RangeError : public Exception {
   public:
@@ -75,6 +82,13 @@
   };
 
 
+  /// Error for problems introduced outside YODA, to put it nicely.
+  class UserError : public Exception {
+  public:
+    UserError(const std::string& what) : Exception(what) {}
+  };
+
+
 }
 
 #endif

Modified: trunk/include/YODA/Histo1D.h
==============================================================================
--- trunk/include/YODA/Histo1D.h	Mon Aug 22 12:23:29 2011	(r319)
+++ trunk/include/YODA/Histo1D.h	Mon Aug 22 12:50:04 2011	(r320)
@@ -292,7 +292,10 @@
 
 
   /// Divide two histograms
-  Scatter2D divide(const Histo1D& numer, const Histo1D& denom);
+  ///
+  /// The
+  Scatter2D divide(const Histo1D& numer, const Histo1D& denom,
+                   ErrorCombination erropt=QUAD);
 
 
   /// Divide two histograms, with an uncorrelated error treatment

Modified: trunk/src/Histo1D.cc
==============================================================================
--- trunk/src/Histo1D.cc	Mon Aug 22 12:23:29 2011	(r319)
+++ trunk/src/Histo1D.cc	Mon Aug 22 12:50:04 2011	(r320)
@@ -110,31 +110,49 @@
   ////////////////////////////////////////
 
 
-  /// Divide two histograms
-  ///
-  /// @todo Add QUAD, BINOMIAL enum for error treatment
-  /// @todo Add named operator
-  Scatter2D divide(const Histo1D& numer, const Histo1D& denom) {
-    //assert(numer._axis == denom._axis);
+  // Divide two histograms
+  Scatter2D divide(const Histo1D& numer, const Histo1D& denom,
+                   ErrorCombination erropt) {
+    //if (!numer._axis.compatible(denom._axis)) ...
+
     Scatter2D tmp;
     for (size_t i = 0; i < numer.numBins(); ++i) {
       const HistoBin1D& b1 = numer.bin(i);
       const HistoBin1D& b2 = denom.bin(i);
-      const HistoBin1D& bL = b1 + b2;
+      const HistoBin1D bL = b1 + b2;
 
-      assert(fuzzyEquals(b1.midpoint(), b2.midpoint()));
+      /// @todo Can't we do this check with a method on Axis1D?
+      if (!fuzzyEquals(b1.midpoint(), b2.midpoint())) {
+        throw BinningError("Axis binnings are not equivalent");
+      }
 
       const double x = bL.focus();
+      /// @todo Use exceptions instead, since this is a user-inducible error
       assert(fuzzyEquals(b1.xMin(), b2.xMin()));
       assert(fuzzyEquals(b1.xMax(), b2.xMax()));
       const double exminus = x - b1.xMin();
       const double explus = b1.xMax() - x;
+      /// @todo Use exceptions instead, since this is a user-inducible error
       assert(exminus >= 0);
       assert(explus >= 0);
       //
       const double y = b1.height() / b2.height();
-      const double ey = y * sqrt( sqr(b1.heightErr()/b1.height()) + sqr(b2.heightErr()/b2.height()) );
+      double ey = -1;
+
+      switch (erropt) {
+      case QUAD:
+        ey = y * sqrt( sqr(b1.heightErr()/b1.height()) + sqr(b2.heightErr()/b2.height()) );
+        break;
+      case BINOMIAL:
+        /// @todo Check that this is correct
+        /// @todo I think this is only valid if the fills of b1 are a subset of the fills of b2. Throw an
+        /// error if Neff(b1) > Neff(b2)
+        ey = std::sqrt( b1.effNumEntries() * (1- b1.effNumEntries()/b2.effNumEntries()) )/b2.effNumEntries();
+        break;
+      }
+
       tmp.addPoint(x, exminus, explus, y, ey, ey);
+
     }
     assert(tmp.numPoints() == numer.numBins());
     return tmp;