[Rivet-svn] r4222 - contrib/lhef2hepmc

[blackhole at projects.hepforge.org]

Author: buckley
Date: Tue Mar 12 15:38:32 2013
New Revision: 4222

Log:
Adding numerical safety check to xsec error calculation, thanks to Simone Pagan Griso.

Modified:
   contrib/lhef2hepmc/ChangeLog
   contrib/lhef2hepmc/lhef2hepmc.cc

Modified: contrib/lhef2hepmc/ChangeLog
==============================================================================
--- contrib/lhef2hepmc/ChangeLog	Sun Mar 10 17:24:45 2013	(r4221)
+++ contrib/lhef2hepmc/ChangeLog	Tue Mar 12 15:38:32 2013	(r4222)
@@ -1,3 +1,8 @@
+2013-03-12  Andy Buckley  <andy.buckley at cern.ch>
+
+	* Adding numerical safety check to xsec error calculation, thanks
+	to Simone Pagan Griso.
+
 2011-08-15  Andy Buckley  <andy at insectnation.org>
 
 	* Adding weights treatment to extract cross-sections appropriate

Modified: contrib/lhef2hepmc/lhef2hepmc.cc
==============================================================================
--- contrib/lhef2hepmc/lhef2hepmc.cc	Sun Mar 10 17:24:45 2013	(r4221)
+++ contrib/lhef2hepmc/lhef2hepmc.cc	Tue Mar 12 15:38:32 2013	(r4222)
@@ -86,10 +86,16 @@
       accumulated_weight_squared += weight*weight;
       event_number += 1;
       xsecval = accumulated_weight/event_number;
-      xsecerr = sqrt((accumulated_weight_squared/event_number - xsecval*xsecval)/(event_number));
-      //cout << "Estimated cross-section = " << xsecval << " +- " << xsecerr << endl;
+      // xsecerr = sqrt((accumulated_weight_squared/event_number - xsecval*xsecval)/event_number);
+      double xsecerr2 = (accumulated_weight_squared/event_number - xsecval*xsecval)/event_number;
+      if (xsecerr2 < 0) {
+        cerr << "WARNING: xsecerr^2 < 0, forcing to zero : " << xsecerr2 << endl;
+        xsecerr2 = 0.0;
+      }
+      xsecerr = sqrt(xsecerr2);
+      // cout << "Estimated cross-section = " << xsecval << " +- " << xsecerr << endl;
     } else  {
-      cout << "IDWTUP = " << idwtup << " value not handled yet. Stopping " << endl;
+      cerr << "IDWTUP = " << idwtup << " value not handled yet. Stopping " << endl;
       exit(-1);
     }
     xsec.set_cross_section(xsecval, xsecerr);