Core: event importer modules

CepGen/EventFilter/EventImporter.h

@@ -0,0 +1,52 @@
+/*
+ *  CepGen: a central exclusive processes event generator
+ *  Copyright (C) 2022-2024  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/>.
+ */
+
+#ifndef CepGen_EventFilter_EventImporter_h
+#define CepGen_EventFilter_EventImporter_h
+
+#include "CepGen/Event/Event.h"
+#include "CepGen/EventFilter/EventHandler.h"
+#include "CepGen/Utils/Value.h"
+
+namespace cepgen {
+  /// Base event importer module
+  /// \author Laurent Forthomme <laurent.forthomme@cern.ch>
+  /// \date Dec 2022
+  class EventImporter : public EventHandler {
+  public:
+    explicit EventImporter(const ParametersList& params) : EventHandler(params) {}
+
+    static ParametersDescription description() {
+      auto desc = EventHandler::description();
+      desc.setDescription("Unnamed event importer");
+      return desc;
+    }
+
+    virtual bool operator>>(Event&) const = 0;           ///< Read the next event
+    const Value& crossSection() const { return xsec_; }  ///< Process cross section and uncertainty, in pb
+
+  protected:
+    /// Specify the process cross section and uncertainty, in pb
+    void setCrossSection(const Value& xsec) { xsec_ = xsec; }
+
+  private:
+    Value xsec_;
+  };
+}  // namespace cepgen
+
+#endif

CepGen/Modules/EventImporterFactory.h

@@ -0,0 +1,41 @@
+/*
+ *  CepGen: a central exclusive processes event generator
+ *  Copyright (C) 2022  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/>.
+ */
+
+#ifndef CepGen_Modules_EventImporterFactory_h
+#define CepGen_Modules_EventImporterFactory_h
+
+#include "CepGen/Modules/ModuleFactory.h"
+
+/** \file */
+
+/// Add a generic export module definition to the factory
+#define REGISTER_EVENT_IMPORTER(name, obj)                                        \
+  namespace cepgen {                                                              \
+    struct BUILDERNM(obj) {                                                       \
+      BUILDERNM(obj)() { EventImporterFactory::get().registerModule<obj>(name); } \
+    };                                                                            \
+    static const BUILDERNM(obj) gIOImporter##obj;                                 \
+  }
+
+namespace cepgen {
+  class EventImporter;
+  /// An event import algorithm factory
+  DEFINE_FACTORY_STR(EventImporterFactory, EventImporter, "Event importers factory");
+}  // namespace cepgen
+
+#endif

CepGen/Modules/ListModules.cpp

@@ -29,6 +29,7 @@
 #include "CepGen/Modules/DerivatorFactory.h"
 #include "CepGen/Modules/DrawerFactory.h"
 #include "CepGen/Modules/EventExporterFactory.h"
+#include "CepGen/Modules/EventImporterFactory.h"
 #include "CepGen/Modules/EventModifierFactory.h"
 #include "CepGen/Modules/FormFactorsFactory.h"
 #include "CepGen/Modules/FunctionalFactory.h"
@@ -77,6 +78,7 @@ namespace cepgen {
       });
       list_modules(AlphaEMFactory::get(), "alpha(EM) evolution algorithms");
       list_modules(AlphaSFactory::get(), "alpha(s) evolution algorithms");
+      list_modules(EventImporterFactory::get(), "Event import modules");
       list_modules(GeneratorWorkerFactory::get(), "Event generation modules");
       list_modules(EventModifierFactory::get(), "Event modification modules");
       list_modules(EventExporterFactory::get(), "Export modules");

CepGen/Modules/ModuleFactory.cpp

@@ -26,6 +26,7 @@
 #include "CepGen/Core/GeneratorWorker.h"
 #include "CepGen/Event/Event.h"
 #include "CepGen/EventFilter/EventExporter.h"
+#include "CepGen/EventFilter/EventImporter.h"
 #include "CepGen/EventFilter/EventModifier.h"
 #include "CepGen/FormFactors/Parameterisation.h"
 #include "CepGen/Integration/AnalyticIntegrator.h"
