Page MenuHomeHEPForge
Files F8723733

No OneTemporary
Actions

View Options

diff --git a/config.yml b/config.yml
index 26e46a7..be28642 100644
--- a/config.yml
+++ b/config.yml
@@ -1,94 +1,94 @@
# number of attempted resummation phase space points for each input event
trials: 10
min extparton pt: 30 # minimum transverse momentum of extremal partons
resummation jets: # resummation jet properties
min pt: 35 # minimum jet transverse momentum
algorithm: antikt # jet algorithm
R: 0.4 # jet R parameter
fixed order jets: # properties of input jets
min pt: 30
# by default, algorithm and R are like for resummation jets
# treatment of he various event classes
# the supported settings are: reweight, keep, discard
# non-FKL events cannot be reweighted
FKL: reweight
unordered: keep
-non-FKL: keep
+non-HEJ: keep
# scale settings similar to original HEJ
#
# Use combinations of max jet pperp, input scales, ht/2,
# and the jet invariant mass and vary all scales by factors
# of 1, sqrt(2), and 2. Discard combinations where mur and muf
# differ by a factor of more than two.
#
# The weight entries in the final events are ordered as follows:
# 0-18: max jet pperp
# 19-37: input scales
# 38-56: ht/2
# 57-75: jet invariant mass
# In each of these groups, the first entry corresponds to the basic
# scale choice. In the following entries, mur and muf are varied with
# the above factors. The entries are ordered lexicographically so that
# mur1 < mur2 or (mur1 == mur2 and muf1 < muf2).
#
# Note that in contrast to HEJ, the central choice for the event is always
# max jet pperp and cannot be configured (yet).
#
# scales: [max jet pperp, input, Ht/2, jet invariant mass]
# scale factors: [0.5, 0.7071, 1, 1.41421, 2]
# max scale ratio: 2.0001
scales: 91.188
# import scale setting functions
#
# import scales:
# lib_my_scales.so: [scale0,scale1]
log correction: false # whether or not to include higher order logs
unweight: false # TODO: whether or not to unweight events
# event output files
#
# the supported formats are
# - Les Houches (suffix .lhe)
# - HepMC (suffix .hepmc3)
# TODO: - ROOT ntuples (suffix .root)
#
# An output file's format is deduced either automatically from the suffix
# or from an explicit specification, e.g.
# - Les Houches: outfile
event output:
- RHEJ.lhe
# - RHEJ_events.hepmc
analysis:
# to use a custom analysis
# plugin: ./src/analysis-plugins/libVBF.so
# output: RHEJ.root
# wtwt cut: # optional cut on (event weight)^2
# to use rivet
# rivet: MC_XS # rivet analysis name
# output: RHEJ # name of the yoda files, ".yoda" and scale suffix will be added
# selection of random number generator and seed
# the choices are
# - mixmax (seed is an integer)
# - ranlux64 (seed is a filename containing parameters)
random generator:
name: mixmax
# seed: 1
# parameters for Higgs-gluon couplings
# this requires compilation with looptools
# Higgs coupling:
# use impact factors: false
# mt: 174
# include bottom: true
# mb: 4.7
diff --git a/doc/sphinx/rHEJ.rst b/doc/sphinx/rHEJ.rst
index 6caec96..6919dc4 100644
--- a/doc/sphinx/rHEJ.rst
+++ b/doc/sphinx/rHEJ.rst
@@ -1,282 +1,281 @@
.. _`Running reversed HEJ`:
Running reversed HEJ
====================
Quick start
-----------
In order to run reversed HEJ, you need a configuration file and a file
containing fixed-order events. A sample configuration is given by the
:file:`config.yml` file distributed together with reversed HEJ. Events
in the Les Houches Event File format can be generated with standard
Monte Carlo generators like `MadGraph5_aMC@NLO
<https://launchpad.net/mg5amcnlo>`_ or `Sherpa
<https://sherpa.hepforge.org/trac/wiki>`_. Reversed HEJ assumes that the
cross section is given by the sum of the event weights. Depending on the
fixed-order generator it may be necessary to adjust the weights in the
Les Houches Event File accordingly.
The processes supported by reversed HEJ are
- Pure multijet production
- Production of a Higgs boson with jets
..
- *TODO* Production of a W boson with jets
- *TODO* Production of a Z boson or photon with jets
where at least two jets are required in each case. For the time being,
only leading-order events are supported.
After generating an event file :file:`events.lhe` adjust the parameters
under the `fixed order jets`_ setting in :file:`config.yml` to the
settings in the fixed-order generation. Resummation can then be added by
running::
rHEJ config.yml events.lhe
Using the default settings, this will produce an output event file
:file:`RHEJ.lhe` with events including high-energy resummation.
.. _`reversed HEJ settings`:
Settings
--------
Reversed HEJ configuration files follow the `YAML <http://yaml.org/>`_
format. The following configuration parameters are supported:
.. _`trials`:
**trials**
High-energy resummation is performed by generating a number of
resummation phase space configurations corresponding to the input
fixed-order event. This parameter specifies how many such
configurations reversed HEJ should try to generate for each input
event. Typical values vary between 10 and 100.
.. _`min extparton pt`:
**min extparton pt**
Specifies the minimum transverse momentum in GeV of the most forward
and the most backward parton. This setting is needed to regulate an
otherwise uncancelled divergence. Its value should be slightly below
the minimum transverse momentum of jets specified by `resummation
jets: min pt`_. See also the `max ext soft pt fraction`_ setting.
.. _`max ext soft pt fraction`:
**max ext soft pt fraction**
Specifies the maximum fraction that soft radiation can contribute to
the transverse momentum of each the most forward and the most backward
jet. Values between around 0.05 and 0.1 are recommended. See also the
`min extparton pt`_ setting.
.. _`fixed order jets`:
**fixed order jets**
This tag collects a number of settings specifying the jet definition
in the event input. The settings should correspond to the ones used in
the fixed-order Monte Carlo that generated the input events.
.. _`fixed order jets: min pt`:
**min pt**
Minimum transverse momentum in GeV of fixed-order jets.
.. _`fixed order jets: algorithm`:
**algorithm**
The algorithm used to define jets. Allowed settings are
:code:`kt`, :code:`cambridge`, :code:`antikt`, :code:`genkt`,
:code:`cambridge for passive`, :code:`genkt for passive`,
:code:`ee kt`, :code:`ee genkt`. See the `FastJet
<http://fastjet.fr/>`_ documentation for a description of these
algorithms.
.. _`fixed order jets: R`:
**R**
The R parameter used in the jet algorithm, roughly corresponding
to the jet radius in the plane spanned by the rapidity and the
azimuthal angle.
.. _`resummation jets`:
**resummation jets**
This tag collects a number of settings specifying the jet definition
in the observed, i.e. resummed events. These settings are optional, by
default the same values as for the `fixed order jets`_ are assumed.
.. _`resummation jets: min pt`:
**min pt**
Minimum transverse momentum in GeV of resummation jets. This
should be around 25% larger than the minimum transverse momentum
of fixed order jets set by `fixed order jets: min pt`_.
.. _`resummation jets: algorithm`:
**algorithm**
The algorithm used to define jets. The reversed HEJ approach to
resummation relies on properties of :code:`antikt` jets, so this
value is strongly recommended. For a list of possible other
values, see the `fixed order jets: algorithm`_ setting.
.. _`resummation jets: R`:
**R**
The R parameter used in the jet algorithm.
.. _`FKL`:
**FKL**
Specifies how to treat FKL events. The possible values are
:code:`reweight` to enable resummation, :code:`keep` to keep the
events as they are up to a possible change of renormalisation and
factorisation scale, and :code:`discard` to discard these events.
.. _`unordered`:
**unordered**
Specifies how to treat events with one emission that does not respect
FKL ordering. The possible values are the same as for the `FKL`_
setting, but :code:`reweight` may not be supported for all process
types.
-.. TODO: rename to non-HEJ
-.. _`non-FKL`:
+.. _`non-HEJ`:
-**non-FKL**
+**non-HEJ**
Specifies how to treat events where no resummation is possible. The
allowed values are :code:`keep` to keep the events as they are up to
a possible change of renormalisation and factorisation scale and
:code:`discard` to discard these events.
.. _`scales`:
**scales**
Specifies the renormalisation and factorisation scales for the output
events. This can either be a single entry or a list :code:`[scale1,
scale2, ...]`. For the case of a list the first entry defines the
central scale. Possible values are fixed numbers to set the scale in
GeV or the following:
- :code:`H_T`: The sum of the scalar transverse momenta of all
final-state particles
- :code:`max jet pperp`: The maximum transverse momentum of all jets
- :code:`jet invariant mass`: Sum of the invariant masses of all jets
- :code:`m_j1j2`: Invariant mass between the two hardest jets.
Scales can be multiplied or divided by an overall factor,
e.g. :code:`H_T/2`.
It is also possible to import scales from an external library, see
:ref:`Custom scales`
.. _`scale factors`:
**scale factors**
A list of numeric factors by which each of the `scales`_ should be
multiplied. Renormalisation and factorisation scales are varied
independently. For example, a list with entries :code:`[0.5, 2]`
would give the four scale choices (0.5μ\ :sub:`r`, 0.5μ\ :sub:`f`);
(0.5μ\ :sub:`r`, 2μ\ :sub:`f`); (2μ\ :sub:`r`, 0.5μ\ :sub:`f`); (2μ\
:sub:`r`, 2μ\ :sub:`f`) in this order. The ordering corresponds to
the order of the final event weights.
.. _`max scale ratio`:
**max scale ratio**
Specifies the maximum factor by which renormalisation and
factorisation scales may difer. For a value of :code:`2` and the
example given for the `scale factors`_ the scale choices
(0.5μ\ :sub:`r`, 2μ\ :sub:`f`) and (2μ\ :sub:`r`, 0.5μ\ :sub:`f`)
will be discarded.
.. _`log correction`:
**log correction**
Whether to include corrections due to the evolution of the strong
coupling constant in the virtual corrections. Allowed values are
:code:`true` and :code:`false`.
.. TODO: unweight
.. _`event output`:
**event output**
Specifies the name of a single event output file or a list of such
files. The file format is either specified explicitly or derived from
the suffix. For example, :code:`events.lhe` or, equivalently
:code:`Les Houches: events.lhe` generates an output event file
:code:`events.lhe` in the Les Houches format. The supported formats
are
- :code:`file.lhe` or :code:`Les Houches: file`: The Les Houches
event file format.
- :code:`file.hepmc` or :code:`HepMC: file`: The HepMC format.
.. _`random generator`:
**random generator**
Sets parameters for random number generation.
.. _`random generator: name`:
**name**
Which random number generator to use. Currently, :code:`mixmax`
and :code:`ranlux64` are implemented. Mixmax is recommended. See
the `CLHEP documentation
<http://proj-clhep.web.cern.ch/proj-clhep/index.html#docu>`_ for
details on the generators.
.. _`random generator: seed`:
**seed**
The seed for random generation. This should be a single number for
mixmax and the name of a state file for ranlux64.
.. _`analysis`:
**analysis**
Name and Setting for the event analyses; either a custom
analysis plugin or rivet. For the first the :code:`plugin` sub-entry
should be set to the analysis file path. All further entries are passed on
to the analysis. To use rivet a list of rivet-analysis have to be
given in :code:`rivet` and prefix for the yoda file has to be set
through :code:`output`. See :ref:`Writing custom analyses` for details.
.. _`Higgs coupling`:
**Higgs coupling**
This collects a number of settings concerning the effective coupling
of the Higgs boson to gluons. This is only relevant for the
production process of a Higgs boson with jets and only supported if
reversed HEJ was compiled with `QCDloop
<https://github.com/scarrazza/qcdloop>`_ support.
.. _`Higgs coupling: use impact factors`:
**use impact factors**
Whether to use impact factors for the coupling to the most forward
and most backward partons. Impact factors imply the infinite
top-quark mass limit.
.. _`Higgs coupling: mt`:
**mt**
The value of the top-quark mass in GeV. If this is not specified,
the limit of an infinite mass is taken.
.. _`Higgs coupling: include bottom`:
**include bottom**
Whether to include the Higgs coupling to bottom quarks.
.. _`Higgs coupling: mb`:
**mb**
The value of the bottom-quark mass in GeV.
diff --git a/include/RHEJ/event_types.hh b/include/RHEJ/event_types.hh
index b3a2331..b7e8a42 100644
--- a/include/RHEJ/event_types.hh
+++ b/include/RHEJ/event_types.hh
@@ -1,44 +1,44 @@
/** \file event_types.hh
* \brief This file details the classification of events.
*
* This file makes use of a macro in order to avoid repetition of EventType names.
* To this end, the Documentation of the EventType enumeration below is incomplete.
* The possible event classifications are stored within that enumeration.
*/
#pragma once
#include "RHEJ/utility.hh"
namespace RHEJ{
// macro definition to avoid repetition of EventType names
-#define RHEJ_EVENT_TYPES(F) F(FKL), F(unordered_backward), F(unordered_forward), F(nonFKL), F(no_2_jets), F(bad_final_state)
+#define RHEJ_EVENT_TYPES(F) F(FKL), F(unordered_backward), F(unordered_forward), F(nonHEJ), F(no_2_jets), F(bad_final_state)
#define RHEJ_AS_ENUM(VAR) VAR
#define RHEJ_AS_STRING(VAR) #VAR
//! Event_Type NameSpace
namespace event_type{
/** \enum EventType
* \brief EventType enumeration gives different possible event types
*
* This enumeration is used to distinguish between different event types.
*/
enum EventType: size_t{
RHEJ_EVENT_TYPES(RHEJ_AS_ENUM), /**< Macro of Possible States */
unob = unordered_backward, /**< Unordered Emission by backwards Jet */
unof = unordered_forward, /**< Unordered Emission by forwards Jet */
first_type = FKL, /**< FKL Event */
last_type = bad_final_state /**< Bad Final State Event */
};
static constexpr auto names = make_array(
RHEJ_EVENT_TYPES(RHEJ_AS_STRING)
);
}
#undef RHEJ_HISTOGRAMS
#undef RHEJ_AS_ENUM
#undef RHEJ_AS_STRING
}
diff --git a/src/Event.cc b/src/Event.cc
index fb6a40c..0886a63 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,343 +1,343 @@
#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;
+ return EventType::nonHEJ;
}
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) {
// 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;
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] == 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);
}
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 = event.type()+1; // event ID
result.AQEDUP = 1./128.; // alpha_EW
//result.AQCDUP = 0.118 // alpha_QCD
// 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(size_t i = 0; i < event.outgoing().size(); ++i){
Sparticle const & out = event.outgoing()[i];
result.IDUP.emplace_back(out.type);
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/YAMLreader.cc b/src/YAMLreader.cc
index 2ed2833..a3e0108 100644
--- a/src/YAMLreader.cc
+++ b/src/YAMLreader.cc
@@ -1,491 +1,491 @@
#include "RHEJ/YAMLreader.hh"
#include <set>
#include <string>
#include <vector>
#include <iostream>
#include <stdexcept>
#include <dlfcn.h>
namespace RHEJ{
namespace{
//! Get YAML tree of supported options
/**
* The configuration file is checked against this tree of options
* in assert_all_options_known.
*/
YAML::Node const & get_supported_options(){
const static YAML::Node supported = [](){
YAML::Node supported;
static const auto opts = {
"trials", "min extparton pt", "max ext soft pt fraction",
- "FKL", "unordered", "non-FKL",
+ "FKL", "unordered", "non-HEJ",
"scales", "scale factors", "max scale ratio", "import scales",
"log correction", "unweight", "event output", "analysis"
};
// add subnodes to "supported" - the assigned value is irrelevant
for(auto && opt: opts) supported[opt] = "";
for(auto && jet_opt: {"min pt", "algorithm", "R"}){
supported["resummation jets"][jet_opt] = "";
supported["fixed order jets"][jet_opt] = "";
}
for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
supported["Higgs coupling"][opt] = "";
}
for(auto && opt: {"name", "seed"}){
supported["random generator"][opt] = "";
}
return supported;
}();
return supported;
}
fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){
using namespace fastjet;
static const std::map<std::string, fastjet::JetAlgorithm> known = {
{"kt", kt_algorithm},
{"cambridge", cambridge_algorithm},
{"antikt", antikt_algorithm},
{"genkt", genkt_algorithm},
{"cambridge for passive", cambridge_for_passive_algorithm},
{"genkt for passive", genkt_for_passive_algorithm},
{"ee kt", ee_kt_algorithm},
{"ee genkt", ee_genkt_algorithm},
{"plugin", plugin_algorithm}
};
const auto res = known.find(algo);
if(res == known.end()){
throw std::invalid_argument("Unknown jet algorithm " + algo);
}
return res->second;
}
EventTreatment to_EventTreatment(std::string const & name){
static const std::map<std::string, EventTreatment> known = {
{"reweight", EventTreatment::reweight},
{"keep", EventTreatment::keep},
{"discard", EventTreatment::discard}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown event treatment " + name);
}
return res->second;
}
} // namespace anonymous
ParticleID to_ParticleID(std::string const & name){
using namespace RHEJ::pid;
static const std::map<std::string, ParticleID> known = {
{"d", d}, {"down", down}, {"u", u}, {"up", up}, {"s", s}, {"strange", strange},
{"c", c}, {"charm", charm}, {"b", b}, {"bottom", bottom}, {"t", t}, {"top", top},
{"e", e}, {"electron", electron}, {"nu_e", nu_e}, {"electron_neutrino", electron_neutrino},
{"mu", mu}, {"muon", muon}, {"nu_mu", nu_mu}, {"muon_neutrino", muon_neutrino},
{"tau", tau}, {"nu_tau", nu_tau}, {"tau_neutrino", tau_neutrino},
{"d_bar", d_bar}, {"u_bar", u_bar}, {"s_bar", s_bar}, {"c_bar", c_bar},
{"b_bar", b_bar}, {"t_bar", t_bar}, {"e_bar", e_bar},
{"nu_e_bar", nu_e_bar}, {"mu_bar", mu_bar}, {"nu_mu_bar", nu_mu_bar},
{"tau_bar", tau_bar}, {"nu_tau_bar", nu_tau_bar},
{"gluon", gluon}, {"g", g}, {"photon", photon}, {"gamma", gamma},
{"Z", Z}, {"Wp", Wp}, {"Wm", Wm}, {"W+", Wp}, {"W-", Wm},
{"h", h}, {"H", h}, {"Higgs", Higgs}, {"higgs", higgs},
{"p", p}, {"proton", proton}, {"p_bar", p_bar}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown particle " + name);
}
return res->second;
}
namespace detail{
void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml){
setting = to_JetAlgorithm(yaml.as<std::string>());
}
void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml){
setting = to_EventTreatment(yaml.as<std::string>());
}
void set_from_yaml(ParticleID & setting, YAML::Node const & yaml){
setting = to_ParticleID(yaml.as<std::string>());
}
} // namespace detail
JetParameters get_jet_parameters(
YAML::Node const & node,
std::string const & entry
){
assert(node);
JetParameters result;
fastjet::JetAlgorithm jet_algo = fastjet::antikt_algorithm;
double R;
set_from_yaml_if_defined(jet_algo, node, entry, "algorithm");
set_from_yaml(R, node, entry, "R");
result.def = fastjet::JetDefinition{jet_algo, R};
set_from_yaml(result.min_pt, node, entry, "min pt");
return result;
}
RNGConfig to_RNGConfig(
YAML::Node const & node,
std::string const & entry
){
assert(node);
RNGConfig result;
set_from_yaml(result.name, node, entry, "name");
set_from_yaml_if_defined(result.seed, node, entry, "seed");
return result;
}
HiggsCouplingSettings get_Higgs_coupling(
YAML::Node const & node,
std::string const & entry
){
assert(node);
static constexpr double mt_max = 2e4;
#ifndef RHEJ_BUILD_WITH_QCDLOOP
if(node[entry]){
throw std::invalid_argument{
"Higgs coupling settings require building Reversed HEJ "
"with QCDloop support"
};
}
#endif
HiggsCouplingSettings settings;
set_from_yaml_if_defined(settings.mt, node, entry, "mt");
set_from_yaml_if_defined(settings.mb, node, entry, "mb");
set_from_yaml_if_defined(settings.include_bottom, node, entry, "include bottom");
set_from_yaml_if_defined(settings.use_impact_factors, node, entry, "use impact factors");
if(settings.use_impact_factors){
if(settings.mt != std::numeric_limits<double>::infinity()){
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
else{
// huge values of the top mass are numerically unstable
settings.mt = std::min(settings.mt, mt_max);
}
return settings;
}
FileFormat to_FileFormat(std::string const & name){
static const std::map<std::string, FileFormat> known = {
{"Les Houches", FileFormat::Les_Houches},
{"HepMC", FileFormat::HepMC}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown file format " + name);
}
return res->second;
}
std::string extract_suffix(std::string const & filename){
size_t separator = filename.rfind('.');
if(separator == filename.npos) return {};
return filename.substr(separator + 1);
}
FileFormat format_from_suffix(std::string const & filename){
const std::string suffix = extract_suffix(filename);
if(suffix == "lhe") return FileFormat::Les_Houches;
if(suffix == "hepmc") return FileFormat::HepMC;
throw std::invalid_argument{
"Can't determine format for output file " + filename
};
}
void assert_all_options_known(
YAML::Node const & conf, YAML::Node const & supported
){
if(!conf.IsMap()) return;
if(!supported.IsMap()) throw invalid_type{"must not have sub-entries"};
for(auto const & entry: conf){
const auto name = entry.first.as<std::string>();
if(! supported[name]) throw unknown_option{name};
/* check sub-options, e.g. 'resummation jets: min pt'
* we don't check analysis sub-options
* those depend on the analysis being used and should be checked there
* similar for "import scales"
*/
if(name != "analysis" && name != "import scales"){
try{
assert_all_options_known(conf[name], supported[name]);
}
catch(unknown_option const & ex){
throw unknown_option{name + ": " + ex.what()};
}
catch(invalid_type const & ex){
throw invalid_type{name + ": " + ex.what()};
}
}
}
}
} // namespace RHEJ
namespace YAML {
Node convert<RHEJ::OutputFile>::encode(RHEJ::OutputFile const & outfile) {
Node node;
node[to_string(outfile.format)] = outfile.name;
return node;
};
bool convert<RHEJ::OutputFile>::decode(Node const & node, RHEJ::OutputFile & out) {
switch(node.Type()){
case NodeType::Map: {
YAML::const_iterator it = node.begin();
out.format = RHEJ::to_FileFormat(it->first.as<std::string>());
out.name = it->second.as<std::string>();
return true;
}
case NodeType::Scalar:
out.name = node.as<std::string>();
out.format = RHEJ::format_from_suffix(out.name);
return true;
default:
return false;
}
}
} // namespace YAML
namespace RHEJ{
namespace detail{
void set_from_yaml(OutputFile & setting, YAML::Node const & yaml){
setting = yaml.as<OutputFile>();
}
}
namespace{
void update_fixed_order_jet_parameters(
JetParameters & fixed_order_jets, YAML::Node const & yaml
){
if(!yaml["fixed order jets"]) return;
set_from_yaml_if_defined(
fixed_order_jets.min_pt, yaml, "fixed order jets", "min pt"
);
fastjet::JetAlgorithm algo = fixed_order_jets.def.jet_algorithm();
set_from_yaml_if_defined(algo, yaml, "fixed order jets", "algorithm");
double R = fixed_order_jets.def.R();
set_from_yaml_if_defined(R, yaml, "fixed order jets", "R");
fixed_order_jets.def = fastjet::JetDefinition{algo, R};
}
// like std::stod, but throw if not the whole string can be converted
double to_double(std::string const & str){
std::size_t pos;
const double result = std::stod(str, &pos);
if(pos < str.size()){
throw std::invalid_argument(str + " is not a valid double value");
}
return result;
}
void import_scale_functions(
std::string const & file,
std::vector<std::string> const & scale_names,
std::unordered_map<std::string, ScaleFunction> & known
) {
using ScaleFunction = double (*)(Event const &);
auto handle = dlopen(file.c_str(), RTLD_NOW);
char * error = dlerror();
if(error != nullptr) throw std::runtime_error{error};
for(auto const & scale: scale_names) {
void * sym = dlsym(handle, scale.c_str());
error = dlerror();
if(error != nullptr) throw std::runtime_error{error};
known.emplace(scale, reinterpret_cast<ScaleFunction>(sym));
}
}
auto get_scale_map(
YAML::Node const & yaml
) {
std::unordered_map<std::string, ScaleFunction> scale_map;
scale_map.emplace("H_T", H_T);
scale_map.emplace("max jet pperp", max_jet_pt);
scale_map.emplace("jet invariant mass", jet_invariant_mass);
scale_map.emplace("m_j1j2", m_j1j2);
if(yaml["import scales"]) {
if(! yaml["import scales"].IsMap()) {
throw invalid_type{"Entry 'import scales' is not a map"};
}
for(auto const & import: yaml["import scales"]) {
const auto file = import.first.as<std::string>();
const auto scale_names =
import.second.IsSequence()
?import.second.as<std::vector<std::string>>()
:std::vector<std::string>{import.second.as<std::string>()};
import_scale_functions(file, scale_names, scale_map);
}
}
return scale_map;
}
// simple (as in non-composite) scale functions
/**
* An example for a simple scale function would be H_T,
* H_T/2 is then composite (take H_T and then divide by 2)
*/
ScaleFunction parse_simple_ScaleFunction(
std::string const & scale_fun,
std::unordered_map<std::string, ScaleFunction> const & known
) {
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
const auto it = known.find(scale_fun);
if(it != end(known)) return it->second;
try{
const double scale = to_double(scale_fun);
return FixedScale{scale};
} catch(std::invalid_argument const &){}
throw std::invalid_argument{"Unknown scale choice: " + scale_fun};
}
std::string trim_front(std::string const & str){
const auto new_begin = std::find_if(
begin(str), end(str), [](char c){ return ! std::isspace(c); }
);
return std::string(new_begin, end(str));
}
std::string trim_back(std::string str){
size_t pos = str.size() - 1;
// use guaranteed wrap-around behaviour to check whether we have
// traversed the whole string
for(; pos < str.size() && std::isspace(str[pos]); --pos) {}
str.resize(pos + 1); // note that pos + 1 can be 0
return str;
}
ScaleFunction parse_ScaleFunction(
std::string const & scale_fun,
std::unordered_map<std::string, ScaleFunction> const & known
){
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
const size_t delim = scale_fun.find_first_of("*/");
if(delim == scale_fun.npos){
return parse_simple_ScaleFunction(scale_fun, known);
}
const std::string first = trim_back(std::string{scale_fun, 0, delim});
const std::string second = trim_front(std::string{scale_fun, delim+1});
double factor;
ScaleFunction fun;
if(scale_fun[delim] == '/'){
factor = 1/to_double(second);
fun = parse_simple_ScaleFunction(first, known);
}
else{
assert(scale_fun[delim] == '*');
try{
factor = to_double(second);
fun = parse_simple_ScaleFunction(first, known);
}
catch(std::invalid_argument const &){
factor = to_double(first);
fun = parse_simple_ScaleFunction(second, known);
}
}
assert(fun != nullptr);
return Product{factor, std::move(fun)};
}
EventTreatMap get_event_treatment(
YAML::Node const & yaml
){
using namespace event_type;
EventTreatMap treat {
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{FKL, EventTreatment::reweight},
{unob, EventTreatment::keep},
{unof, EventTreatment::keep},
- {nonFKL, EventTreatment::keep}
+ {nonHEJ, EventTreatment::keep}
};
set_from_yaml(treat.at(FKL), yaml, "FKL");
set_from_yaml(treat.at(unob), yaml, "unordered");
treat.at(unof) = treat.at(unob);
- set_from_yaml(treat.at(nonFKL), yaml, "non-FKL");
- if(treat[nonFKL] == EventTreatment::reweight){
- throw std::invalid_argument{"Cannot reweight non-FKL events"};
+ set_from_yaml(treat.at(nonHEJ), yaml, "non-HEJ");
+ if(treat[nonHEJ] == EventTreatment::reweight){
+ throw std::invalid_argument{"Cannot reweight non-HEJ events"};
}
return treat;
}
Config to_Config(YAML::Node const & yaml){
try{
assert_all_options_known(yaml, get_supported_options());
}
catch(unknown_option const & ex){
throw unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
}
Config config;
config.resummation_jets = get_jet_parameters(yaml, "resummation jets");
config.fixed_order_jets = config.resummation_jets;
update_fixed_order_jet_parameters(config.fixed_order_jets, yaml);
set_from_yaml(config.min_extparton_pt, yaml, "min extparton pt");
config.max_ext_soft_pt_fraction = std::numeric_limits<double>::infinity();
set_from_yaml_if_defined(
config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction"
);
set_from_yaml(config.trials, yaml, "trials");
set_from_yaml(config.log_correction, yaml, "log correction");
set_from_yaml(config.unweight, yaml, "unweight");
config.treat = get_event_treatment(yaml);
set_from_yaml_if_defined(config.output, yaml, "event output");
config.rng = to_RNGConfig(yaml, "random generator");
set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
config.scales = to_ScaleConfig(yaml);
config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
return config;
}
} // namespace anonymous
ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
ScaleConfig config;
auto scale_funs = get_scale_map(yaml);
std::vector<std::string> scales;
set_from_yaml(scales, yaml, "scales");
config.base.reserve(scales.size());
std::transform(
begin(scales), end(scales), std::back_inserter(config.base),
[scale_funs](auto const & entry){
return parse_ScaleFunction(entry, scale_funs);
}
);
set_from_yaml_if_defined(config.factors, yaml, "scale factors");
config.max_ratio = std::numeric_limits<double>::infinity();
set_from_yaml_if_defined(config.max_ratio, yaml, "max scale ratio");
return config;
}
Config load_config(std::string const & config_file){
try{
return to_Config(YAML::LoadFile(config_file));
}
catch(...){
std::cerr << "Error reading " << config_file << ":\n ";
throw;
}
}
} // namespace RHEJ
diff --git a/t/check_res.cc b/t/check_res.cc
index a68442d..00bc897 100644
--- a/t/check_res.cc
+++ b/t/check_res.cc
@@ -1,88 +1,88 @@
#include <iostream>
#include "LHEF/LHEF.h"
#include "RHEJ/stream.hh"
#include "RHEJ/EventReweighter.hh"
#include "RHEJ/Ranlux64.hh"
namespace{
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
const fastjet::JetDefinition Born_jet_def{jet_def};
constexpr double Born_jetptmin = 30;
constexpr double extpartonptmin = 30;
constexpr double max_ext_soft_pt_fraction =
std::numeric_limits<double>::infinity();
constexpr double jetptmin = 35;
constexpr bool log_corr = false;
using EventTreatment = RHEJ::EventTreatment;
using namespace RHEJ::event_type;
RHEJ::EventTreatMap treat{
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
- {nonFKL, EventTreatment::discard},
+ {nonHEJ, EventTreatment::discard},
{unof, EventTreatment::discard},
{unob, EventTreatment::discard},
{FKL, EventTreatment::reweight}
};
};
int main(int argn, char** argv) {
if(argn == 5 && std::string(argv[4]) == "uno"){
--argn;
treat[unof] = EventTreatment::reweight;
treat[unob] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
if(argn != 4){
std::cerr << "Usage: check_res eventfile xsection tolerance [uno]";
return EXIT_FAILURE;
}
const double xsec_ref = std::stod(argv[2]);
const double tolerance = std::stod(argv[3]);
RHEJ::istream in{argv[1]};
LHEF::Reader reader{in};
RHEJ::PhaseSpacePointConfig psp_conf;
psp_conf.jet_param = RHEJ::JetParameters{jet_def, jetptmin};
psp_conf.min_extparton_pt = extpartonptmin;
psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;
RHEJ::MatrixElementConfig ME_conf;
ME_conf.jet_param = psp_conf.jet_param;
ME_conf.log_correction = log_corr;
ME_conf.Higgs_coupling = RHEJ::HiggsCouplingSettings{};
RHEJ::EventReweighterConfig conf;
conf.psp_config = std::move(psp_conf);
conf.ME_config = std::move(ME_conf);
conf.jet_param = psp_conf.jet_param;
conf.treat = treat;
reader.readEvent();
const bool has_Higgs = std::find(
begin(reader.hepeup.IDUP),
end(reader.hepeup.IDUP),
25
) != end(reader.hepeup.IDUP);
const double mu = has_Higgs?125.:91.188;
RHEJ::ScaleGenerator scale_gen{{RHEJ::FixedScale{mu}}, {}, 1.};
RHEJ::Ranlux64 ran{};
RHEJ::EventReweighter rhej{reader.heprup, std::move(scale_gen), conf, ran};
double xsec = 0.;
do{
RHEJ::Event ev{
RHEJ::UnclusteredEvent{reader.hepeup},
Born_jet_def, Born_jetptmin
};
auto resummed_events = rhej.reweight(ev, 10);
for(auto const & ev: resummed_events) xsec += ev.central().weight;
} while(reader.readEvent());
if(std::abs(xsec - xsec_ref) > tolerance){
std::cerr << "cross section is off: "
<< xsec << " != " << xsec_ref
<< " +- " << tolerance << '\n';
return EXIT_FAILURE;
}
}
diff --git a/t/jet_config.yml b/t/jet_config.yml
index 3d19554..ced0d0c 100644
--- a/t/jet_config.yml
+++ b/t/jet_config.yml
@@ -1,26 +1,26 @@
trials: 10
min extparton pt: 30
resummation jets:
min pt: 35
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
FKL: reweight
unordered: discard
-non-FKL: discard
+non-HEJ: discard
log correction: false
unweight: false
scales: 91.188
random generator:
name: ranlux64
event output:
- tst.lhe
diff --git a/t/jet_config_with_import.yml b/t/jet_config_with_import.yml
index 93e0ffa..467f765 100644
--- a/t/jet_config_with_import.yml
+++ b/t/jet_config_with_import.yml
@@ -1,29 +1,29 @@
trials: 10
min extparton pt: 30
resummation jets:
min pt: 35
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
FKL: reweight
unordered: discard
-non-FKL: discard
+non-HEJ: discard
log correction: false
unweight: false
scales: softest_jet_pt
event output:
- tst.lhe
random generator:
name: ranlux64
import scales:
./libscales.so: softest_jet_pt

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:14 PM (21 h, 57 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242434
Default Alt Text
(49 KB)

Event Timeline