diff --git a/src/RivetAnalysis.cc b/src/RivetAnalysis.cc index e3f6fb5..19f5a34 100644 --- a/src/RivetAnalysis.cc +++ b/src/RivetAnalysis.cc @@ -1,122 +1,123 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #include "HEJ/RivetAnalysis.hh" #ifdef HEJ_BUILD_WITH_RIVET #include <ostream> #include <stddef.h> #include "yaml-cpp/yaml.h" #include "Rivet/AnalysisHandler.hh" #include "HepMC/GenEvent.h" #include "HEJ/Event.hh" #include "HEJ/exceptions.hh" #endif namespace HEJ{ std::unique_ptr<Analysis> RivetAnalysis::create(YAML::Node const & config){ return std::unique_ptr<Analysis>{new RivetAnalysis{config}}; } } #ifdef HEJ_BUILD_WITH_RIVET namespace HEJ { RivetAnalysis::RivetAnalysis(YAML::Node const & config): output_name_{config["output"].as<std::string>()}, first_event_(true) { // read in analyses const auto & name_node = config["rivet"]; switch(name_node.Type()){ case YAML::NodeType::Scalar: analyses_names_.push_back(name_node.as<std::string>()); break; case YAML::NodeType::Sequence: for(YAML::const_iterator it = name_node.begin(); it != name_node.end(); ++it){ analyses_names_.push_back(it->as<std::string>()); } break; default: throw std::invalid_argument{ "No Analysis was provided to rivet. " "Either give an analysis or deactivate rivet." }; } } void RivetAnalysis::initialise(LHEF::HEPRUP const & heprup){ heprup_ = heprup; } // it is a bit ugly that we can't do this directly in `initialise` void RivetAnalysis::init(Event const & event){ rivet_runs_.reserve(event.variations().size()+1); rivet_runs_.push_back( {std::make_unique<Rivet::AnalysisHandler>(), "", HepMC2Interface(heprup_)} ); rivet_runs_.back().handler->addAnalyses(analyses_names_); if( !event.variations().empty() ){ for(auto const & vari : event.variations()){ rivet_runs_.push_back( {std::make_unique<Rivet::AnalysisHandler>(), "."+to_simple_string(*vari.description), HepMC2Interface(heprup_)} ); rivet_runs_.back().handler->addAnalyses(analyses_names_); } } } void RivetAnalysis::fill(Event const & event, Event const &){ if(first_event_){ first_event_=false; init(event); } HepMC::GenEvent hepmc_kin(rivet_runs_[0].hepmc.init_kinematics(event)); + rivet_runs_[0].hepmc.add_variation(hepmc_kin, event.variations()); for(size_t i = 0; i < rivet_runs_.size(); ++i){ auto & run = rivet_runs_[i]; run.hepmc.set_central(hepmc_kin, event, i-1); // -1: first = central run.handler->analyze(hepmc_kin); } } void RivetAnalysis::finalise(){ for(auto const & run: rivet_runs_){ run.handler->finalize(); run.handler->writeData(output_name_+run.name+std::string(".yoda")); } } } // namespace HEJ #else // no rivet => create empty analysis namespace Rivet { class AnalysisHandler {}; } namespace HEJ { RivetAnalysis::RivetAnalysis(YAML::Node const &) { throw std::invalid_argument( "Failed to create RivetAnalysis: " "HEJ 2 was built without rivet support" ); } void RivetAnalysis::initialise(LHEF::HEPRUP const &){} void RivetAnalysis::init(Event const &){} void RivetAnalysis::fill(Event const &, Event const &){} void RivetAnalysis::finalise(){} } // namespace HEJ #endif