diff --git a/src/RivetAnalysis.cc b/src/RivetAnalysis.cc index 8e6478a..eb94b4d 100644 --- a/src/RivetAnalysis.cc +++ b/src/RivetAnalysis.cc @@ -1,108 +1,110 @@ #include "RHEJ/RivetAnalysis.hh" #ifdef RHEJ_BUILD_WITH_RIVET #include <ostream> #include "RHEJ/Event.hh" #include "yaml-cpp/yaml.h" #include "Rivet/AnalysisHandler.hh" #endif namespace RHEJ{ std::unique_ptr<Analysis> RivetAnalysis::create(YAML::Node const & config){ return std::unique_ptr<Analysis>{new RivetAnalysis{config}}; } } #ifdef RHEJ_BUILD_WITH_RIVET namespace RHEJ { 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::init(Event const & event){ rivet_runs_.push_back( {std::make_unique<Rivet::AnalysisHandler>(), "", HepMCInterface()} ); rivet_runs_.back().handler->addAnalyses(analyses_names_); if( !event.variations().empty() ){ for(auto const & vari : event.variations()){ std::ostringstream name; // calculate name ratio of the scales and use them for the output file name name << ".Scale" << vari.description->scale_name << "_MuR" << vari.description->mur_factor << "_MuF" << vari.description->muf_factor; /// @TODO replace this by proper weight name rivet_runs_.push_back( {std::make_unique<Rivet::AnalysisHandler>(), name.str(), HepMCInterface()} ); 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)); for(size_t i = 0; i < rivet_runs_.size(); ++i){ auto & run = rivet_runs_[i]; - run.handler->analyze(run.hepmc(event, i-1)); // -1: first = central + 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 RHEJ #else // no rivet => create empty analysis namespace Rivet { class AnalysisHandler {}; } namespace RHEJ { RivetAnalysis::RivetAnalysis(YAML::Node const &) { throw std::invalid_argument( "Failed to create RivetAnalysis: " "Reversed HEJ was built without rivet support" ); } void RivetAnalysis::init(Event const &){} void RivetAnalysis::fill(Event const &, Event const &){} void RivetAnalysis::finalise(){} } // namespace RHEJ #endif