Merge branch 'sync_rivet_with_21.6' into 'master'

make Rivet_i compatible with Rivet/YODA 3.1.2/1.8.3

See merge request atlas/athena!36981

Generators/Rivet_i/examples/convert2root

+#! /usr/bin/env python
+
+from array import array
+import ROOT as rt
+import yoda, sys
+
+fName = str(sys.argv[1])
+yodaAOs = yoda.read(fName)
+rtFile = rt.TFile(fName[:fName.find('.yoda')] + '.root', 'recreate')
+for name in yodaAOs:
+  yodaAO = yodaAOs[name];  rtAO = None
+  if 'Histo1D' in str(yodaAO):
+    rtAO = rt.TH1D(name, '', yodaAO.numBins(), array('d', yodaAO.xEdges()))
+    rtAO.Sumw2(); rtErrs = rtAO.GetSumw2()
+    for i in range(rtAO.GetNbinsX()):
+      rtAO.SetBinContent(i + 1, yodaAO.bin(i).sumW())
+      rtErrs.AddAt(yodaAO.bin(i).sumW2(), i+1)
+  elif 'Scatter2D' in str(yodaAO):
+    rtAO = rt.TGraphAsymmErrors(yodaAO.numPoints())
+    for i in range(yodaAO.numPoints()):
+      x = yodaAO.point(i).x(); y = yodaAO.point(i).y()
+      xLo, xHi = yodaAO.point(i).xErrs()
+      yLo, yHi = yodaAO.point(i).yErrs()
+      rtAO.SetPoint(i, x, y)
+      rtAO.SetPointError(i, xLo, xHi, yLo, yHi)
+  else:
+    continue
+  rtAO.Write(name)
+rtFile.Close()
+

Generators/Rivet_i/examples/grid_jO.py

+theApp.EvtMax = -1
+
+import AthenaPoolCnvSvc.ReadAthenaPool
+
+from AthenaCommon.AlgSequence import AlgSequence
+job = AlgSequence()
+
+from Rivet_i.Rivet_iConf import Rivet_i
+rivet = Rivet_i()
+import os
+rivet.AnalysisPath = os.environ['PWD']
+
+rivet.Analyses += [ 'MC_JETS' ]
+rivet.RunName = ''
+rivet.HistoFile = 'MyOutput.yoda.gz'
+rivet.CrossSection = 1.0
+#rivet.IgnoreBeamCheck = True
+job += rivet
+

Generators/Rivet_i/examples/jobOptions.rivet-readevgenfile.py

-# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-
-## Example job option script to run Rivet inside Athena
-## using events read in from an EVGEN POOL/ROOT file.
-##
-## Author: James Monk <jmonk@hep.ucl.ac.uk>
-## Author: Andy Buckley <andy.buckley@cern.ch>
-
-include("GeneratorUtils/StdAnalysisSetup.py")
-theApp.EvtMax = 1000
-
-## Specify input evgen files
-import glob, AthenaPoolCnvSvc.ReadAthenaPool
-svcMgr.EventSelector.InputCollections = ["/afs/cern.ch/atlas/offline/ProdData/16.6.X/16.6.7.6/minbias-pythia8-7000.evgen.pool.root"]
-
-## Now set up Rivet
-from Rivet_i.Rivet_iConf import Rivet_i
-topAlg += Rivet_i("Rivet")
-topAlg.Rivet.Analyses = ["EXAMPLE", "MC_GENERIC"]
-topAlg.Rivet.Analyses += ["ATLAS_2012_I1082936"]
-
-## Configure ROOT file output
-from AthenaCommon.AppMgr import ServiceMgr as svcMgr
-from GaudiSvc.GaudiSvcConf import THistSvc
-svcMgr += THistSvc()
-svcMgr.THistSvc.Output = ["Rivet DATAFILE='Rivet.root' OPT='RECREATE'"]
-#svcMgr.MessageSvc.OutputLevel = ERROR

Generators/Rivet_i/examples/jobOptions.rivet.py

-# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-
-## Example job option script to run an event generator
-## and Rivet inside Athena
-##
-## Author: James Monk <jmonk@hep.ucl.ac.uk>
-## Author: Andy Buckley <andy.buckley@cern.ch>
-
-include("GeneratorUtils/StdEvgenSetup.py")
-theApp.EvtMax = 1000
-
-## Configure and add an event generator to the alg seq
-from Pythia8_i.Pythia8_iConf import Pythia8_i
-topAlg += Pythia8_i("Pythia8")
-topAlg.Pythia8.CollisionEnergy = 7000.0
-topAlg.Pythia8.Commands += ['HardQCD:all = on', 'PhaseSpace:pTHatMin = 30.0']
-
-## Now set up Rivet
-from Rivet_i.Rivet_iConf import Rivet_i
-topAlg += Rivet_i("Rivet")
-topAlg.Rivet.Analyses = ["MC_GENERIC"]
-#topAlg.Rivet.Analyses += ["ATLAS_2012_I1082936"]
-#topAlg.Rivet.DoRootHistos = False
-#topAlg.Rivet.OutputLevel = DEBUG
-
-## Configure ROOT file output
-from AthenaCommon.AppMgr import ServiceMgr as svcMgr
-from GaudiSvc.GaudiSvcConf import THistSvc
-svcMgr += THistSvc()
-svcMgr.THistSvc.Output = ["Rivet DATAFILE='Rivet.root' OPT='RECREATE'"]
-#svcMgr.MessageSvc.OutputLevel = ERROR

Generators/Rivet_i/examples/local_jO.py

+theApp.EvtMax = -1
+
+import AthenaPoolCnvSvc.ReadAthenaPool
+svcMgr.EventSelector.InputCollections = [ 'EVNT.root' ]
+
+from AthenaCommon.AlgSequence import AlgSequence
+job = AlgSequence()
+
+from Rivet_i.Rivet_iConf import Rivet_i
+rivet = Rivet_i()
+import os
+rivet.AnalysisPath = os.environ['PWD']
+
+rivet.Analyses += [ 'MC_JETS' ]
+rivet.RunName = ''
+rivet.HistoFile = 'MyOutput.yoda.gz'
+rivet.CrossSection = 1.0
+#rivet.IgnoreBeamCheck = True
+job += rivet
+

Generators/Rivet_i/src/Rivet_i.cxx

@@ -10,7 +10,6 @@
 #include "AtlasHepMC/GenEvent.h"
 
 #include "GeneratorObjects/McEventCollection.h"
-#include "GenInterfaces/IHepMCWeightSvc.h"
 #include "AthenaKernel/errorcheck.h"
 #include "PathResolver/PathResolver.h"
 
@@ -36,15 +35,14 @@
 
 Rivet_i::Rivet_i(const std::string& name, ISvcLocator* pSvcLocator) :
   AthAlgorithm(name, pSvcLocator),
-  m_hepMCWeightSvc("HepMCWeightSvc", name),
   m_analysisHandler(0),
   m_init(false)
 {
   // Options
   declareProperty("McEventKey", m_genEventKey="GEN_EVENT");
   declareProperty("Analyses", m_analysisNames);
-  declareProperty("CrossSection", m_crossSection=-1.0);
-  declareProperty("CrossSectionUncertainty", m_crossSection_uncert=-1.0);
+  declareProperty("CrossSection", m_crossSection=0.0);
+  declareProperty("CrossSectionUncertainty", m_crossSection_uncert=0.0);
   declareProperty("Stream", m_stream="/Rivet");
   declareProperty("RunName", m_runname="");
   declareProperty("HistoFile", m_file="Rivet.yoda");
@@ -53,8 +51,8 @@ Rivet_i::Rivet_i(const std::string& name, ISvcLocator* pSvcLocator) :
   declareProperty("IgnoreBeamCheck", m_ignorebeams=false);
   declareProperty("DoRootHistos", m_doRootHistos=false);
   declareProperty("SkipWeights", m_skipweights=false);
-  //declareProperty("MatchWeights", m_matchWeights="");
-  //declareProperty("UnmatchWeights", m_unmatchWeights="");
+  declareProperty("MatchWeights", m_matchWeights="");
+  declareProperty("UnmatchWeights", m_unmatchWeights="");
   declareProperty("WeightCap", m_weightcap=-1.0);
 }
 
