Rivet: first attempt on a Rivet interfacing plugin

- test for LPAIR/CMS run 1 analysis - test to retrieve the engine
cepgen · Feb 10, 2023 · f988085 · f988085
1 parent e4e15ba
commit f988085
Show file tree

Hide file tree

Showing 4 changed files with 175 additions and 5 deletions.
diff --git a/CepGenAddOns/RivetWrapper/CMakeLists.txt b/CepGenAddOns/RivetWrapper/CMakeLists.txt
@@ -1,15 +1,23 @@
 #--- searching for YODA
-find_library(YODA YODA HINTS $ENV{YODA_DIR} PATH_SUFFIXES lib)
-find_path(YODA_INCLUDE YODA HINTS $ENV{YODA_DIR} PATH_SUFFIXES include)
+find_library(YODA YODA HINTS $ENV{YODA_DIR} ${YODA_DIR} PATH_SUFFIXES lib)
+find_path(YODA_INCLUDE YODA HINTS $ENV{YODA_DIR} ${YODA_DIR} PATH_SUFFIXES include)
+#--- searching for Fastjet
+find_library(FASTJET NAMES fastjet HINTS $ENV{FASTJET_DIR} ${FASTJET_DIR} PATH_SUFFIXES lib)
+find_path(FASTJET_INCLUDE NAMES fastjet HINTS $ENV{FASTJET_DIR} ${FASTJET_DIR} PATH_SUFFIXES include)
+find_library(FASTJET_CONTRIB NAMES fastjetcontribfragile HINTS $ENV{FASTJET_CONTRIB_DIR} ${FASTJET_CONTRIB_DIR} PATH_SUFFIXES lib)
+#--- searching for Rivet
+find_library(RIVET NAMES Rivet HINTS $ENV{RIVET_DIR} ${RIVET_DIR} PATH_SUFFIXES lib)
+find_path(RIVET_INCLUDE NAMES Rivet HINTS $ENV{RIVET_DIR} ${RIVET_DIR} PATH_SUFFIXES include)
 