@@ -137,6 +138,7 @@ namespace cepgen {
   template class ModuleFactory<Coupling, std::string>;
   template class ModuleFactory<utils::Derivator, std::string>;
   template class ModuleFactory<utils::Drawer, std::string>;
+  template class ModuleFactory<EventImporter, std::string>;
   template class ModuleFactory<EventModifier, std::string>;
   template class ModuleFactory<EventExporter, std::string>;
   template class ModuleFactory<formfac::Parameterisation, std::string>;

CepGenAddOns/HepMC2Wrapper/CMakeLists.txt

@@ -12,6 +12,7 @@ endif()
 cepgen_build(CepGenHepMC2 SOURCES *.cpp
     EXT_LIBS ${HEPMC_LIB}
     EXT_HEADERS ${HEPMC_INCLUDE}
+    TESTS test/*.cc
     INSTALL_COMPONENT hepmc2)
 cpack_add_component(hepmc2
     DISPLAY_NAME "CepGen HepMC wrappers library"

CepGenAddOns/HepMC2Wrapper/HepMC2EventInterface.cpp

@@ -33,8 +33,10 @@
 
 namespace HepMC {
   CepGenEvent::CepGenEvent(const cepgen::Event& evt) : GenEvent(Units::GEV, Units::MM) {
-    set_alphaQCD(evt.metadata.at("alphaS"));
-    set_alphaQED(evt.metadata.at("alphaEM"));
+    if (!evt.metadata.empty()) {
+      set_alphaQCD(evt.metadata.at("alphaS"));
+      set_alphaQED(evt.metadata.at("alphaEM"));
+    }
 
     weights().push_back(1.);  // unweighted events
 
@@ -50,7 +52,7 @@ namespace HepMC {
       auto part = new GenParticle(pmom, part_orig.integerPdgId(), (int)part_orig.status());
       part->set_generated_mass(cepgen::PDG::get().mass(part_orig.pdgId()));
       part->suggest_barcode(idx);
-      assoc_map_[idx].reset(part);
+      assoc_map_[idx] = part;
 
       switch (part_orig.role()) {
         case cepgen::Particle::IncomingBeam1:
@@ -100,7 +102,7 @@ namespace HepMC {
             if (!vprod) {
               vprod = new GenVertex();
               for (const auto& id : ids)
-                vprod->add_particle_in(assoc_map_.at(id).get());
+                vprod->add_particle_in(assoc_map_.at(id));
               add_vertex(vprod);
             }
             vprod->add_particle_out(part);
@@ -120,4 +122,10 @@ namespace HepMC {
       set_event_scale(evt.oneWithRole(cepgen::Particle::Role::Intermediate).momentum().mass());
     set_signal_process_vertex(vcm);
   }
+
+  CepGenEvent::operator cepgen::Event() const {
+    cepgen::Event evt;
+    print();
+    return evt;
+  }
 }  // namespace HepMC

CepGenAddOns/HepMC2Wrapper/HepMC2EventInterface.h

@@ -36,9 +36,11 @@ namespace HepMC {
   public:
     /// Construct an event interface from a CepGen Event object
     CepGenEvent(const cepgen::Event& ev);
+    /// Extract a CepGen Event object from a HepMC2 GenEvent object
+    operator cepgen::Event() const;
 
   private:
-    std::unordered_map<unsigned short, std::shared_ptr<GenParticle> > assoc_map_;
+    std::unordered_map<unsigned short, GenParticle*> assoc_map_;
   };
 }  // namespace HepMC
 #endif

CepGenAddOns/HepMC2Wrapper/HepMC2Handler.cpp

@@ -1,6 +1,6 @@
 /*
  *  CepGen: a central exclusive processes event generator
- *  Copyright (C) 2013-2023  Laurent Forthomme
+ *  Copyright (C) 2016-2024  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
@@ -59,16 +59,20 @@ namespace cepgen {
       event.set_cross_section(*xs_);
       event.set_event_number(event_num_++);
       output_->write_event(&event);
+      CG_DEBUG("HepMC2Handler").log([&event](auto& log) {
+        log << "\n";
+        event.print(log.stream());
+      });
     }
     void setCrossSection(const Value& cross_section) override {
       xs_->set_cross_section((double)cross_section, cross_section.uncertainty());
     }
 
   private:
     /// Writer object
-    std::unique_ptr<T> output_;
+    const std::unique_ptr<T> output_;
     /// Generator cross section and error
-    std::shared_ptr<GenCrossSection> xs_;
+    const std::shared_ptr<GenCrossSection> xs_;
   };
 }  // namespace cepgen
 

CepGenAddOns/HepMC2Wrapper/HepMC2Importer.cpp

@@ -0,0 +1,67 @@
+/*
+ *  CepGen: a central exclusive processes event generator
+ *  Copyright (C) 2022-2024  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 <HepMC/GenEvent.h>
+#include <HepMC/IO_GenEvent.h>
+#include <HepMC/Version.h>
+
+#include <memory>
+
+#include "CepGen/Core/Exception.h"
+#include "CepGen/EventFilter/EventImporter.h"
+#include "CepGen/Modules/EventImporterFactory.h"
+#include "CepGenAddOns/HepMC2Wrapper/HepMC2EventInterface.h"
+
+namespace cepgen {
+  /// Handler for the HepMC file output
+  /// \author Laurent Forthomme <laurent.forthomme@cern.ch>
+  /// \date Dec 2022
+  class HepMC2Importer final : public EventImporter {
+  public:
+    /// Class constructor
+    explicit HepMC2Importer(const ParametersList& params)
+        : EventImporter(params), reader_(new HepMC::IO_GenEvent(steer<std::string>("filename"), std::ios::in)) {
+      if (!reader_)
+        throw CG_FATAL("HepMC2Importer") << "Failed to initialise HepMCv2 reader.";
+      CG_INFO("HepMC2Importer") << "Interfacing module initialised "
+                                << "for HepMC version " << HEPMC_VERSION << " and HepMC ASCII file '"
+                                << steer<std::string>("filename") << "' with I/O state " << reader_->rdstate() << ".";
+    }
+
+    bool operator>>(Event& evt) const override {
+      HepMC::GenEvent event;
+      if (!reader_->fill_next_event(&event))
+        return false;
+      CG_DEBUG("HepMC2Importer").log([&event](auto& log) { event.print(log.stream()); });
+      evt = Event(static_cast<const HepMC::CepGenEvent&>(event));
+      return true;
+    }
+
+    static ParametersDescription description() {
+      auto desc = EventImporter::description();
+      desc.setDescription("HepMC2 ASCII file importer module");
+      desc.add<std::string>("filename", "input.hepmc").setDescription("Input filename");
+      return desc;
+    }
+
+  private:
+    void initialise() override {}
+    const std::unique_ptr<HepMC::IO_GenEvent> reader_;
+  };
+}  // namespace cepgen
+REGISTER_EVENT_IMPORTER("hepmc2", HepMC2Importer)

CepGenAddOns/HepMC2Wrapper/test/extract_event.cc

@@ -0,0 +1,59 @@
+#include "CepGen/Event/Event.h"
+#include "CepGen/EventFilter/EventExporter.h"
+#include "CepGen/EventFilter/EventImporter.h"
+#include "CepGen/Generator.h"
+#include "CepGen/Modules/EventExporterFactory.h"
+#include "CepGen/Modules/EventImporterFactory.h"
+#include "CepGen/Utils/ArgumentsParser.h"
+#include "CepGen/Utils/Test.h"
+
+using namespace std;
+
+int main(int argc, char* argv[]) {
+  cepgen::ArgumentsParser(argc, argv).parse();
+  cepgen::initialise();
+
+  auto evt = cepgen::Event();
+  {
+    auto p0 =
+        cepgen::Particle(cepgen::Particle::Role::IncomingBeam1, 2212, cepgen::Particle::Status::PrimordialIncoming);
+    p0.momentum().setP(0.5, 1., 1.5, 2.);
+    evt.addParticle(p0);
+    auto p1 =
+        cepgen::Particle(cepgen::Particle::Role::IncomingBeam2, 2212, cepgen::Particle::Status::PrimordialIncoming);
+    p1.momentum().setP(1., 2., 3., 4.);
+    evt.addParticle(p1);
+    auto p2 = cepgen::Particle(cepgen::Particle::Role::OutgoingBeam1, 2212, cepgen::Particle::Status::FinalState);
+    p2.momentum().setP(2., 4., 6., 8.);
+    evt.addParticle(p2);
+    p2.addMother(p0);
+    auto p3 = cepgen::Particle(cepgen::Particle::Role::OutgoingBeam2, 2212, cepgen::Particle::Status::FinalState);
+    p3.momentum().setP(4., 8., 12., 16.);
+    evt.addParticle(p3);
+    p3.addMother(p1);
+    auto p4 = cepgen::Particle(cepgen::Particle::Role::CentralSystem, 12, cepgen::Particle::Status::FinalState);
+    p4.momentum().setP(8., 16., 24., 32.);
+    evt.addParticle(p4);
+    p4.addMother(p0);
+    p4.addMother(p1);
+  }
+
+  auto temp_file = "/tmp/test_hepmc.out";
+  {
+    auto hepmc_out = cepgen::EventExporterFactory::get().build(
+        "hepmc2", cepgen::ParametersList().set<string>("filename", temp_file));
+    (*hepmc_out) << evt;
+  }
+  {
+    auto hepmc_in = cepgen::EventImporterFactory::get().build(
+        "hepmc2", cepgen::ParametersList().set<string>("filename", temp_file));
+    cepgen::Event evt_in;
+    CG_TEST_EQUAL(((*hepmc_in) >> evt_in), true, "Event re-import [HepMC2]");
+    CG_TEST_EQUAL(evt_in.size(), evt.size(), "Event re-import size");
+    for (const auto& part : evt_in.particles()) {
+      //CG_TEST_EQUAL(part.pdgId(), evt[part.id()].pdgId(), "Event re-import");
+    }
+  }
+
+  CG_TEST_SUMMARY;
+}