@@ -123,7 +121,7 @@ StatusCode Rivet_i::initialize() {
   m_analysisHandler = new Rivet::AnalysisHandler(m_runname);
   assert(m_analysisHandler);
   m_analysisHandler->setIgnoreBeams(m_ignorebeams); //< Whether to do beam ID/energy consistency checks
-  m_analysisHandler->skipMultiWeights(m_skipweights); //< Whether to skip weights or not
+  m_analysisHandler->skipMultiWeights(m_skipweights); //< Only run on the nominal weight
   //m_analysisHandler->selectMultiWeights(m_matchWeights); //< Only run on a subset of the multi-weights
   //m_analysisHandler->deselectMultiWeights(m_unmatchWeights); //< Veto a subset of the multi-weights
   if (m_weightcap>0) m_analysisHandler->setWeightCap(m_weightcap);
@@ -208,11 +206,11 @@ StatusCode Rivet_i::execute() {
 StatusCode Rivet_i::finalize() {
   ATH_MSG_INFO("Rivet_i finalizing");
 
-  // Set xsec in Rivet
-  // Set xsec in Rivet
-  double custom_xs = m_crossSection > 0 ? m_crossSection : 1.0;
-  double custom_xserr = m_crossSection_uncert > 0 ? m_crossSection_uncert : 0.0;
-  m_analysisHandler->setCrossSection({custom_xs, custom_xserr});
+  // Setting cross-section in Rivet
+  // If no user-specified cross-section available,
+  // set AMI cross-section at plotting time 
+  double custom_xs = m_crossSection != 0 ? m_crossSection : 1.0;
+  m_analysisHandler->setCrossSection({custom_xs, m_crossSection_uncert});
 
   // Call Rivet finalize
   m_analysisHandler->finalize();
@@ -256,45 +254,36 @@ const HepMC::GenEvent* Rivet_i::checkEvent(const HepMC::GenEvent* event) {
 
   // weight-name cleaning
 #ifdef HEPMC3
-  std::vector<std::string>  w_names = event->weight_names();
-  std::vector<std::pair<std::string,std::string> > w_subs = {
-    {" nominal ",""},
-    {" set = ","_"},
-    {" = ","_"},
-    {"=",""},
-    {",",""},
-    {".",""},
-    {":",""},
-    {" ","_"},
-    {"#","num"},
-    {"\n","_"},
-    {"/","over"}
-  };
-  for (std::string& wname : w_names) {
-    for (const auto& sub : w_subs) {
-      size_t start_pos = wname.find(sub.first);
-      while (start_pos != std::string::npos) {
-        wname.replace(start_pos, sub.first.length(), sub.second);
-        start_pos = wname.find(sub.first);
+  std::vector<std::string>  w_wnames = event->weight_names();
+  if (w_names.size()) {
+    std::vector<std::pair<std::string,std::string> > w_subs = {
+      {" nominal ",""},
+      {" set = ","_"},
+      {" = ","_"},
+      {"=",""},
+      {",",""},
+      {".",""},
+      {":",""},
+      {" ","_"},
+      {"#","num"},
+      {"\n","_"},
+      {"/","over"}
+    };
+    for (std::string& wname : w_names) {
+      for (const auto& sub : w_subs) {
+        size_t start_pos = wname.find(sub.first);
+        while (start_pos != std::string::npos) {
+          wname.replace(start_pos, sub.first.length(), sub.second);
+          start_pos = wname.find(sub.first);
+        }
       }
     }
+    modEvent->run_info()->set_weight_names(w_names);
   }
-  modEvent->run_info()->set_weight_names(w_names);
 #else
   const HepMC::WeightContainer& old_wc = event->weights();
-  std::ostringstream stream;
-  old_wc.print(stream);
-  std::string str =  stream.str();
-  // if it only has one element, 
-  // then it doesn't use named weights
-  // --> no need for weight-name cleaning
-  if (str.size() > 1) {
-    std::vector<std::string> orig_order(m_hepMCWeightSvc->weightNames().size());
-    for (const auto& item : m_hepMCWeightSvc->weightNames()) {
-      orig_order[item.second] = item.first;
-    }
-    std::map<std::string, double> new_name_to_value;
-    std::map<std::string, std::string> old_name_to_new_name;
+  std::vector<std::string> old_wnames = old_wc.weight_names();
+  if (old_wnames.size()) {
     HepMC::WeightContainer& new_wc = modEvent->weights();
     new_wc.clear();
     std::vector<std::pair<std::string,std::string> > w_subs = {
@@ -311,33 +300,66 @@ const HepMC::GenEvent* Rivet_i::checkEvent(const HepMC::GenEvent* event) {
       {"/","over"}
     };
     std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses ()
-    for (std::sregex_iterator i = std::sregex_iterator(str.begin(), str.end(), re);
-         i != std::sregex_iterator(); ++i ) {
-      std::smatch m = *i;
-      std::vector<std::string> temp = ::split(m.str(), "[,]");
-      if (temp.size() == 2 || temp.size() == 3) {
-        std::string wname = temp[0];
-        if (temp.size() == 3)  wname += "," + temp[1];
-        std::string old_name = std::string(wname);
-        double value = old_wc[wname];
-        for (const auto& sub : w_subs) {
-          size_t start_pos = wname.find(sub.first);
-          while (start_pos != std::string::npos) {
-            wname.replace(start_pos, sub.first.length(), sub.second);
-            start_pos = wname.find(sub.first);
-          }
+    // TEMP from Rivet dev branch
+    std::vector<std::std::regex> select_patterns, deselect_patterns;
+    if (m_matchWeights != "") {
+      // Compile regex from each string in the comma-separated list
+      for (const std::string& pattern : split(m_matchWeights, ",")) {
+        select_patterns.push_back( std::regex(pattern) );
+      }
+    }
+    if (m_unmatchWeights != "") {
+      // Compile regex from each string in the comma-separated list
+      for (const std::string& pattern : split(m_unmatchWeights, ",")) {
+        deselect_patterns.push_back( std::regex(pattern) );
+      }
+    }
+    // END OF TEMP
+
+    std::map<std::string, double> new_name_to_value;
+    std::map<std::string, std::string> old_name_to_new_name;
+    for (const std::string& old_name : old_wnames) {
+      std::string wname = std::string(old_name);
+      double value = old_wc[old_name];
+      for (const auto& sub : w_subs) {
+        size_t start_pos = wname.find(sub.first);
+        while (start_pos != std::string::npos) {
+          wname.replace(start_pos, sub.first.length(), sub.second);
+          start_pos = wname.find(sub.first);
+        }
+      }
+      // Pulling some logic from the Rivet development branch
+      // until we have a release with this patch:
+
+      // Check if weight name matches a supplied string/regex and filter to select those only
+      bool match = select_patterns.empty();
+      for (const std::regex& re : select_patterns) {
+        if ( std::regex_match(wname, re) ) {
+          match = true;
+          break;
         }
-        new_name_to_value[wname] = value;
-        old_name_to_new_name[old_name] = wname;
       }
+      // Check if the remaining weight names match supplied string/regexes and *de*select accordingly
+      bool unmatch = false;
+      for (const std::regex& re : deselect_patterns) {
+        if ( std::regex_match(wname, re) ) { unmatch = true; break; }
+      }
+      if (!match || unmatch) continue;
+
+      // end of borrowing logic from the Rivet development branch
+      new_name_to_value[wname] = value;
+      old_name_to_new_name[old_name] = wname;
     }
-    for (const std::string& old_name : orig_order) {
-      const std::string& new_name = old_name_to_new_name[old_name];
+    auto itEnd = old_name_to_new_name.end();
+    for (const string& old_name : old_wnames) {
+      if (old_name_to_new_name.find(old_name) == itEnd)  continue;
+      const string& new_name = old_name_to_new_name[old_name];
       new_wc[ new_name ] = new_name_to_value[new_name];
     }
     // end of weight-name cleaning
   }
 #endif
+
 #ifdef HEPMC3
   if (
   false//!modEvent->valid_beam_particles()//FIXME!

Generators/Rivet_i/src/Rivet_i.h

@@ -6,22 +6,20 @@
 #define RIVET_I_H
 
 #include "AthenaBaseComps/AthAlgorithm.h"
-#include "GaudiKernel/ServiceHandle.h"
 
 #include "Rivet/AnalysisHandler.hh"
-#include "AtlasHepMC/GenEvent.h"
 
 #include <vector>
 #include <string>
 
 class ISvcLocator;
-class IHepMCWeightSvc;
 //class ITHistSvc;
 
 
 /// Interface to the Rivet analysis package
 /// @author James Monk <jmonk@cern.ch>
 /// @author Andy Buckley <andy.buckley@cern.ch>
+/// @author Christian Gutschow <chris.g@cern.ch>
 class Rivet_i : public AthAlgorithm {
 public:
 
@@ -56,9 +54,6 @@ private:
   /// A pointer to the THistSvc
   //ServiceHandle<ITHistSvc> m_histSvc;
 
-  /// A pointer to the HepMCWeightSvc
-  ServiceHandle<IHepMCWeightSvc> m_hepMCWeightSvc;
-
   /// The stream name for storing the output plots under (default "/Rivet")
   std::string m_stream;
 
@@ -111,10 +106,10 @@ private:
   bool m_skipweights;
 
   /// String of weight names (or regex) to select multiweights
-  //std::string m_matchWeights;
+  std::string m_matchWeights;
 
   /// String of weight names (or regex) to veto multiweights
-  //std::string m_unmatchWeights;
+  std::string m_unmatchWeights;
 
   ///Weight cap to set allowed maximum for weights 
   double m_weightcap;