[Rivet-svn] r2018 - in trunk: . include/LWH

[blackhole at projects.hepforge.org]

Author: hoeth
Date: Thu Nov  5 17:05:54 2009
New Revision: 2018

Log:
Fixing error calculation in histogram division

Modified:
   trunk/ChangeLog
   trunk/include/LWH/HistogramFactory.h

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Thu Nov  5 15:24:13 2009	(r2017)
+++ trunk/ChangeLog	Thu Nov  5 17:05:54 2009	(r2018)
@@ -1,3 +1,7 @@
+2009-11-05  Hendrik Hoeth <hendrik.hoeth at cern.ch>
+
+	* Fixing histogram division code.
+
 2009-11-04  Hendrik Hoeth <hendrik.hoeth at cern.ch>
 
 	* New analysis STAR_2006_S6500200 (pion and proton pT spectra

Modified: trunk/include/LWH/HistogramFactory.h
==============================================================================
--- trunk/include/LWH/HistogramFactory.h	Thu Nov  5 15:24:13 2009	(r2017)
+++ trunk/include/LWH/HistogramFactory.h	Thu Nov  5 17:05:54 2009	(r2018)
@@ -642,30 +642,21 @@
         x->setErrorPlus(binwidth/2.0);
         x->setErrorMinus(binwidth/2.0);
 
-        double yval(0), yerrup(0), yerrdown(0);
+        double yval(0), yerr(0);
         if ( h1.binHeight(i) == 0 || h2.binHeight(j) == 0 ) {
           /// @todo Bad way of handling div by zero!
           yval = 0.0;
-          yerrup = 0.0;
-          yerrdown = 0.0;
+          yerr = 0.0;
         } else {
           yval = h1.binHeight(i) / h2.binHeight(j);
-          double h1val = h1.binHeight(i);
-          double h1min = h1val - h1.binRms(i);
-          double h1max = h1val + h1.binRms(i);
-          double h2val = h2.binHeight(j);
-          double h2min = h2.binHeight(j) - h2.binRms(j);
-          double h2max = h2.binHeight(j) + h2.binRms(j);
-          double h2invval = 1 / h2val;
-          double h2invmin = 1 / h2max;
-          double h2invmax = 1 / h2min;
-          yerrup = (h1val/h2val) * sqrt(pow((h1max-h1val)/h1val, 2) + pow((h2invmax-h2invval)/h2invval, 2));
-          yerrdown = (h1val/h2val) * sqrt(pow((h1val-h1min)/h1val, 2) + pow((h2invval-h2invmin)/h2invval, 2));
+          double y1relerr = h1.binError(i)/h1.binHeight(i);
+          double y2relerr = h2.binError(j)/h2.binHeight(j);
+          yerr = yval * sqrt(pow(y1relerr, 2) + pow(y2relerr, 2));
         }
         IMeasurement* y = point->coordinate(1);
         y->setValue(yval);
-        y->setErrorPlus(yerrup);
-        y->setErrorMinus(yerrdown);
+        y->setErrorPlus(yerr);
+        y->setErrorMinus(yerr);
       }
     }
     if ( !tree->insert(path, h) ) return 0;