diff --git a/include/RHEJ/HepMCInterface.hh b/include/RHEJ/HepMCInterface.hh
index 6b49471..d5c1c6d 100644
--- a/include/RHEJ/HepMCInterface.hh
+++ b/include/RHEJ/HepMCInterface.hh
@@ -1,47 +1,48 @@
/** \file
* \brief Header file for the HepMCInterface
*/
#pragma once
#include <cstddef>
#include "LHEF/LHEF.h"
namespace HepMC{
class GenEvent;
class GenCrossSection;
class GenRunInfo;
}
namespace RHEJ{
class Event;
//! This class converts the Events into HepMC::GenEvents
/**
* \details The output is depended on the HepMC version rHEJ is compiled with,
* both HepMC 2 and HepMC 3 are supported. If reversed HEJ is compiled
* without HepMC calling this interface will throw an error.
*
* This interface will also keep track of the cross section of all the events that
* being fed into it.
*/
class HepMCInterface{
public:
HepMCInterface();
/**
* \brief main function to convert an event into HepMC::GenEvent
*
* \param event Event to convert
- * \param weight_index optional input to choose a specific weight
+ * \param weight_index optional selection of specific weight
+ * (negative value gives central weight)
*/
- HepMC::GenEvent operator()(Event const & event, size_t weight_index = 0);
+ HepMC::GenEvent operator()(Event const & event, ssize_t weight_index = -1);
private:
size_t event_count_;
double tot_weight_;
double tot_weight2_;
void add_weight(double wt);
HepMC::GenCrossSection cross_section() const;
};
}
diff --git a/src/HepMCInterface.cc b/src/HepMCInterface.cc
index 7160ef3..8727221 100644
--- a/src/HepMCInterface.cc
+++ b/src/HepMCInterface.cc
@@ -1,144 +1,148 @@
#include "RHEJ/HepMCInterface.hh"
#include <cassert>
#ifdef RHEJ_BUILD_WITH_HepMC_VERSION
#include "RHEJ/Event.hh"
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenCrossSection.h"
#include "LHEF/LHEF.h"
namespace RHEJ{
namespace {
HepMC::FourVector to_FourVector(Particle const & sp){
return {sp.px(), sp.py(), sp.pz(), sp.E()};
}
constexpr int status_in = -1;
constexpr int status_decayed = 3;
constexpr int status_out = 1;
template<class HepMCClass, typename... Args>
auto make_ptr(Args&&... args){
#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
return HepMC::make_shared<HepMCClass>(std::forward<Args>(args)...);
#else
return new HepMCClass(std::forward<Args>(args)...);
#endif
}
} // namespace anonymous
HepMCInterface::HepMCInterface():
event_count_(0.), tot_weight_(0.), tot_weight2_(0.)
{}
void HepMCInterface::add_weight(double const wt){
++event_count_;
tot_weight_ += wt;
tot_weight2_ += wt * wt;
}
HepMC::GenCrossSection HepMCInterface::cross_section() const {
HepMC::GenCrossSection xs;
#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
xs.set_cross_section(tot_weight_, sqrt(tot_weight2_), event_count_);
/// @TODO add number of attempted events
#else // HepMC 2
xs.set_cross_section(tot_weight_, sqrt(tot_weight2_));
#endif
return xs;
}
- HepMC::GenEvent HepMCInterface::operator()(Event const & event, size_t const weight_index){
+ HepMC::GenEvent HepMCInterface::operator()(Event const & event,
+ ssize_t const weight_index){
HepMC::GenEvent out_ev{HepMC::Units::GEV, HepMC::Units::MM};
auto vx = make_ptr<HepMC::GenVertex>();
for(auto const & in: event.incoming()){
vx->add_particle_in(
make_ptr<HepMC::GenParticle>(
to_FourVector(in), static_cast<int>(in.type), status_in
)
);
}
- for(int i=0; i< (int) event.outgoing().size(); ++i){
+ for(size_t i=0; i < event.outgoing().size(); ++i){
auto const & out = event.outgoing()[i];
auto particle = make_ptr<HepMC::GenParticle>(
to_FourVector(out), static_cast<int>(out.type), status_out
);
const int status = event.decays().count(i)?status_decayed:status_out;
particle->set_status(status);
if( status == status_decayed){
auto vx_decay = make_ptr<HepMC::GenVertex>();
vx_decay->add_particle_in(particle);
for( auto const & out: event.decays().at(i)){
vx_decay->add_particle_out(
make_ptr<HepMC::GenParticle>(
to_FourVector(out), static_cast<int>(out.type), status_out
)
);
}
out_ev.add_vertex(vx_decay);
}
vx->add_particle_out(particle);
}
out_ev.add_vertex(vx);
// weights
- EventParameters central_parameters;
- if(event.variations().size() == 0){
- central_parameters = event.central();
- }
- else{
- central_parameters = event.variations(weight_index);
- }
- out_ev.weights().push_back(central_parameters.weight);
+ EventParameters event_param;
+ if(weight_index < 0)
+ event_param = event.central();
+ else if ( (size_t) weight_index < event.variations().size())
+ event_param = event.variations(weight_index);
+ else
+ throw std::invalid_argument{
+ "HepMCInterface tried to access a weight outside of the variation range."
+ };
+
+ out_ev.weights().push_back(event_param.weight);
for(auto const & var: event.variations()){
out_ev.weights().push_back(var.weight);
}
/// @TODO add name list for weights
// cross section
- add_weight(central_parameters.weight);
+ add_weight(event_param.weight);
#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
out_ev.set_cross_section(
HepMC::make_shared<HepMC::GenCrossSection>(cross_section()) );
#else // HepMC 2
out_ev.set_cross_section( cross_section() );
out_ev.set_signal_process_id(event.type()+1); // "+1": conistent with lhe
- out_ev.set_event_scale(central_parameters.mur);
+ out_ev.set_event_scale(event_param.mur);
#endif
out_ev.set_event_number(event_count_);
/// @TODO add alphaQCD (need function) and alphaQED
/// @TODO output pdf (currently not avaiable from event alone)
return out_ev;
}
}
#else // no HepMC => empty class
namespace HepMC {
class GenEvent {};
class GenCrossSection {};
}
namespace RHEJ{
HepMCInterface::HepMCInterface(){
throw std::invalid_argument(
"Failed to create HepMCInterface: "
"Reversed HEJ was built without HepMC support"
);
}
HepMC::GenEvent HepMCInterface::operator()(Event const &, size_t)
{return HepMC::GenEvent();}
void HepMCInterface::add_weight(double) {}
HepMC::GenCrossSection HepMCInterface::cross_section() const
{return HepMC::GenCrossSection();}
}
#endif
diff --git a/src/RivetAnalysis.cc b/src/RivetAnalysis.cc
index c995def..82ec78e 100644
--- a/src/RivetAnalysis.cc
+++ b/src/RivetAnalysis.cc
@@ -1,109 +1,109 @@
#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){
- if(event.variations().empty()){
- rivet_runs_.push_back(
- {std::make_unique<Rivet::AnalysisHandler>(), "", HepMCInterface()}
- );
- rivet_runs_.back().handler->addAnalyses(analyses_names_);
- }
- else{
+ rivet_runs_.push_back(
+ {std::make_unique<Rivet::AnalysisHandler>(), "", HepMCInterface()}
+ );
+ rivet_runs_.back().handler->addAnalyses(analyses_names_);
+ if( !event.variations().empty() ){
const auto & central = event.central();
for(size_t i = 0; i < event.variations().size(); ++i){
const auto & vari = event.variations(i);
std::ostringstream name;
// calculate name ratio of the scales and use them for the output file name
name << ".Scale" << i << "_MuR" << vari.mur/central.mur
<< "_MuF" << vari.muf/central.muf;
+ /// @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);
}
- for(auto & run: rivet_runs_){
- run.handler->analyze(run.hepmc(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
}
}
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