Page MenuHomeHEPForge
Files F7878356

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/HepMCInterface.hh b/include/HEJ/HepMCInterface.hh
index 94c788e..d32bebd 100644
--- a/include/HEJ/HepMCInterface.hh
+++ b/include/HEJ/HepMCInterface.hh
@@ -1,72 +1,81 @@
/** \file
* \brief Header file for the HepMCInterface
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
+#include <array>
#include <sys/types.h>
#include <vector>
+#include "HEJ/PDG_codes.hh"
+
namespace HepMC{
class GenCrossSection;
class GenEvent;
}
+namespace LHEF{
+ class HEPRUP;
+}
+
namespace HEJ{
class Event;
class EventParameters;
//! This class converts the Events into HepMC::GenEvents
/**
* \details The output is depended on the HepMC version HEJ is compiled with,
* both HepMC 2 and HepMC 3 are supported. If HEJ 2 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();
+ HepMCInterface(LHEF::HEPRUP const &);
/**
* \brief main function to convert an event into HepMC::GenEvent
*
* \param event Event to convert
* \param weight_index optional selection of specific weight
* (negative value gives central weight)
*/
HepMC::GenEvent operator()(Event const & event, ssize_t weight_index = -1);
/**
* \brief initialise the event kinematics (everything but the weights)
*
* \param event Event to convert
* \param weight_index optional selection of specific weight
* (negative value gives central weight)
*/
HepMC::GenEvent init_kinematics(Event const & event);
/**
* \brief Sets the central value from \p event to \p out_ev
*
* \param out_ev HepMC::GenEvent to write to
* \param event Event to convert
* \param weight_index optional selection of specific weight
* (negative value gives "central")
*/
void set_central(HepMC::GenEvent & out_ev, Event const & event,
ssize_t weight_index = -1);
/**
* \brief Add the event \p variations to \p out_ev
*/
void add_variation(HepMC::GenEvent & out_ev,
std::vector<EventParameters> const & variations);
private:
+ std::array<ParticleID,2> beam_particle_;
+ std::array<double,2> beam_energy_;
size_t event_count_;
double tot_weight_;
double tot_weight2_;
HepMC::GenCrossSection cross_section() const;
};
}
diff --git a/include/HEJ/RivetAnalysis.hh b/include/HEJ/RivetAnalysis.hh
index cac5738..9bb70a5 100644
--- a/include/HEJ/RivetAnalysis.hh
+++ b/include/HEJ/RivetAnalysis.hh
@@ -1,68 +1,73 @@
/** \file
* \brief HEJ 2 interface to rivet analyses
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#include <vector>
+#include "LHEF/LHEF.h"
+
#include "HEJ/Analysis.hh"
#include "HEJ/HepMCInterface.hh"
namespace Rivet {
class AnalysisHandler;
}
namespace YAML {
class Node;
}
namespace HEJ {
/**
* @brief Class representing a Rivet analysis
*
* This class inherits from Analysis and can therefore be used
* like any other HEJ 2 analysis.
*/
class RivetAnalysis: public HEJ::Analysis {
public:
static std::unique_ptr<Analysis> create(YAML::Node const & config);
//! Constructor
/**
* @param config Configuration parameters
*
* config["analysis"] should be the name of a single Rivet analysis or
* a list of Rivet analyses. config["output"] is the prefix for
* the .yoda output files.
*/
RivetAnalysis(YAML::Node const & config);
+ void initialise(LHEF::HEPRUP const &) override;
//! Pass an event to the underlying Rivet analysis
void fill(HEJ::Event const & event, HEJ::Event const &) override;
bool pass_cuts(HEJ::Event const &, HEJ::Event const &) override
{return true;} //< no additional cuts are applied
void finalise() override;
private:
std::vector<std::string> analyses_names_;
const std::string output_name_;
+ LHEF::HEPRUP heprup_;
/// struct to organise the infos per rivet run/scale setting
struct rivet_info {
std::unique_ptr<Rivet::AnalysisHandler> handler;
std::string name;
HEJ::HepMCInterface hepmc;
};
std::vector<rivet_info> rivet_runs_;
/**
* \internal
- * @brief Calculates the scale variation from the first event for the output file
+ * @brief Calculates the scale variation from the first event for the output
+ * file
*/
void init(HEJ::Event const & event);
bool first_event_;
};
}
diff --git a/src/HepMCInterface.cc b/src/HepMCInterface.cc
index 0eac318..01cb081 100644
--- a/src/HepMCInterface.cc
+++ b/src/HepMCInterface.cc
@@ -1,199 +1,208 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/HepMCInterface.hh"
#include "HEJ/exceptions.hh"
#ifdef HEJ_BUILD_WITH_HepMC_VERSION
#include <math.h>
#include <utility>
#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 {
HepMC::FourVector to_FourVector(Particle const & sp){
return {sp.px(), sp.py(), sp.pz(), sp.E()};
}
constexpr int status_beam = 4;
constexpr int status_in = 11;
constexpr int status_decayed = 2;
constexpr int status_out = 1;
template<class HepMCClass, typename... Args>
auto make_ptr(Args&&... args){
#if HEJ_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():
+ HepMCInterface::HepMCInterface(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 HepMCInterface::cross_section() const {
HepMC::GenCrossSection xs;
#if HEJ_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::init_kinematics(Event const & event) {
#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
using GenParticlePtr = HepMC::GenParticlePtr;
#else
using GenParticlePtr = HepMC::GenParticle*;
#endif
HepMC::GenEvent out_ev{HepMC::Units::GEV, HepMC::Units::MM};
//! @TODO don't hardcode beam energy
#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
static
#endif
const std::array<GenParticlePtr,2> beam {
- make_ptr<HepMC::GenParticle>(HepMC::FourVector(0,0,-3500,3500),
- ParticleID::p,status_beam),
- make_ptr<HepMC::GenParticle>(HepMC::FourVector(0,0, 3500,3500),
- ParticleID::p,status_beam)
+ make_ptr<HepMC::GenParticle>(
+ HepMC::FourVector(0,0,-beam_energy_[0],beam_energy_[0]),
+ beam_particle_[0],status_beam ),
+ make_ptr<HepMC::GenParticle>(
+ HepMC::FourVector(0,0, beam_energy_[1],beam_energy_[1]),
+ beam_particle_[1],status_beam )
};
out_ev.set_beam_particles(beam[0],beam[1]);
auto vx = make_ptr<HepMC::GenVertex>();
for(size_t i=0; i<event.incoming().size(); ++i){
auto particle{ make_ptr<HepMC::GenParticle>(
to_FourVector(event.incoming()[i]),
static_cast<int>(event.incoming()[i].type), status_in
) };
auto vx_beam{ make_ptr<HepMC::GenVertex>() };
vx_beam->add_particle_in(beam[i]);
vx_beam->add_particle_out(particle);
out_ev.add_vertex(vx_beam);
vx->add_particle_in(particle);
}
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);
return out_ev;
}
void HepMCInterface::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{
"HepMCInterface 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;
}
#if HEJ_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(event_param.mur);
#endif
++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 HepMCInterface::add_variation(HepMC::GenEvent & out_ev,
std::vector<EventParameters> const & varis
) {
for(auto const & var: varis){
out_ev.weights().push_back(var.weight);
}
/// @TODO add name list for weights
}
HepMC::GenEvent HepMCInterface::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 HepMC => empty class
namespace HepMC {
class GenEvent {};
class GenCrossSection {};
}
namespace HEJ{
- HepMCInterface::HepMCInterface(){
+ HepMCInterface::HepMCInterface(LHEF::HEPRUP const &){
throw std::invalid_argument(
"Failed to create HepMCInterface: "
"HEJ 2 was built without HepMC support"
);
}
HepMC::GenEvent HepMCInterface::operator()(Event const &, ssize_t)
{return HepMC::GenEvent();}
HepMC::GenEvent HepMCInterface::init_kinematics(Event const &)
{return HepMC::GenEvent();}
void HepMCInterface::add_variation(HepMC::GenEvent &,
std::vector<EventParameters> const &){}
void HepMCInterface::set_central(HepMC::GenEvent &, Event const &, ssize_t) {}
HepMC::GenCrossSection HepMCInterface::cross_section() const
{return HepMC::GenCrossSection();}
}
#endif
diff --git a/src/HepMCWriter.cc b/src/HepMCWriter.cc
index c23014b..938a558 100644
--- a/src/HepMCWriter.cc
+++ b/src/HepMCWriter.cc
@@ -1,148 +1,148 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/HepMCWriter.hh"
#include <cassert>
#include "LHEF/LHEF.h"
#ifdef HEJ_BUILD_WITH_HepMC_VERSION
#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
#include "HepMC/LHEFAttributes.h"
#include "HepMC/WriterAscii.h"
#include "HEJ/Version.hh"
#else
#include "HepMC/IO_GenEvent.h"
#endif
#include <utility>
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/HepMCInterface.hh"
#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
namespace {
void reset_weight_info(LHEF::HEPRUP & heprup){
heprup.IDWTUP = 2;
// use placeholders for unknown init block values
// we can overwrite them after processing all events
heprup.XSECUP = {0.};
heprup.XERRUP = {0.};
heprup.XMAXUP = {0.};
}
HepMC::shared_ptr<HepMC::GenRunInfo> init_runinfo(LHEF::HEPRUP && heprup){
reset_weight_info(heprup);
HepMC::GenRunInfo runinfo;
auto hepr = HepMC::make_shared<HepMC::HEPRUPAttribute>();
hepr->heprup = heprup;
runinfo.add_attribute(std::string("HEPRUP"), hepr);
for (int i = 0, N = hepr->heprup.generators.size(); i < N; ++i ){
HepMC::GenRunInfo::ToolInfo tool;
tool.name = hepr->heprup.generators[i].name;
tool.version = hepr->heprup.generators[i].version;
tool.description = hepr->heprup.generators[i].contents;
runinfo.tools().push_back(tool);
}
return HepMC::make_shared<HepMC::GenRunInfo>(runinfo);
}
} // namespace anonymous
#endif // HepMC 3
namespace HEJ{
struct HepMCWriter::HepMCWriterImpl{
HepMCInterface hepmc_;
HepMCWriterImpl & operator=(HepMCWriterImpl const & other) = delete;
HepMCWriterImpl(HepMCWriterImpl const & other) = delete;
HepMCWriterImpl & operator=(HepMCWriterImpl && other) = delete;
HepMCWriterImpl(HepMCWriterImpl && other) = delete;
#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
HepMC::WriterAscii writer_;
HepMCWriterImpl(
std::string const & file, LHEF::HEPRUP && heprup
):
- hepmc_(),
+ hepmc_(heprup),
writer_{file, init_runinfo(std::move(heprup))}
{}
~HepMCWriterImpl(){
writer_.close();
}
#else // HepMC 2
HepMC::IO_GenEvent writer_;
HepMCWriterImpl(
- std::string const & file, LHEF::HEPRUP &&
+ std::string const & file, LHEF::HEPRUP && heprup
):
- hepmc_(),
+ hepmc_(heprup),
writer_{file}
{}
#endif
void write(Event const & ev){
auto out_ev = hepmc_(ev);
#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
writer_.write_event(out_ev);
#else // HepMC 2
writer_.write_event(&out_ev);
#endif
}
};
void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
delete p;
}
HepMCWriter::HepMCWriter(std::string const & file, LHEF::HEPRUP heprup):
impl_{std::unique_ptr<HepMCWriterImpl, HepMCWriterImplDeleter>{
new HepMCWriterImpl(file, std::move(heprup))
}}
{}
void HepMCWriter::write(Event const & ev){
impl_->write(ev);
}
} // namespace HEJ
#else // no HepMC
namespace HEJ{
class HepMCWriter::HepMCWriterImpl{};
HepMCWriter::HepMCWriter(std::string const &, LHEF::HEPRUP){
throw std::invalid_argument(
"Failed to create HepMC writer: "
"HEJ 2 was built without HepMC support"
);
}
void HepMCWriter::write(Event const &){
assert(false);
}
void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
delete p;
}
}
#endif
diff --git a/src/RivetAnalysis.cc b/src/RivetAnalysis.cc
index cb07fd9..9f1fef2 100644
--- a/src/RivetAnalysis.cc
+++ b/src/RivetAnalysis.cc
@@ -1,143 +1,148 @@
/**
* \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."
};
}
}
namespace {
void replace(
std::string & str,
char to_replace, std::string const & replacement
) {
for(
auto pos = str.find(to_replace);
pos != str.npos;
pos = str.find(to_replace, pos)
) {
str.replace(pos, 1, replacement);
pos += replacement.size();
}
}
// remove "special" characters from scale name
// so that we can more easily use it as part of a file name
std::string sanitise_scalename(std::string scalename) {
replace(scalename, '/', "_over_");
replace(scalename, '*', "_times_");
return scalename;
}
}
+ void RivetAnalysis::initialise(LHEF::HEPRUP const & heprup){
+ heprup_ = heprup;
+ }
+
void RivetAnalysis::init(Event const & event){
+ rivet_runs_.reserve(event.variations().size()+1);
rivet_runs_.push_back(
- {std::make_unique<Rivet::AnalysisHandler>(), "", HepMCInterface()}
+ {std::make_unique<Rivet::AnalysisHandler>(), "", HepMCInterface(heprup_)}
);
rivet_runs_.back().handler->addAnalyses(analyses_names_);
if( !event.variations().empty() ){
- rivet_runs_.reserve(event.variations().size()+1);
for(auto const & vari : event.variations()){
std::ostringstream name;
name << ".Scale" << sanitise_scalename(vari.description->scale_name)
<< "_MuR" << vari.description->mur_factor
<< "_MuF" << vari.description->muf_factor;
rivet_runs_.push_back(
- {std::make_unique<Rivet::AnalysisHandler>(), name.str(), HepMCInterface()}
+ {std::make_unique<Rivet::AnalysisHandler>(), name.str(),
+ HepMCInterface(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));
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::init(Event const &){}
void RivetAnalysis::fill(Event const &, Event const &){}
void RivetAnalysis::finalise(){}
} // namespace HEJ
#endif

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 5:57 PM (1 d, 15 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805469
Default Alt Text
(20 KB)

Event Timeline