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

[blackhole at projects.hepforge.org]

Author: hoeth
Date: Mon Nov 30 12:57:23 2009
New Revision: 2112

Log:
Fixing division by zero in Profile1D bin errors for bins with just a single entry.

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

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Mon Nov 30 12:10:02 2009	(r2111)
+++ trunk/ChangeLog	Mon Nov 30 12:57:23 2009	(r2112)
@@ -1,3 +1,8 @@
+2009-11-30  Hendrik Hoeth <hendrik.hoeth at cern.ch>
+
+	* Fixing division by zero in Profile1D bin errors for
+	bins with just a single entry.
+
 2009-11-24  Hendrik Hoeth <hendrik.hoeth at cern.ch>
 
 	* First working version of STAR_2006_S6860818

Modified: trunk/include/LWH/Profile1D.h
==============================================================================
--- trunk/include/LWH/Profile1D.h	Mon Nov 30 12:10:02 2009	(r2111)
+++ trunk/include/LWH/Profile1D.h	Mon Nov 30 12:57:23 2009	(r2112)
@@ -246,12 +246,12 @@
     bool first = true;
     for ( int i = 3; i < ax->bins() + 2; ++i ) {
       if (sumw[i] > 0.) {
-	double yw = sumyw[i]/sumw[i];
-	if (first) {
-	  minw = yw;
-	  first = false;
-	}
-	else if (yw < minw) minw = yw;
+        double yw = sumyw[i]/sumw[i];
+        if (first) {
+          minw = yw;
+          first = false;
+        }
+        else if (yw < minw) minw = yw;
       }
     }
     return minw;
@@ -267,12 +267,12 @@
     bool first = true;
     for ( int i = 3; i < ax->bins() + 2; ++i ) {
       if (sumw[i] > 0.) {
-	double yw = sumyw[i]/sumw[i];
-	if (first) {
-	  maxw = yw;
-	  first = false;
-	}
-	else if (yw > maxw) maxw = yw;
+        double yw = sumyw[i]/sumw[i];
+        if (first) {
+          maxw = yw;
+          first = false;
+        }
+        else if (yw > maxw) maxw = yw;
       }
     }
     return maxw;
@@ -352,6 +352,9 @@
    */
   double binError(int index) const {
     if (sumw[index+2] > 0.0) {
+      if ((sumw[index+2]*sumw[index+2] - sumw2[index+2]) == 0) {
+        return sumyw[index+2]/sumw[index+2];
+      }
       double binErr2 = sumy2w[index+2]*sumw[index+2] - sumyw[index+2]*sumyw[index+2];
       binErr2 /= sumw[index+2]*sumw[index+2] - sumw2[index+2];
       binErr2 /= sumw[index+2]; //< s_hat ~ s/sqrt(N)
@@ -417,8 +420,8 @@
    */
   bool add(const Profile1D & h) {
     if ( ax->upperEdge() != h.ax->upperEdge() ||
-	 ax->lowerEdge() != h.ax->lowerEdge() ||
-	 ax->bins() != h.ax->bins() ) return false;
+         ax->lowerEdge() != h.ax->lowerEdge() ||
+         ax->bins() != h.ax->bins() ) return false;
     for ( int i = 0; i < ax->bins() + 2; ++i ) {
       sum[i] += h.sum[i];
       sumw[i] += h.sumw[i];
@@ -481,7 +484,7 @@
     if ( vax ) {
       os << ">\n";
       for ( int i = 0, N = ax->bins() - 1; i < N; ++i )
-	os << "      <binBorder value=\"" << ax->binUpperEdge(i) << "\"/>\n";
+        os << "      <binBorder value=\"" << ax->binUpperEdge(i) << "\"/>\n";
       os << "    </axis>\n";
     } else {
       os << "/>\n";
@@ -519,8 +522,8 @@
        << " \"" << title() << " \"" << std::endl;
     for ( int i = 2; i < ax->bins() + 2; ++i )
       if ( sum[i] && binError(i)>0.)
-	os << binMean(i - 2) << " "
-	   << binHeight(i) << " " << binError(i) << " " << sum[i] << std::endl;
+        os << binMean(i - 2) << " "
+           << binHeight(i) << " " << binError(i) << " " << sum[i] << std::endl;
     os << std::endl;
     return true;
   }
@@ -546,7 +549,7 @@
       nbins = vax->bins();
       double* bins = new double[nbins+1];
       for (int i=0; i<nbins; ++i) {
-	bins[i] = vax->binEdges(i).first;
+        bins[i] = vax->binEdges(i).first;
       }
       bins[nbins] = vax->binEdges(nbins-1).second; //take last bin right border
       prof1d = new TProfile(name.c_str(), title().c_str(), nbins, bins);
@@ -557,14 +560,14 @@
     double entries = 0;
     for ( int i = 0; i < nbins + 2; ++i ) {
       if ( sum[i] && binError(i)>0.) {
-	//i==0: underflow->RootBin(0), i==1: overflow->RootBin(NBins+1)
-	entries = entries + sum[i];
-	int j=i;
-	if (i==0) j=0; //underflow
-	else if (i==1) j=nbins+1; //overflow
-	if (i>=2) j=i-1; //normal bin entries
-	prof1d->SetBinContent(j, binHeight(i));
-	prof1d->SetBinError(j, binError(i));
+        //i==0: underflow->RootBin(0), i==1: overflow->RootBin(NBins+1)
+        entries = entries + sum[i];
+        int j=i;
+        if (i==0) j=0; //underflow
+        else if (i==1) j=nbins+1; //overflow
+        if (i>=2) j=i-1; //normal bin entries
+        prof1d->SetBinContent(j, binHeight(i));
+        prof1d->SetBinError(j, binError(i));
       }
     }