From 37ef0a7c574aeee03c22e1b423a9e2cfac121bc6 Mon Sep 17 00:00:00 2001
From: Vitaliano Ciulli <vitaliano.ciulli@cern.ch>
Date: Sun, 14 Aug 2016 01:12:56 +0200
Subject: [PATCH] Fix setting cross section in Rivet analaysis

---
 .../RivetInterface/interface/RivetAnalyzer.h  |  1 +
 .../RivetInterface/plugins/RivetAnalyzer.cc   | 52 +++++++++++--------
 2 files changed, 31 insertions(+), 22 deletions(-)

diff --git a/GeneratorInterface/RivetInterface/interface/RivetAnalyzer.h b/GeneratorInterface/RivetInterface/interface/RivetAnalyzer.h
index 4f7418c01c27a..cbdf824f7f9ac 100644
--- a/GeneratorInterface/RivetInterface/interface/RivetAnalyzer.h
+++ b/GeneratorInterface/RivetInterface/interface/RivetAnalyzer.h
@@ -50,6 +50,7 @@ class RivetAnalyzer : public edm::EDAnalyzer
   std::string              _outFileName;
   bool                     _doFinalize;
   bool                     _produceDQM;
+  double                   _xsection;     
 
   DQMStore *dbe;
   std::vector<MonitorElement *> _mes;
diff --git a/GeneratorInterface/RivetInterface/plugins/RivetAnalyzer.cc b/GeneratorInterface/RivetInterface/plugins/RivetAnalyzer.cc
index 7508f25d0e82b..b15f91047be2e 100644
--- a/GeneratorInterface/RivetInterface/plugins/RivetAnalyzer.cc
+++ b/GeneratorInterface/RivetInterface/plugins/RivetAnalyzer.cc
@@ -19,7 +19,8 @@ _outFileName(pset.getParameter<std::string>("OutputFile")),
 //decide whether to finlaize tthe plots or not.
 //deciding not to finalize them can be useful for further harvesting of many jobs
 _doFinalize(pset.getParameter<bool>("DoFinalize")),
-_produceDQM(pset.getParameter<bool>("ProduceDQMOutput"))
+_produceDQM(pset.getParameter<bool>("ProduceDQMOutput")),
+_xsection(-1.)
 {
   //retrive the analysis name from paarmeter set
   std::vector<std::string> analysisNames = pset.getParameter<std::vector<std::string> >("AnalysisNames");
@@ -47,11 +48,10 @@ _produceDQM(pset.getParameter<bool>("ProduceDQMOutput"))
   std::set< AnaHandle, CmpAnaHandle >::const_iterator ibeg = analyses.begin();
   std::set< AnaHandle, CmpAnaHandle >::const_iterator iend = analyses.end();
   std::set< AnaHandle, CmpAnaHandle >::const_iterator iana; 
-  double xsection = -1.;
-  xsection = pset.getParameter<double>("CrossSection");
+  _xsection = pset.getParameter<double>("CrossSection");
   for (iana = ibeg; iana != iend; ++iana){
     if ((*iana)->needsCrossSection())
-    (*iana)->setCrossSection(xsection);
+      (*iana)->setCrossSection(_xsection);
   }
   if (_produceDQM){
     // book stuff needed for DQM
@@ -89,26 +89,34 @@ void RivetAnalyzer::analyze(const edm::Event& iEvent,const edm::EventSetup& iSet
   // get HepMC GenEvent
   const HepMC::GenEvent *myGenEvent = evt->GetEvent();
   
-  //if you want to use an external weight we have to clone the GenEvent and change the weight  
-  if ( _useExternalWeight ){
-    
+  //if you want to use an external weight or set the cross section we have to clone the GenEvent and change the weight  
+  if ( _useExternalWeight || _xsection > 0 ){
     HepMC::GenEvent * tmpGenEvtPtr = new HepMC::GenEvent( *(evt->GetEvent()) );
-    if (tmpGenEvtPtr->weights().size() == 0) {
-      throw cms::Exception("RivetAnalyzer") << "Original weight container has 0 size ";
-    }
-    if (tmpGenEvtPtr->weights().size() > 1) {
-      edm::LogWarning("RivetAnalyzer") << "Original event weight size is " << tmpGenEvtPtr->weights().size() << ". Will change only the first one ";  
-    }
+
+    if (_xsection > 0){
+      HepMC::GenCrossSection xsec;
+      xsec.set_cross_section(_xsection);
+      tmpGenEvtPtr->set_cross_section(xsec);
+    } 
+
+    if ( _useExternalWeight ){
+      if (tmpGenEvtPtr->weights().size() == 0) {
+	throw cms::Exception("RivetAnalyzer") << "Original weight container has 0 size ";
+      }
+      if (tmpGenEvtPtr->weights().size() > 1) {
+	edm::LogWarning("RivetAnalyzer") << "Original event weight size is " << tmpGenEvtPtr->weights().size() << ". Will change only the first one ";  
+      }
     
-    if(!_useLHEweights){
-      edm::Handle<GenEventInfoProduct> genEventInfoProduct;
-      iEvent.getByToken(_genEventInfoCollection, genEventInfoProduct);
-      tmpGenEvtPtr->weights()[0] = genEventInfoProduct->weight();
-    }else{
-      edm::Handle<LHEEventProduct> lheEventHandle;
-      iEvent.getByToken(_LHECollection,lheEventHandle);
-      const LHEEventProduct::WGT& wgt = lheEventHandle->weights().at(_LHEweightNumber);
-      tmpGenEvtPtr->weights()[0] = wgt.wgt;
+      if(!_useLHEweights){
+	edm::Handle<GenEventInfoProduct> genEventInfoProduct;
+	iEvent.getByToken(_genEventInfoCollection, genEventInfoProduct);
+	tmpGenEvtPtr->weights()[0] = genEventInfoProduct->weight();
+      }else{
+	edm::Handle<LHEEventProduct> lheEventHandle;
+	iEvent.getByToken(_LHECollection,lheEventHandle);
+	const LHEEventProduct::WGT& wgt = lheEventHandle->weights().at(_LHEweightNumber);
+	tmpGenEvtPtr->weights()[0] = wgt.wgt;
+      }
     }
     myGenEvent = tmpGenEvtPtr;