-if(NOT YODA)
+if(NOT YODA OR NOT FASTJET OR NOT RIVET)
   return()
 endif()
 
 #----- build the object
 
 cepgen_build(CepGenRivet SOURCES *.cpp
-    EXT_LIBS ${YODA}
-    EXT_HEADERS ${YODA_INCLUDE}
+    EXT_LIBS ${RIVET} ${FASTJET} ${FASTJET_CONTRIB} ${YODA} CepGenHepMC2
+    EXT_HEADERS ${RIVET_INCLUDE} ${FASTJET_INCLUDE} ${YODA_INCLUDE}
+    TESTS test/*.cc
     OPTIONS "-Wno-deprecated-copy"
     INSTALL_COMPONENT rivet)
diff --git a/CepGenAddOns/RivetWrapper/RivetAnalysisHandler.cpp b/CepGenAddOns/RivetWrapper/RivetAnalysisHandler.cpp
@@ -0,0 +1,105 @@
+/*
+ *  CepGen: a central exclusive processes event generator
+ *  Copyright (C) 2019-2023  Laurent Forthomme
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <Rivet/Analysis.hh>
+#include <Rivet/AnalysisHandler.hh>
+#include <Rivet/Tools/RivetPaths.hh>
+
+#include "CepGen/Core/Exception.h"
+#include "CepGen/Core/ParametersList.h"
+#include "CepGen/Event/Event.h"
+#include "CepGen/EventFilter/EventExporter.h"
+#include "CepGen/Modules/EventExporterFactory.h"
+#include "CepGen/Parameters.h"
+#include "CepGen/Process/Process.h"
+#include "CepGen/Utils/String.h"
+#include "CepGenAddOns/HepMC2Wrapper/HepMC2EventInterface.h"
+
+namespace cepgen {
+  /**
+     * \brief Handler for the Rivet analysis framework
+     * \author Laurent Forthomme <laurent.forthomme@cern.ch>
+     * \date Aug 2019
+     */
+  class RivetAnalysisHandler final : public EventExporter {
+  public:
+    explicit RivetAnalysisHandler(const ParametersList&);
+    ~RivetAnalysisHandler();
+
+    static ParametersDescription description() {
+      auto desc = EventExporter::description();
+      desc.setDescription("Rivet analysis handler");
+      desc.add<std::string>("filename", "output.rivet.yoda");
+      desc.add<std::vector<std::string> >("analyses", {});
+      return desc;
+    }
+
+    void initialise() override;
+    void setCrossSection(double xsec, double err_xsec) override { xsec_ = xsec, err_xsec_ = err_xsec; }
+    void operator<<(const Event&) override;
+
+  private:
+    void* enginePtr() override { return (void*)rivet_.get(); }
+
+    std::unique_ptr<Rivet::AnalysisHandler> rivet_;
+    const std::string filename_;
+    const std::vector<std::string> analyses_;
+
+    double xsec_{0.}, err_xsec_{0.};
+  };
+
+  RivetAnalysisHandler::RivetAnalysisHandler(const ParametersList& params)
+      : EventExporter(params),
+        rivet_(new Rivet::AnalysisHandler("CepGen")),
+        filename_(steer<std::string>("filename")),
+        analyses_(steer<std::vector<std::string> >("analyses")) {
+    if (analyses_.empty())
+      throw CG_FATAL("RivetAnalysisHandler") << "At least one analysis is required!";
+    for (const auto& path : params.get<std::vector<std::string> >("paths"))
+      Rivet::addAnalysisLibPath(path);
+    rivet_->addAnalyses(analyses_);
+    if (analyses_.size() != rivet_->analysesMap().size())
+      throw CG_FATAL("RivetAnalysisHandler") << "Rivet failed to find all analyses requested!\n\t"
+                                             << "You may used `rivet --list-analyses` to dump a full list.";
+  }
+
+  RivetAnalysisHandler::~RivetAnalysisHandler() {
+    rivet_->setCrossSection({xsec_, err_xsec_});
+    rivet_->finalize();
+    rivet_->writeData(filename_);
+  }
+
+  void RivetAnalysisHandler::initialise() {
+    if (!runParameters().hasProcess())
+      throw CG_FATAL("RivetAnalysisHandler") << "No process defined!";
+    if (!runParameters().process().hasEvent())
+      throw CG_FATAL("RivetAnalysisHandler")
+          << "Process \"" << runParameters().processName() << "\" has no event content!";
+    rivet_->init(HepMC::CepGenEvent(runParameters().process().event()));
+  }
+
+  void RivetAnalysisHandler::operator<<(const Event& ev) {
+    HepMC::CepGenEvent hepmc_evt(ev);
+    try {
+      rivet_->analyze(hepmc_evt);
+    } catch (const YODA::Exception& err) {
+      CG_WARNING("RivetAnalysisHandler") << "Rivet/YODA encountered the following exception:\n\t" << err.what();
+    }
+  }
+}  // namespace cepgen
+REGISTER_EXPORTER("rivet", RivetAnalysisHandler);
diff --git a/CepGenAddOns/RivetWrapper/test/rivet_wrapper.cc b/CepGenAddOns/RivetWrapper/test/rivet_wrapper.cc
@@ -0,0 +1,26 @@
+#include <Rivet/Analysis.hh>
+#include <Rivet/AnalysisHandler.hh>
+
+#include "CepGen/EventFilter/EventExporter.h"
+#include "CepGen/Generator.h"
+#include "CepGen/Modules/EventExporterFactory.h"
+#include "CepGen/Utils/Test.h"
+
+using namespace std;
+
+int main() {
+  cepgen::initialise();
+
+  vector<string> analyses = {"CMS_2011_I954992", "OPAL_1998_I474012"};
+
+  auto* rivet_wrp = cepgen::EventExporterFactory::get()
+                        .build("rivet", cepgen::ParametersList().set("analyses", analyses))
+                        .release();  // do not call the destructor (Rivet will not be initialised)
+  auto* rivet = rivet_wrp->engine<Rivet::AnalysisHandler>();
+  CG_TEST_EQUAL(rivet->analysisNames(), analyses, "List of analyses");
+  auto analysis = rivet->analysis(analyses.at(0));
+  CG_TEST_EQUAL(analysis->experiment(), "CMS", "Analysis experiment");
+  CG_TEST_EQUAL(analysis->collider(), "LHC", "Analysis collider");
+
+  CG_TEST_SUMMARY;
+}
diff --git a/CepGenAddOns/RivetWrapper/test/test_lpair.cc b/CepGenAddOns/RivetWrapper/test/test_lpair.cc
@@ -0,0 +1,31 @@
+#include <Rivet/Analysis.hh>
+#include <Rivet/AnalysisHandler.hh>
+
+#include "CepGen/EventFilter/EventExporter.h"
+#include "CepGen/Generator.h"
+#include "CepGen/Modules/EventExporterFactory.h"
+#include "CepGen/Modules/ProcessFactory.h"
+#include "CepGen/Parameters.h"
+#include "CepGen/Process/Process.h"
+
+int main() {
+  cepgen::Generator gen;
+
+  auto proc = cepgen::proc::ProcessFactory::get().build("lpair", cepgen::ParametersList());
+  gen.parametersRef().setProcess(std::move(proc));
+
+  auto rivet_wrp = cepgen::EventExporterFactory::get().build(
+      "rivet", cepgen::ParametersList().set<std::vector<std::string> >("analyses", {"CMS_2011_I954992"}));
+  auto* rivet = rivet_wrp->engine<Rivet::AnalysisHandler>();
+  gen.parametersRef().addEventExporter(std::move(rivet_wrp));
+
+  cepgen::Kinematics kin;
+  kin.incomingBeams().positive().setPdgId(2212);
+  kin.incomingBeams().negative().setPdgId(2212);
+  kin.incomingBeams().setSqrtS(7.e3);
+  gen.parametersRef().process().setKinematics(kin);
+
+  gen.integrate();
+
+  return 0;
+}