diff --git a/src/HepMC2Interface.cc b/src/HepMC2Interface.cc
index 959f03f..89cba28 100644
--- a/src/HepMC2Interface.cc
+++ b/src/HepMC2Interface.cc
@@ -1,164 +1,168 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/HepMC2Interface.hh"
#include "HEJ/exceptions.hh"
#ifdef HEJ_BUILD_WITH_HepMC2
#include <math.h>
#include <utility>
#include "HEJ/detail/HepMCInterface_common.hh"
#include "HEJ/Event.hh"
#include "HEJ/Particle.hh"
#include "LHEF/LHEF.h"
#include "HepMC/GenCrossSection.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
namespace HEJ{
namespace detail_HepMC {
template<>
struct HepMCVersion<2> {
using GenEvent = HepMC::GenEvent;
using Beam = std::array<HepMC::GenParticle*,2>;
};
template<>
auto make_particle_ptr<2> (
Particle const & sp, int status
) {
return new HepMC::GenParticle(
to_FourVector<HepMC::FourVector>(sp),
static_cast<int> (sp.type),
status
);
}
template<>
auto make_vx_ptr<2>() {
return new HepMC::GenVertex();
}
}
HepMC2Interface::HepMC2Interface(LHEF::HEPRUP const & heprup):
event_count_(0.), tot_weight_(0.), tot_weight2_(0.)
{
beam_particle_[0] = static_cast<ParticleID>(heprup.IDBMUP.first);
beam_particle_[1] = static_cast<ParticleID>(heprup.IDBMUP.second);
beam_energy_[0] = heprup.EBMUP.first;
beam_energy_[1] = heprup.EBMUP.second;
}
HepMC::GenCrossSection HepMC2Interface::cross_section() const {
HepMC::GenCrossSection xs;
xs.set_cross_section(tot_weight_, sqrt(tot_weight2_));
return xs;
}
HepMC::GenEvent HepMC2Interface::init_kinematics(Event const & event) const {
const std::array<HepMC::GenParticle*,2> beam {
new HepMC::GenParticle(
HepMC::FourVector(0,0,-beam_energy_[0],beam_energy_[0]),
beam_particle_[0], detail_HepMC::status_beam ),
new HepMC::GenParticle(
HepMC::FourVector(0,0, beam_energy_[1],beam_energy_[1]),
beam_particle_[1], detail_HepMC::status_beam )
};
- return detail_HepMC::HepMC_init_kinematics<2>(
+ auto hepmc_ev{ detail_HepMC::HepMC_init_kinematics<2>(
event, beam, HepMC::GenEvent{ HepMC::Units::GEV, HepMC::Units::MM }
- );
+ ) };
+ hepmc_ev.weights().push_back( event.central().weight );
+ for(auto const & var: event.variations()){
+ hepmc_ev.weights().push_back( var.weight );
+ // no weight name for HepMC2 since rivet3 seem to mix them up
+ // (could be added via hepmc_ev.weights()[name]=weight)
+ }
+ return hepmc_ev;
}
void HepMC2Interface::set_central(
HepMC::GenEvent & out_ev, Event const & event, ssize_t const weight_index
){
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{
"HepMC2Interface tried to access a weight outside of the variation range."
};
const double wt = event_param.weight;
tot_weight_ += wt;
tot_weight2_ += wt * wt;
- if(out_ev.weights().size() == 0){
- out_ev.weights().push_back(wt);
- } else { // central always on first
- out_ev.weights()[0] = wt;
- }
+ // central always on first
+ assert(out_ev.weights().size() == event.variations().size()+1);
+ out_ev.weights()[0] = wt;
out_ev.set_cross_section( cross_section() );
out_ev.set_signal_process_id(event.type());
out_ev.set_event_scale(event_param.mur);
++event_count_;
out_ev.set_event_number(event_count_);
/// @TODO add alphaQCD (need function) and alphaQED
/// @TODO output pdf (currently not avaiable from event alone)
}
void HepMC2Interface::add_variation(HepMC::GenEvent & out_ev,
std::vector<EventParameters> const & varis
) const {
- for(auto const & var: varis){
- out_ev.weights().push_back(var.weight);
- }
- /// @TODO add name list for weights
+ assert(out_ev.weights().size() == varis.size()+1);
+ for(size_t i=0; i<varis.size(); ++i)
+ out_ev.weights()[i+1] = varis[i].weight; // 0 for central
}
HepMC::GenEvent HepMC2Interface::operator()(Event const & event,
ssize_t const weight_index
){
HepMC::GenEvent out_ev(init_kinematics(event));
set_central(out_ev, event, weight_index);
add_variation(out_ev, event.variations());
return out_ev;
}
}
#else // no HepMC2 => empty class
namespace HepMC {
class GenEvent {};
class GenCrossSection {};
}
namespace HEJ{
HepMC2Interface::HepMC2Interface(LHEF::HEPRUP const &){
throw std::invalid_argument(
"Failed to create HepMC2Interface: "
"HEJ 2 was built without HepMC2 support"
);
}
HepMC::GenEvent HepMC2Interface::operator()(Event const &, ssize_t)
{return HepMC::GenEvent();}
HepMC::GenEvent HepMC2Interface::init_kinematics(Event const &) const
{return HepMC::GenEvent();}
void HepMC2Interface::add_variation(HepMC::GenEvent &,
std::vector<EventParameters> const &) const {}
void HepMC2Interface::set_central(HepMC::GenEvent &, Event const &, ssize_t){}
HepMC::GenCrossSection HepMC2Interface::cross_section() const
{return HepMC::GenCrossSection();}
}
#endif
namespace HEJ{
HepMC2Interface::~HepMC2Interface() = default;
}
diff --git a/src/RivetAnalysis.cc b/src/RivetAnalysis.cc
index 21fe34c..6bd3917 100644
--- a/src/RivetAnalysis.cc
+++ b/src/RivetAnalysis.cc
@@ -1,165 +1,165 @@
/**
* \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 "Rivet/Config/RivetConfig.hh"
#ifdef RIVET_ENABLE_HEPMC_3
#include "HepMC3/GenEvent.h"
#include "HEJ/HepMC3Interface.hh"
#else // rivet with HepMC 2
#include "HepMC/GenEvent.h"
#include "HEJ/HepMC2Interface.hh"
#endif
#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 {
#ifdef RIVET_ENABLE_HEPMC_3
using HepMCInterface=HepMC3Interface;
#else
using HepMCInterface=HepMC2Interface;
#endif
struct RivetAnalysis::rivet_info {
rivet_info(std::unique_ptr<Rivet::AnalysisHandler> && h,
std::string && s, HEJ::HepMCInterface && m):
handler{std::move(h)}, name{std::move(s)}, hepmc{std::move(m)}
{}
// AnalysisHandler doesn't allow a copy constructor -> use ptr
std::unique_ptr<Rivet::AnalysisHandler> handler;
std::string name;
HEJ::HepMCInterface hepmc;
};
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){
#ifdef HEJ_USE_RIVET2
rivet_runs_.reserve(event.variations().size()+1);
#else
(void) event; // shut up compiler
#endif
rivet_runs_.emplace_back( std::make_unique<Rivet::AnalysisHandler>(),
std::string(""), HepMCInterface(heprup_)
);
rivet_runs_.back().handler->addAnalyses(analyses_names_);
#ifdef HEJ_USE_RIVET2
//! scale variation in rivet 2 through multiple analyses
if( !event.variations().empty() ){
for(auto const & vari : event.variations()){
rivet_runs_.emplace_back( std::make_unique<Rivet::AnalysisHandler>(),
"."+to_simple_string(*vari.description), HepMCInterface(heprup_)
);
rivet_runs_.back().handler->addAnalyses(analyses_names_);
}
}
#endif
}
void RivetAnalysis::fill(Event const & event, Event const &){
if(first_event_){
first_event_=false;
init(event);
}
- auto hepmc_kin(rivet_runs_[0].hepmc(event));
- rivet_runs_[0].handler->analyze(hepmc_kin);
+ auto hepmc_event(rivet_runs_.front().hepmc(event));
+ rivet_runs_.front().handler->analyze(hepmc_event);
#ifdef HEJ_USE_RIVET2
for(size_t i = 1; i < rivet_runs_.size(); ++i){
auto & run = rivet_runs_[i];
- run.hepmc.set_central(hepmc_kin, event, i-1);
- run.handler->analyze(hepmc_kin);
+ run.hepmc.set_central(hepmc_event, event, i-1);
+ run.handler->analyze(hepmc_event);
}
#endif
}
void RivetAnalysis::finalise(){
for(auto & 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 {
struct RivetAnalysis::rivet_info{};
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
namespace HEJ {
RivetAnalysis::~RivetAnalysis() = default;
}