[Rivet-svn] r3289 - contrib/lhef2hepmc

[blackhole at projects.hepforge.org]

Author: buckley
Date: Mon Aug 15 14:56:59 2011
New Revision: 3289

Log:
Adding weights treatment to extract cross-sections appropriate for POWHEG, code supplied by Simone Alioli and Emanuele Re. Is this definitely usable for all LHEF generators?

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

Modified: contrib/lhef2hepmc/ChangeLog
==============================================================================
--- contrib/lhef2hepmc/ChangeLog	Mon Aug 15 12:36:24 2011	(r3288)
+++ contrib/lhef2hepmc/ChangeLog	Mon Aug 15 14:56:59 2011	(r3289)
@@ -1,6 +1,12 @@
+2011-08-15  Andy Buckley  <andy at insectnation.org>
+
+	* Adding weights treatment to extract cross-sections appropriate
+	for POWHEG, code supplied by Simone Alioli and Emanuele Re. Is
+	this definitely usable for all LHEF generators?
+
 2011-08-09  Andy Buckley  <andy at insectnation.org>
 
-	* Adding cross-section info passing.
+	* Adding direct cross-section info passing.
 
 2011-05-29  Andy Buckley  <andy at insectnation.org>
 

Modified: contrib/lhef2hepmc/LHEF.h
==============================================================================
--- contrib/lhef2hepmc/LHEF.h	Mon Aug 15 12:36:24 2011	(r3288)
+++ contrib/lhef2hepmc/LHEF.h	Mon Aug 15 14:56:59 2011	(r3289)
@@ -1545,9 +1545,12 @@
     // the event block. Save any inbetween lines. Exit if we didn't
     // find an event.
     while ( getline() && currentLine.find("<event") == std::string::npos &&
-      currentLine.find("</eventgroup>") == std::string::npos )
-      outsideBlock += currentLine + "\n";
-
+	    currentLine.find("</eventgroup>") == std::string::npos )  {
+        outsideBlock += currentLine + "\n";
+    // exit on end of Les Houches events tag found
+        if ( currentLine.find("</LesHouchesEvents>") != std::string::npos )  return false;
+    } 
+   
     if ( currentLine.find("<eventgroup") != std::string::npos ) {
       std::vector<XMLTag*> tags = XMLTag::findXMLTags(currentLine + 
 						       "</eventgroup>");

Modified: contrib/lhef2hepmc/lhef2hepmc.cc
==============================================================================
--- contrib/lhef2hepmc/lhef2hepmc.cc	Mon Aug 15 12:36:24 2011	(r3288)
+++ contrib/lhef2hepmc/lhef2hepmc.cc	Mon Aug 15 14:56:59 2011	(r3289)
@@ -16,6 +16,15 @@
 using namespace std;
 using namespace HepMC;
 
+// In case IDWTUP=+/-4 one has to keep track of the accumulated weights and event numbers
+// to evaluate the cross section on-the-fly. The last evaluation is the one used.
+// Better to be sure that crossSection() is never used to fill the histograms, but only in the finalization stage, by reweighting the histograms with crossSection()/sumOfWeights()
+
+double accumulated_weight = 0.0;
+double accumulated_weight_squared = 0.0;
+int event_number = 0;
+
+
 int main(int argc, char** argv) {
 
   // Look for a help argument
@@ -58,14 +67,31 @@
   while (reader->readEvent()) {
     GenEvent evt;
     evt.use_units(Units::GEV, Units::MM);
-    evt.weights().push_back(reader->hepeup.weight());
+    const double weight = reader->hepeup.weight();
+    evt.weights().push_back(weight);
 
     // Cross-section
     #ifdef HEPMC_HAS_CROSS_SECTION
     HepMC::GenCrossSection xsec;
-    const double xsecval = reader->heprup.XSECUP[0];
-    const double xsecerr = reader->heprup.XSECUP[1];
-    //cout << "Read cross-section = " << xsecval << " +- " << xsecerr << endl;
+    const int idwtup = reader->heprup.IDWTUP;
+
+    double xsecval = -1.0;
+    double xsecerr = -1.0;
+    if (abs(idwtup) == 3) {
+      xsecval = reader->heprup.XSECUP[0];
+      xsecerr = reader->heprup.XSECUP[1];
+      //cout << "Read cross-section = " << xsecval << " +- " << xsecerr << endl;
+    } else if (abs(idwtup) == 4) {
+      accumulated_weight += weight;
+      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;
+    } else  {
+      cout << "IDWTUP = " << idwtup << " value not handled yet. Stopping " << endl;
+      exit(-1);
+    }
     xsec.set_cross_section(xsecval, xsecerr);
     evt.set_cross_section(xsec);
     #endif