Page MenuHomeHEPForge
Files F7877900

No OneTemporary
Actions

View Options

diff --git a/include/RHEJ/Event.hh b/include/RHEJ/Event.hh
index fa795d7..77dd373 100644
--- a/include/RHEJ/Event.hh
+++ b/include/RHEJ/Event.hh
@@ -1,129 +1,136 @@
/** \file Event.hh
* \brief Details the parameters of an Event.
*
* Contains the EventParameters struct and the Event struct. Also
* contains the ClusteredEvent class.
*/
#pragma once
#include <string>
+#include <unordered_map>
#include "RHEJ/utility.hh"
#include "RHEJ/event_types.hh"
#include "LHEF/LHEF.h"
#include "fastjet/JetDefinition.hh"
#include "fastjet/ClusterSequence.hh"
namespace RHEJ{
/** \struct EventParameters Event.hh "include/RHEJ/Event.hh"
* \brief A struct containing the parameters of an event
*/
struct EventParameters{
double mur; /**< Value of the Renormalisation Scale */
double muf; /**< Value of the Factorisation Scale */
double weight; /**< Event Weight */
};
/** \struct UnclusteredEvent Event.hh "include/RHEJ/Event.hh"
* \brief Event Struct which contains Particle Content.
*/
struct UnclusteredEvent{
//! Default Constructor
UnclusteredEvent() = default;
//! Constructor from LesHouches event information
UnclusteredEvent(LHEF::HEPEUP const & hepeup);
std::array<Sparticle, 2> incoming; /**< Incoming Particles */
std::vector<Sparticle> outgoing; /**< Outgoing Particles */
+ //! Particle decays in the format {outgoing index, decay products}
+ std::unordered_map<int, std::vector<Sparticle>> decays;
//! Central parameter (e.g. scale) choice
EventParameters central;
std::vector<EventParameters> variations; /**< For parameter variation */
};
/** \class Event Event.hh "include/RHEJ/Event.hh"
* \brief Class for Events
*
* This is the default reversed HEJ event class.
* It contains kinematic information including jet clustering,
* parameter (e.g. scale) settings and the event weight.
*/
class Event{
public:
//! Default Event Constructor
Event() = default;
//! Event Constructor adding jet clustering to an unclustered event
Event(
UnclusteredEvent ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
);
std::vector<fastjet::PseudoJet> jets() const;
UnclusteredEvent const & unclustered() const {
return ev_;
}
EventParameters const & central() const{
return ev_.central;
}
EventParameters & central(){
return ev_.central;
}
std::array<Sparticle, 2> const & incoming() const{
return ev_.incoming;
}
std::vector<Sparticle> const & outgoing() const{
return ev_.outgoing;
}
+ std::unordered_map<int, std::vector<Sparticle>> const & decays() const{
+ return ev_.decays;
+ }
+
std::vector<EventParameters> const & variations() const{
return ev_.variations;
}
std::vector<EventParameters> & variations(){
return ev_.variations;
}
EventParameters const & variations(size_t i) const{
return ev_.variations[i];
}
EventParameters & variations(size_t i){
return ev_.variations[i];
}
std::vector<int> particle_jet_indices(
std::vector<fastjet::PseudoJet> const & jets
) const{
return cs_.particle_jet_indices(jets);
}
fastjet::JetDefinition const & jet_def() const{
return cs_.jet_def();
}
double min_jet_pt() const{
return min_jet_pt_;
}
event_type::EventType type() const{
return type_;
}
private:
UnclusteredEvent ev_;
fastjet::ClusterSequence cs_;
double min_jet_pt_;
event_type::EventType type_;
};
double shat(Event const & ev);
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP *);
}
diff --git a/include/RHEJ/exceptions.hh b/include/RHEJ/exceptions.hh
index fbe70c4..f1ab441 100644
--- a/include/RHEJ/exceptions.hh
+++ b/include/RHEJ/exceptions.hh
@@ -1,44 +1,53 @@
/** \file exceptions.hh
* \brief Custom exception classes
*/
#pragma once
#include <stdexcept>
namespace RHEJ{
//! Exception indicating wrong option type
/**
* This exception is thrown if a configuration option has
* the wrong type (e.g. 'trials' is not set to a number)
*/
struct invalid_type: std::invalid_argument {
explicit invalid_type(std::string const & what):
std::invalid_argument{what} {};
explicit invalid_type(char const * what):
std::invalid_argument{what} {};
};
//! Exception indicating unknown option
/**
* This exception is thrown if an unknown configuration option
* is set (e.g. the 'trials' setting is misspelt as 'trails')
*/
struct unknown_option: std::invalid_argument {
explicit unknown_option(std::string const & what):
std::invalid_argument{what} {};
explicit unknown_option(char const * what):
std::invalid_argument{what} {};
};
//! Exception indicating missing option setting
/**
* This exception is thrown if a mandatory configuration option
* (e.g. 'trials') is not set.
*/
struct missing_option: std::logic_error {
explicit missing_option(std::string const & what):
std::logic_error{what} {};
explicit missing_option(char const * what):
std::logic_error{what} {};
};
+
+ //! Exception indicating functionality that has not been implemented yet
+ struct not_implemented: std::logic_error {
+ explicit not_implemented(std::string const & what):
+ std::logic_error{what} {};
+ explicit not_implemented(char const * what):
+ std::logic_error{what} {};
+ };
+
}
diff --git a/src/Event.cc b/src/Event.cc
index 173e2c2..d4a4fb0 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,299 +1,342 @@
#include "RHEJ/Event.hh"
#include "RHEJ/debug.hh"
namespace RHEJ{
namespace{
constexpr int status_in = -1;
+ constexpr int status_decayed = 2;
constexpr int status_out = 1;
// helper functions to determine event type
// check if there is at most one photon, W, H, Z in the final state
// and all the rest are quarks or gluons
bool final_state_ok(std::vector<Sparticle> const & outgoing){
bool has_AWZH_boson = false;
for(auto const & out: outgoing){
if(is_AWZH_boson(out.type)){
if(has_AWZH_boson) return false;
has_AWZH_boson = true;
}
else if(! is_parton(out.type)) return false;
}
return true;
}
template<class Iterator>
Iterator remove_AWZH(Iterator begin, Iterator end){
return std::remove_if(
begin, end, [](Sparticle const & p){return is_AWZH_boson(p);}
);
}
template<class Iterator>
bool valid_outgoing(Iterator begin, Iterator end){
return std::distance(begin, end) >= 2
&& std::is_sorted(begin, end, rapidity_less{})
&& std::count_if(
begin, end, [](Sparticle const & s){return is_AWZH_boson(s);}
) < 2;
}
// Note that this changes the outgoing range!
template<class ConstIterator, class Iterator>
bool is_FKL(
ConstIterator begin_incoming, ConstIterator end_incoming,
Iterator begin_outgoing, Iterator end_outgoing
){
assert(std::distance(begin_incoming, end_incoming) == 2);
assert(std::distance(begin_outgoing, end_outgoing) >= 2);
// One photon, W, H, Z in the final state is allowed.
// Remove it for remaining tests,
end_outgoing = remove_AWZH(begin_outgoing, end_outgoing);
// Test if this is a standard FKL configuration.
return
(begin_incoming->type == begin_outgoing->type)
&& ((end_incoming-1)->type == (end_outgoing-1)->type)
&& std::all_of(
begin_outgoing + 1, end_outgoing - 1,
[](Sparticle const & p){ return p.type == pid::gluon; }
);
}
bool is_FKL(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> outgoing
){
assert(std::is_sorted(begin(incoming), end(incoming), pz_less{}));
assert(valid_outgoing(begin(outgoing), end(outgoing)));
return is_FKL(
begin(incoming), end(incoming),
begin(outgoing), end(outgoing)
);
}
bool has_2_jets(Event const & event){
return event.jets().size() >= 2;
}
/**
* \brief Checks whether event is unordered backwards
* @param ev Event
* @returns Is Event Unordered Backwards
*
* Checks there is more than 3 constuents in the final state
* Checks there is more than 3 jets
* Checks the most backwards parton is a gluon
* Checks the most forwards jet is not a gluon
* Checks the rest of the event is FKL
* Checks the second most backwards is not a different boson
* Checks the unordered gluon actually forms a jet
*/
bool is_unordered_backward(Event const & ev){
auto const & in = ev.incoming();
auto const & out = ev.outgoing();
assert(std::is_sorted(begin(in), end(in), pz_less{}));
assert(valid_outgoing(begin(out), end(out)));
if(out.size() < 3) return false;
if(ev.jets().size() < 3) return false;
if(in.front().type == pid::gluon) return false;
if(out.front().type != pid::gluon) return false;
// When skipping the unordered emission
// the remainder should be a regular FKL event,
// except that the (new) first outgoing particle must not be a A,W,Z,H.
const auto FKL_begin = next(begin(out));
if(is_AWZH_boson(*FKL_begin)) return false;
if(!is_FKL(in, {FKL_begin, end(out)})) return false;
// check that the unordered gluon forms an extra jet
const auto jets = sorted_by_rapidity(ev.jets());
const auto indices = ev.particle_jet_indices({jets.front()});
return indices[0] >= 0 && indices[1] == -1;
}
/**
* \brief Checks for a forward unordered gluon emission
* @param ev Event
* @returns Is the event a forward unordered emission
*
* \see is_unordered_backward
*/
bool is_unordered_forward(Event const & ev){
auto const & in = ev.incoming();
auto const & out = ev.outgoing();
assert(std::is_sorted(begin(in), end(in), pz_less{}));
assert(valid_outgoing(begin(out), end(out)));
if(out.size() < 3) return false;
if(ev.jets().size() < 3) return false;
if(in.back().type == pid::gluon) return false;
if(out.back().type != pid::gluon) return false;
// When skipping the unordered emission
// the remainder should be a regular FKL event,
// except that the (new) last outgoing particle must not be a A,W,Z,H.
const auto FKL_end = prev(end(out));
if(is_AWZH_boson(*prev(FKL_end))) return false;
if(!is_FKL(in, {begin(out), FKL_end})) return false;
// check that the unordered gluon forms an extra jet
const auto jets = sorted_by_rapidity(ev.jets());
const auto indices = ev.particle_jet_indices({jets.back()});
return indices.back() >= 0 && indices[indices.size()-2] == -1;
}
using event_type::EventType;
EventType classify(Event const & ev){
if(! final_state_ok(ev.outgoing())) return EventType::bad_final_state;
if(! has_2_jets(ev)) return EventType::no_2_jets;
if(is_FKL(ev.incoming(), ev.outgoing())) return EventType::FKL;
if(is_unordered_backward(ev)){
return EventType::unordered_backward;
}
if(is_unordered_forward(ev)){
return EventType::unordered_forward;
}
return EventType::nonFKL;
}
+ Sparticle extract_particle(LHEF::HEPEUP const & hepeup, int i){
+ return Sparticle{
+ static_cast<ParticleID>(hepeup.IDUP[i]),
+ fastjet::PseudoJet{
+ hepeup.PUP[i][0], hepeup.PUP[i][1],
+ hepeup.PUP[i][2], hepeup.PUP[i][3]
+ }
+ };
+ }
+
+ bool is_decay_product(std::pair<int, int> const & mothers){
+ if(mothers.first == 0) return false;
+ return mothers.second == 0 || mothers.first == mothers.second;
+ }
+
}
UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup):
central(EventParameters{
hepeup.scales.mur, hepeup.scales.muf, hepeup.weight()
})
{
size_t in_idx = 0;
for (int i = 0; i < hepeup.NUP; ++i) {
- // Only consider final state particles
- if (abs(hepeup.ISTUP[i]) != status_out) continue;
+ // skip decay products
+ // we will add them later on, but we have to ensure that
+ // the decayed particle is added before
+ if(is_decay_product(hepeup.MOTHUP[i])) continue;
- // Save particle information
- Sparticle particle{
- static_cast<ParticleID>(hepeup.IDUP[i]),
- fastjet::PseudoJet{
- hepeup.PUP[i][0], hepeup.PUP[i][1],
- hepeup.PUP[i][2], hepeup.PUP[i][3]
- }
- };
+ auto particle = extract_particle(hepeup, i);
+ // needed to identify mother particles for decay products
+ particle.p.set_user_index(i+1);
- if (hepeup.ISTUP[i]>0) outgoing.emplace_back(std::move(particle));
- else{
- if(in_idx > incoming.size()) {
- throw std::invalid_argument{
- "Event has too many incoming particles"
- };
- }
- incoming[in_idx++] = std::move(particle);
+ if(hepeup.ISTUP[i] == status_in){
+ if(in_idx > incoming.size()) {
+ throw std::invalid_argument{
+ "Event has too many incoming particles"
+ };
+ }
+ incoming[in_idx++] = std::move(particle);
}
+ else outgoing.emplace_back(std::move(particle));
}
std::sort(
begin(incoming), end(incoming),
[](Sparticle o1, Sparticle o2){return o1.p.pz()<o2.p.pz();}
);
std::sort(begin(outgoing), end(outgoing), rapidity_less{});
+
+ // add decay products
+ for (int i = 0; i < hepeup.NUP; ++i) {
+ if(!is_decay_product(hepeup.MOTHUP[i])) continue;
+ const int mother_id = hepeup.MOTHUP[i].first;
+ const auto mother = std::find_if(
+ begin(outgoing), end(outgoing),
+ [mother_id](Sparticle const & particle){
+ return particle.p.user_index() == mother_id;
+ }
+ );
+ if(mother == end(outgoing)){
+ throw std::invalid_argument{"invalid decay product parent"};
+ }
+ const int mother_idx = std::distance(begin(outgoing), mother);
+ assert(mother_idx > 0);
+ decays[mother_idx].emplace_back(extract_particle(hepeup, i));
+ }
}
Event::Event(
UnclusteredEvent ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
):
ev_{std::move(ev)},
cs_{to_PseudoJet(filter_partons(ev_.outgoing)), jet_def},
min_jet_pt_{min_jet_pt}
{
type_ = classify(*this);
}
std::vector<fastjet::PseudoJet> Event::jets() const{
return cs_.inclusive_jets(min_jet_pt_);
}
/**
* \brief Returns the invarient mass of the event
* @param ev Event
* @returns s hat
*
* Makes use of the FastJet PseudoJet function m2().
* Applies this function to the sum of the incoming partons.
*/
double shat(Event const & ev){
return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
}
namespace{
// colour flow according to Les Houches standard
// TODO: stub
std::vector<std::pair<int, int>> colour_flow(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing
){
std::vector<std::pair<int, int>> result(
incoming.size() + outgoing.size()
);
for(auto & col: result){
col = std::make_pair(-1, -1);
}
return result;
}
}
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){
LHEF::HEPEUP result;
result.heprup = heprup;
result.weights = {{event.central().weight, nullptr}};
for(auto const & var: event.variations()){
result.weights.emplace_back(var.weight, nullptr);
}
- const size_t num_particles =
- event.incoming().size() + event.outgoing().size();
+ size_t num_particles = event.incoming().size() + event.outgoing().size();
+ for(auto const & decay: event.decays()) num_particles += decay.second.size();
result.NUP = num_particles;
// the following entries are pretty much meaningless
result.IDPRUP = 1;
result.AQEDUP = 1./128.;
// end meaningless part
result.XWGTUP = event.central().weight;
result.SCALUP = event.central().muf;
result.scales.muf = event.central().muf;
result.scales.mur = event.central().mur;
result.scales.SCALUP = event.central().muf;
result.pdfinfo.p1 = event.incoming().front().type;
result.pdfinfo.p2 = event.incoming().back().type;
result.pdfinfo.scale = event.central().muf;
for(Sparticle const & in: event.incoming()){
result.IDUP.emplace_back(in.type);
result.ISTUP.emplace_back(status_in);
result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()});
result.MOTHUP.emplace_back(0, 0);
}
- for(Sparticle const & out: event.outgoing()){
+ for(size_t i = 0; i < event.outgoing().size(); ++i){
+ Sparticle const & out = event.outgoing()[i];
result.IDUP.emplace_back(out.type);
- result.ISTUP.emplace_back(status_out);
+ const int status = event.decays().count(i)?status_decayed:status_out;
+ result.ISTUP.emplace_back(status);
result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
result.MOTHUP.emplace_back(1, 2);
}
result.ICOLUP = colour_flow(
event.incoming(), filter_partons(event.outgoing())
);
if(result.ICOLUP.size() < num_particles){
const size_t AWZH_boson_idx = std::find_if(
begin(event.outgoing()), end(event.outgoing()),
[](Sparticle const & s){ return is_AWZH_boson(s); }
) - begin(event.outgoing()) + event.incoming().size();
assert(AWZH_boson_idx <= result.ICOLUP.size());
result.ICOLUP.insert(
begin(result.ICOLUP) + AWZH_boson_idx,
std::make_pair(0,0)
);
}
+ for(auto const & decay: event.decays()){
+ for(auto const out: decay.second){
+ result.IDUP.emplace_back(out.type);
+ result.ISTUP.emplace_back(status_out);
+ result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
+ const int mother_idx = 1 + event.incoming().size() + decay.first;
+ result.MOTHUP.emplace_back(mother_idx, mother_idx);
+ result.ICOLUP.emplace_back(0,0);
+ }
+ }
assert(result.ICOLUP.size() == num_particles);
static constexpr double unknown_spin = 9.; //per Les Houches accord
result.VTIMUP = std::vector<double>(num_particles, unknown_spin);
result.SPINUP = result.VTIMUP;
return result;
}
}
diff --git a/src/HepMCWriter.cc b/src/HepMCWriter.cc
index 25e848f..1402632 100644
--- a/src/HepMCWriter.cc
+++ b/src/HepMCWriter.cc
@@ -1,108 +1,113 @@
#include "RHEJ/HepMCWriter.hh"
#include <cassert>
#ifdef RHEJ_BUILD_WITH_HepMC
#include "HepMC/LHEFAttributes.h"
#include "HepMC/GenVertex.h"
#include "HepMC/GenParticle.h"
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
+#include "RHEJ/exceptions.hh"
+
namespace RHEJ{
namespace {
void add_generator_tag(LHEF::HEPRUP & heprup){
heprup.generators.emplace_back(LHEF::XMLTag{});
heprup.generators.back().name = "HEJ";
heprup.generators.back().version = "0.0.1";
}
void reset_weight_info(LHEF::HEPRUP & heprup){
heprup.IDWTUP = -4;
// 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::FourVector to_FourVector(Sparticle const & sp){
return {sp.px(), sp.py(), sp.pz(), sp.E()};
}
HepMC::shared_ptr<HepMC::GenRunInfo> & add_HEPRUP(
HepMC::shared_ptr<HepMC::GenRunInfo> & runinfo,
LHEF::HEPRUP heprup
){
add_generator_tag(heprup);
reset_weight_info(heprup);
auto hepr = HepMC::make_shared<HepMC::HEPRUPAttribute>();
hepr->heprup = std::move(heprup);
runinfo->add_attribute("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 runinfo;
}
}
HepMCWriter::HepMCWriter(std::string const & file, LHEF::HEPRUP heprup):
runinfo_{HepMC::make_shared<HepMC::GenRunInfo>()},
writer_{file, add_HEPRUP(runinfo_, std::move(heprup))}
{}
HepMCWriter::~HepMCWriter(){
writer_.close();
}
void HepMCWriter::write(Event const & ev){
+ if(!ev.decays().empty()){
+ throw not_implemented{"HepMC output for events with decays"};
+ }
auto vx = HepMC::make_shared<HepMC::GenVertex>();
for(auto const & in: ev.incoming()){
vx->add_particle_in(
HepMC::make_shared<HepMC::GenParticle>(
to_FourVector(in), static_cast<int>(in.type), -1
)
);
}
for(auto const & out: ev.outgoing()){
vx->add_particle_out(
HepMC::make_shared<HepMC::GenParticle>(
to_FourVector(out), static_cast<int>(out.type), +1
)
);
}
HepMC::GenEvent out_ev{runinfo_, HepMC::Units::GEV, HepMC::Units::MM};
out_ev.add_vertex(vx);
out_ev.weights().reserve(ev.variations().size() + 1u);
out_ev.weights().emplace_back(ev.central().weight);
for(auto const & var: ev.variations()){
out_ev.weights().emplace_back(var.weight);
}
writer_.write_event(out_ev);
}
}
#else
namespace RHEJ{
HepMCWriter::HepMCWriter(std::string const &, LHEF::HEPRUP){
throw std::invalid_argument(
"Failed to create HepMC writer: "
"Reversed HEJ was built without HepMC support"
);
}
void HepMCWriter::write(Event const &){
assert(false);
}
HepMCWriter::~HepMCWriter() = default;
}
#endif

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 4:41 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805206
Default Alt Text
(22 KB)

Event Timeline