Page MenuHomeHEPForge
Files F8725707

No OneTemporary
Actions

View Options

diff --git a/include/RHEJ/EventReweighter.hh b/include/RHEJ/EventReweighter.hh
index 6852cc7..c5ddf24 100644
--- a/include/RHEJ/EventReweighter.hh
+++ b/include/RHEJ/EventReweighter.hh
@@ -1,194 +1,194 @@
/** \file EventReweighter.hh
* \brief Main class for reweighting events according to the reversed HEJ method
*
* EventReweighter Class is the main class used within RHEJ. It reweights the
* resummation events.
*/
#pragma once
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "LHEF/LHEF.h"
#include "RHEJ/PDF.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
-#include "RHEJ/event_classes.hh"
+#include "RHEJ/event_types.hh"
#include "RHEJ/config.hh"
#include "RHEJ/PhaseSpacePoint.hh"
#include "RHEJ/MatrixElement.hh"
namespace RHEJ{
class ScaleGenerator;
/** \struct Beam EventReweighter.hh "include/RHEJ/EventReweighter.hh"
* \brief Beam Struct. Contains beam energy and incoming particle content.
*
* Contains the value of the energy of the beam and the identity of the two incoming Particles
*/
struct Beam{
double E;
std::array<ParticleID, 2> type;
};
/** \class EventReweighter EventReweighter.hh "include/RHEJ/EventReweighter.hh"
* \brief EventReweighter Main class for reweighting events in RHEJ.
*
* This is the class which reweights the resummation phase space points once they have been
* generated by PhaseSpacePoint. This class also handles how many phase space points are generated
* for every fixed order event which is processed.
*/
class EventReweighter{
- using EventClass = event_class::EventClass;
+ using EventType = event_type::EventType;
public:
EventReweighter(
Beam beam, /**< Beam Energy */
int pdf_id, /**< PDF ID */
fastjet::JetDefinition jet_def, /**< Jet Definition */
double jetptmin, /**< Jet Pt minimum */
EventTreatMap treat, /**< Treat??? */
double extpartonptmin, /**< External Parton Pt Minimum */
double max_ext_soft_pt_fraction, /**< External Parton Pt Maximum */
bool log_corr, /**< Log Correct On or Off */
HiggsCouplingSettings Higgs_coupling
);
EventReweighter(
LHEF::HEPRUP const & heprup, /**< HEPRUP???? */
fastjet::JetDefinition jet_def, /**< Jet Definition */
double jetptmin, /**< Jet Pt minimum */
EventTreatMap treat, /**< Treat??? */
double extpartonptmin, /**< External Parton Pt Minimum */
double max_ext_soft_pt_fraction, /**< External Parton Pt Maximum */
bool log_corr, /**< Log Correct On or Off */
HiggsCouplingSettings Higgs_coupling
);
PDF const & pdf() const; /**< PDF Used */
std::vector<Event> reweight(
Event const & ev, /**< Event For Reweighing */
int num_events /**< Number of Events */
);
std::vector<Event> reweight(
Event const & ev, /**< Event For Reweighing */
int num_events, /**< Number of Events */
ScaleGenerator const & gen_scales /**< Generated Scales see config.hh */
);
- //! The EventClass for the Last Event
- EventClass last_event_class() const;
+ //! The EventType for the Last Event
+ EventType last_event_type() const;
private:
/** \struct EventFactors EventReweighter.hh "include/RHEJ/EventReweighter.hh"
* \brief EventFactors Struct containing central value and variations
*
* Contains the Central Value and Variation around
*/
struct EventFactors{
double central; /**< Central Value for the Factor */
std::vector<double> variations; /**< Vector of Variations from central value */
};
template<typename... T>
PDF const & pdf(T&& ...);
std::vector<Event> gen_res_events(
Event const & ev, int num_events
);
std::vector<Event> rescale(
Event const & Born_ev, std::vector<Event> events
) const;
/**
* \brief Do the Jets pass the resummation Cuts?
*
* @param ev Event in Question
* @returns 0 or 1 depending on if ev passes Jet Cuts
*/
bool jets_pass_resummation_cuts(Event const & ev) const;
/**
* \brief pdf_factors Function
*
* @param ev Event in Question
* @returns EventFactor due to PDFs
*
* Calculates the Central value and the variation due
* to the PDF choice made.
*/
EventFactors pdf_factors(Event const & ev) const;
/**
* \brief matrix_elements Function
*
* @param ev Event in question
* @returns EventFactor due to MatrixElements
*
* Calculates the Central value and the variation due
* to the Matrix Element.
*/
EventFactors matrix_elements(Event const & ev) const;
/**
* \brief Scale-dependent part of fixed-order matrix element
*
* @param ev Event in question
* @returns EventFactor scale variation due to FO-ME.
*
* This is only called to compute the scale variation for events where
* we don't do resummation (e.g. non-FKL).
* Since at tree level the scale dependence is just due to alpha_s,
* it is enough to return the alpha_s(mur) factors in the matrix element.
* The rest drops out in the ratio of (output event ME)/(input event ME),
* so we never have to compute it.
*/
EventFactors fixed_order_scale_ME(Event const & ev) const;
/**
* \brief Computes the tree level matrix element
*
* @param ev Event in Question
* @returns HEJ approximation to Tree level Matrix Element
*
* This computes the HEJ approximation to the tree level FO
* Matrix element which is used within the LO weighting process.
*/
double tree_matrix_element(Event const & ev) const;
/**
* \brief Generation of Resummation Phase Space Points
*
* @param ev Event In Question
* @returns Resummation PhaseSpacePoint
*
* Where is this even called? I can't find it with a grep -r in top directory???
*/
PhaseSpacePoint gen_PhaseSpacePoint(Event const & ev) const;
- EventClass last_event_class_; /**< Class of the last event */
+ EventType last_event_type_; /**< Class of the last event */
double extpartonptmin_; /**< External Parton Minimum Pt */
double max_ext_soft_pt_fraction_; /**< External Parton Maximum Pt fraction */
double jetptmin_; /**< Jet Pt Minimum threshold */
double E_beam_; /**< Energy of Beam */
fastjet::JetDefinition jet_def_; /**< Jet Definition being used */
PDF pdf_; /**< Relevant PDF */
MatrixElement MEt2_; /**< Matrix Element Object */
EventTreatMap treat_;
};
template<typename... T>
PDF const & EventReweighter::pdf(T&&... t){
return pdf_ = PDF{std::forward<T>(t)...};
}
}
diff --git a/include/RHEJ/config.hh b/include/RHEJ/config.hh
index 5295d5c..6f506d2 100644
--- a/include/RHEJ/config.hh
+++ b/include/RHEJ/config.hh
@@ -1,343 +1,343 @@
/** \file config.hh
* \brief The file which handles the configuration file parameters
*
* Contains the JetParameters Struct and ScaleConfig Struct. Also
* contains the ScaleGenerator Class and the EventTreatment class.
* Config struct is also defined here. The configuration files
* parameters are read and then stored within this objects.
*/
#pragma once
#include <string>
#include <memory>
#include "fastjet/JetDefinition.hh"
#include "RHEJ/ScaleFunction.hh"
#include "RHEJ/optional.hh"
#include "RHEJ/output_formats.hh"
-#include "RHEJ/event_classes.hh"
+#include "RHEJ/event_types.hh"
#include "RHEJ/HiggsCouplingSettings.hh"
#include "RHEJ/exceptions.hh"
#include "RHEJ/utility.hh"
#include "yaml-cpp/yaml.h"
namespace RHEJ{
/**! \struct JetParameters config.hh "include/RHEJ/config.hh"
* \brief Contains Jet Definition and Minimum Pt to be a Jet.
*
* Contains the Fastjet definition of a Jet and a threshold (minimum) value for pt
*/
struct JetParameters{
fastjet::JetDefinition def; /**< Jet Definition */
double min_pt; /**< Minimum Jet Pt */
};
/**! \struct ScaleConfig config.hh "include/RHEJ/config.hh"
* \brief Sets up the possible choices for scale variation?
*
* There is a vector which contains the possible scale choices and
* also vector of doubles with the scale factors along with a max
* scale ratio.
*/
struct ScaleConfig{
std::vector<std::unique_ptr<ScaleFunction>> scales; /**< Vector of possible Scale choices */
std::vector<double> scale_factors; /**< Vector of possible Scale Factors */
double max_scale_ratio; /**< Maximum ratio for the scales */
};
/**! \class ScaleGenerator config.hh "include/RHEJ/config.hh"
* \brief A Class used to generate scales
*
* This class has a clustered event class and scale config struct
* within it.
*/
class ScaleGenerator{
public:
ScaleGenerator() = default;
explicit ScaleGenerator(ScaleConfig && settings):
settings_{std::move(settings)}
{}
Event operator()(Event ev) const;
ScaleConfig const & settings() const;
private:
ScaleConfig settings_;
};
/**
* \brief Enumeration which deals with how to treat events.
*
* The program will decide on how to treat an event based on
* the value of this enumeration.
*/
enum class EventTreatment{
reweight, /**< Reweigh the event */
keep, /**< Keep the event */
discard, /**< Discard the event */
};
/**
- * \brief An ordered Map with EventClass keys and mapped values EventTreatment
+ * \brief An ordered Map with EventType keys and mapped values EventTreatment
*
- * If called upon with EventTreatMap.at(EventClass) it will return the treatment which has
- * been specified for that EventClass.
+ * If called upon with EventTreatMap.at(EventType) it will return the treatment which has
+ * been specified for that EventType.
*/
- using EventTreatMap = std::map<event_class::EventClass, EventTreatment>;
+ using EventTreatMap = std::map<event_type::EventType, EventTreatment>;
/**! \struct Config config.hh "include/RHEJ/config.hh"
* \brief Config Struct for user input parameters.
*
* The struct which handles the input card given by the user.
*
* \internal To add a new option:
* 1. Add a member to the Config struct.
* 2. Add the option name to the "supported" Node in
* get_supported_options.
* 3. Initialise the new Config member in to_Config.
* The functions set_from_yaml (for mandatory options) and
* set_from_yaml_if_defined (non-mandatory) may be helpful.
* 4. Add a new entry (with short description) to config.yaml
*/
struct Config {
ScaleGenerator scale_gen; /**< Scale */
JetParameters resummation_jets; /**< Resummation Jets */
JetParameters fixed_order_jets; /**< Fixed-Order Jets */
double min_extparton_pt; /**< Min External Parton Pt*/
double max_ext_soft_pt_fraction; /**< Min External Parton Pt Fraction */
int trials; /**< Number of resummation points to generate per FO */
bool log_correction; /**< Log Correct On or Off */
bool unweight; /**< Unweight On or Off */
std::vector<OutputFile> output; /**< Output File Type */
optional<std::string> RanLux_init; /**< Ranlux File Choice*/
EventTreatMap treat; /**< Event Treat Map? */
YAML::Node analysis_parameters; /**< Analysis Parameters */
//! Settings for Effective Higgs-Gluon Coupling
HiggsCouplingSettings Higgs_coupling;
};
//! Load configuration from file
/**
* @param config_file Name of the YAML configuration file
* @returns The reversed HEJ configuration
*/
Config load_config(std::string const & config_file);
//! Set option using the corresponding YAML entry
/**
* @param setting Option variable to be set
* @param yaml Root of the YAML configuration
* @param names Name of the entry
*
* If the entry does not exist or has the wrong type or format
* an exception is thrown.
*
* For example
* @code
* set_from_yaml(foobar, yaml, "foo", "bar")
* @endcode
* is equivalent to
* @code
* foobar = yaml["foo"]["bar"].as<decltype(foobar)>()
* @endcode
* with improved diagnostics on errors.
*
* @see set_from_yaml_if_defined
*/
template<typename T, typename... YamlNames>
void set_from_yaml(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
);
//! Set option using the corresponding YAML entry, if present
/**
* @param setting Option variable to be set
* @param yaml Root of the YAML configuration
* @param names Name of the entry
*
* This function works similar to set_from_yaml, but does not
* throw any exception if the requested YAML entry does not exist.
*
* @see set_from_yaml
*/
template<typename T, typename... YamlNames>
void set_from_yaml_if_defined(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
);
JetParameters get_jet_parameters(
YAML::Node const & node, std::string const & entry
);
HiggsCouplingSettings get_Higgs_coupling(
YAML::Node const & node, std::string const & entry
);
ScaleConfig to_ScaleConfig(YAML::Node const & yaml);
ParticleID to_ParticleID(std::string const & particle_name);
//! Check whether all options in configuration are supported
/**
* @param conf Configuration to be checked
* @param supported Tree of supported options
*
* If conf contains an entry that does not appear in supported
* an unknown_option exception is thrown. Sub-entries of "analysis"
* (if present) are not checked.
*
* @see unknown_option
*/
void assert_all_options_known(
YAML::Node const & conf, YAML::Node const & supported
);
namespace detail{
void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml);
void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml);
void set_from_yaml(ParticleID & setting, YAML::Node const & yaml);
inline
void set_from_yaml(YAML::Node & setting, YAML::Node const & yaml){
setting = yaml;
}
template<typename Scalar>
void set_from_yaml(Scalar & setting, YAML::Node const & yaml){
assert(yaml);
if(!yaml.IsScalar()){
throw invalid_type{"value is not a scalar"};
}
try{
setting = yaml.as<Scalar>();
}
catch(...){
throw invalid_type{
"value " + yaml.as<std::string>()
+ " cannot be converted to a " + type_string(setting)
};
}
}
template<typename T>
void set_from_yaml(optional<T> & setting, YAML::Node const & yaml){
T tmp;
set_from_yaml(tmp, yaml);
setting = tmp;
}
template<typename T>
void set_from_yaml(std::vector<T> & setting, YAML::Node const & yaml){
assert(yaml);
// special case: treat a single value like a vector with one element
if(yaml.IsScalar()){
setting.resize(1);
return set_from_yaml(setting.front(), yaml);
}
if(yaml.IsSequence()){
setting.resize(yaml.size());
for(size_t i = 0; i < setting.size(); ++i){
set_from_yaml(setting[i], yaml[i]);
}
return;
}
throw invalid_type{""};
}
template<typename T, typename FirstName, typename... YamlNames>
void set_from_yaml(
T & setting,
YAML::Node const & yaml, FirstName const & name,
YamlNames && ... names
){
if(!yaml[name]) throw missing_option{""};
set_from_yaml(
setting,
yaml[name], std::forward<YamlNames>(names)...
);
}
template<typename T>
void set_from_yaml_if_defined(T & setting, YAML::Node const & yaml){
return set_from_yaml(setting, yaml);
}
template<typename T, typename FirstName, typename... YamlNames>
void set_from_yaml_if_defined(
T & setting,
YAML::Node const & yaml, FirstName const & name,
YamlNames && ... names
){
if(!yaml[name]) return;
set_from_yaml_if_defined(
setting,
yaml[name], std::forward<YamlNames>(names)...
);
}
}
template<typename T, typename... YamlNames>
void set_from_yaml(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
){
try{
detail::set_from_yaml(setting, yaml, names...);
}
catch(invalid_type const & ex){
throw invalid_type{
"In option " + join(": ", names...) + ": value " + ex.what()
+ " cannot be converted to a " + type_string(setting)
};
}
catch(missing_option const &){
throw missing_option{
"No entry for mandatory option " + join(": ", names...)
};
}
catch(std::invalid_argument const & ex){
throw missing_option{
"In option " + join(": ", names...) + ":"
" invalid value " + ex.what()
};
}
}
template<typename T, typename... YamlNames>
void set_from_yaml_if_defined(
T & setting,
YAML::Node const & yaml, YamlNames const & ... names
){
try{
detail::set_from_yaml_if_defined(setting, yaml, names...);
}
catch(invalid_type const & ex){
throw invalid_type{
"In option " + join(": ", names...) + ": " + ex.what()
};
}
catch(std::invalid_argument const & ex){
throw missing_option{
"In option " + join(": ", names...) + ":"
" invalid value " + ex.what()
};
}
}
}
namespace YAML {
template<>
struct convert<RHEJ::OutputFile> {
static Node encode(RHEJ::OutputFile const & outfile);
static bool decode(Node const & node, RHEJ::OutputFile & out);
};
}
diff --git a/include/RHEJ/event_classes.hh b/include/RHEJ/event_classes.hh
deleted file mode 100644
index 30adb8d..0000000
--- a/include/RHEJ/event_classes.hh
+++ /dev/null
@@ -1,49 +0,0 @@
-/** \file event_classes.hh
- * \brief This file details the classification of events.
- *
- * This file makes use of a macro in order to avoid repetition of EventClass names.
- * To this end, the Documentation of the EventClass 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 EventClass names
-#define RHEJ_EVENT_CLASSES(F) F(FKL), F(unordered_backward), F(unordered_forward), F(nonFKL), F(no_2_jets), F(bad_final_state)
-#define RHEJ_AS_ENUM(VAR) VAR
-#define RHEJ_AS_STRING(VAR) #VAR
-
- //! Event_Class NameSpace
- namespace event_class{
- /** \enum EventClass
- * \brief EventClass enumeration gives different possible event classes
- *
- * This enumeration is used to distinguish between different event types.
- */
- enum EventClass: size_t{
- RHEJ_EVENT_CLASSES(RHEJ_AS_ENUM), /**< Macro of Possible States */
- unob = unordered_backward, /**< Unordered Emission by backwards Jet */
- unof = unordered_forward, /**< Unordered Emission by forwards Jet */
- first_class = FKL, /**< FKL Event */
- last_class = bad_final_state /**< Bad Final State Event */
- };
-
-
- static constexpr auto names = make_array(
- RHEJ_EVENT_CLASSES(RHEJ_AS_STRING)
- );
- }
-#undef RHEJ_HISTOGRAMS
-#undef RHEJ_AS_ENUM
-#undef RHEJ_AS_STRING
-
- class Event;
-
- //! A Function used to classify a clustered Event
- event_class::EventClass classify(Event const & ev);
-
-}
diff --git a/include/RHEJ/event_types.hh b/include/RHEJ/event_types.hh
new file mode 100644
index 0000000..034dffe
--- /dev/null
+++ b/include/RHEJ/event_types.hh
@@ -0,0 +1,49 @@
+/** \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_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
+
+ class Event;
+
+ //! A Function used to typeify a clustered Event
+ event_type::EventType classify(Event const & ev);
+
+}
diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index 8bd2c2f..ff3cf3c 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,323 +1,323 @@
#include "RHEJ/EventReweighter.hh"
#include <string>
#include <unordered_map>
#include "RHEJ/PhaseSpacePoint.hh"
#include "RHEJ/PDG_codes.hh"
#include "RHEJ/config.hh" // ScaleGenerator definition
#include "RHEJ/debug.hh"
namespace RHEJ{
- using EventClass = event_class::EventClass;
+ using EventType = event_type::EventType;
namespace {
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
UnclusteredEvent result;
result.incoming = psp.incoming();
std::sort(
begin(result.incoming), end(result.incoming),
[](Sparticle o1, Sparticle o2){return o1.p.pz()<o2.p.pz();}
);
assert(result.incoming.size() == 2);
result.outgoing = psp.outgoing();
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.central.mur = NaN;
result.central.muf = NaN;
result.central.weight = psp.weight();
return result;
}
} // anonymous namespace
EventReweighter::EventReweighter(
LHEF::HEPRUP const & heprup,
fastjet::JetDefinition jet_def, double jetptmin,
EventTreatMap treat,
double extpartonptmin, double max_ext_soft_pt_fraction,
bool log_corr,
HiggsCouplingSettings Higgs_coupling
):
EventReweighter{
RHEJ::Beam{
heprup.EBMUP.first,
{{
static_cast<RHEJ::ParticleID>(heprup.IDBMUP.first),
static_cast<RHEJ::ParticleID>(heprup.IDBMUP.second)
}}
},
heprup.PDFSUP.first,
jet_def, jetptmin,
std::move(treat),
extpartonptmin, max_ext_soft_pt_fraction,
log_corr,
std::move(Higgs_coupling)
}
{
if(heprup.EBMUP.second != E_beam_){
throw std::invalid_argument(
"asymmetric beam: " + std::to_string(E_beam_)
+ " ---> <--- " + std::to_string(heprup.EBMUP.second)
);
};
if(heprup.PDFSUP.second != pdf_.id()){
throw std::invalid_argument(
"conflicting PDF ids: " + std::to_string(pdf_.id())
+ " vs. " + std::to_string(heprup.PDFSUP.second)
);
}
}
EventReweighter::EventReweighter(
Beam beam,
int pdf_id,
fastjet::JetDefinition jet_def, double jetptmin,
EventTreatMap treat,
double extpartonptmin, double max_ext_soft_pt_fraction,
bool log_corr,
HiggsCouplingSettings Higgs_coupling
):
extpartonptmin_{extpartonptmin},
max_ext_soft_pt_fraction_{max_ext_soft_pt_fraction},
jetptmin_{jetptmin},
E_beam_{beam.E},
jet_def_{jet_def},
pdf_{pdf_id, beam.type.front(), beam.type.back()},
MEt2_{
[this](double mu){ return pdf_.Halphas(mu); },
jet_def, jetptmin,
log_corr, std::move(Higgs_coupling)
},
treat_{std::move(treat)}
{}
PDF const & EventReweighter::pdf() const{
return pdf_;
}
- EventClass EventReweighter::last_event_class() const{
- return last_event_class_;
+ EventType EventReweighter::last_event_type() const{
+ return last_event_type_;
}
namespace{
static_assert(
std::numeric_limits<double>::has_infinity, "infinity not supported"
);
// there is no simple way to create a vector of move-only objects
// this is a clumsy work-around
ScaleGenerator make_identity(){
std::vector<std::unique_ptr<ScaleFunction>> scale_fun;
scale_fun.emplace_back(new InputScales);
return ScaleGenerator{
ScaleConfig{
std::move(scale_fun), {}, std::numeric_limits<double>::infinity()
}
};
}
const ScaleGenerator identity{make_identity()};
}
std::vector<Event> EventReweighter::reweight(
Event const & input_ev, int num_events
){
return reweight(input_ev, num_events, identity);
}
std::vector<Event> EventReweighter::reweight(
Event const & input_ev, int num_events,
ScaleGenerator const & gen_scales
){
auto res_events = gen_res_events(input_ev, num_events);
if(res_events.empty()) return {};
for(auto & event: res_events) event = gen_scales(event);
return rescale(input_ev, std::move(res_events));
}
/**
* \brief main generation/reweighting function: generate phase space points and divide out Born factors
*/
std::vector<Event> EventReweighter::gen_res_events(
Event const & ev,
int phase_space_points
){
assert(ev.variations().empty());
- last_event_class_ = classify(ev);
- switch(treat_.at(last_event_class_)){
+ last_event_type_ = classify(ev);
+ switch(treat_.at(last_event_type_)){
case EventTreatment::discard: return {};
case EventTreatment::keep:
if(! jets_pass_resummation_cuts(ev)) return {};
else return {ev};
default:;
}
const double Born_shat = shat(ev);
std::vector<Event> resummation_events;
for(int psp_number = 0; psp_number < phase_space_points; ++psp_number){
PhaseSpacePoint psp{
ev,
jet_def_, jetptmin_,
extpartonptmin_, max_ext_soft_pt_fraction_
};
if(psp.weight() == 0.) continue;
resummation_events.emplace_back(
to_UnclusteredEvent(std::move(psp)), jet_def_, jetptmin_
);
auto & new_event = resummation_events.back();
assert(new_event.variations().empty());
new_event.central().mur = ev.central().mur;
new_event.central().muf = ev.central().muf;
const double resum_shat = shat(new_event);
new_event.central().weight *= ev.central().weight*Born_shat*Born_shat/
(phase_space_points*resum_shat*resum_shat);
}
return resummation_events;
}
std::vector<Event> EventReweighter::rescale(
Event const & Born_ev,
std::vector<Event> events
) const{
const double Born_pdf = pdf_factors(Born_ev).central;
const double Born_ME = tree_matrix_element(Born_ev);
for(auto & cur_event: events){
const auto pdf = pdf_factors(cur_event);
assert(pdf.variations.size() == cur_event.variations().size());
const auto ME = matrix_elements(cur_event);
assert(ME.variations.size() == cur_event.variations().size());
cur_event.central().weight *= pdf.central*ME.central/(Born_pdf*Born_ME);
for(size_t i = 0; i < cur_event.variations().size(); ++i){
cur_event.variations(i).weight *=
pdf.variations[i]*ME.variations[i]/(Born_pdf*Born_ME);
}
}
return events;
};
bool EventReweighter::jets_pass_resummation_cuts(
Event const & ev
) const{
const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
fastjet::ClusterSequence cs{out_as_PseudoJet, jet_def_};
return cs.inclusive_jets(jetptmin_).size() == ev.jets().size();
}
EventReweighter::EventFactors
EventReweighter::pdf_factors(Event const & ev) const{
auto const & a = ev.incoming().front();
auto const & b = ev.incoming().back();
const double xa = a.p.e()/E_beam_;
const double xb = b.p.e()/E_beam_;
EventFactors result;
std::unordered_map<double, double> known_pdf;
result.central =
pdf_.pdfpt(0,xa,ev.central().muf,a.type)*
pdf_.pdfpt(1,xb,ev.central().muf,b.type);
known_pdf.emplace(ev.central().muf, result.central);
result.variations.reserve(ev.variations().size());
for(auto const & param: ev.variations()){
const double muf = param.muf;
auto cur_pdf = known_pdf.find(muf);
if(cur_pdf == known_pdf.end()){
cur_pdf = known_pdf.emplace(
muf,
pdf_.pdfpt(0,xa,muf,a.type)*pdf_.pdfpt(1,xb,muf,b.type)
).first;
}
result.variations.emplace_back(cur_pdf->second);
}
assert(result.variations.size() == ev.variations().size());
return result;
}
EventReweighter::EventFactors
EventReweighter::matrix_elements(Event const & ev) const{
- assert(treat_.count(last_event_class_) > 0);
- if(treat_.find(last_event_class_)->second == EventTreatment::keep){
+ assert(treat_.count(last_event_type_) > 0);
+ if(treat_.find(last_event_type_)->second == EventTreatment::keep){
return fixed_order_scale_ME(ev);
}
// precompute overall kinematic factor
const double ME_kin = MEt2_.tree_kin(ev.incoming(), ev.outgoing(), true);
EventFactors result;
std::unordered_map<double, double> known_ME;
result.central = MEt2_(
ev.central().mur,
ev.incoming(), ev.outgoing(),
true
);
known_ME.emplace(ev.central().mur, result.central);
result.variations.reserve(ev.variations().size());
for(auto const & param: ev.variations()){
const double mur = param.mur;
auto cur_ME = known_ME.find(mur);
if(cur_ME == known_ME.end()){
const double ME = MEt2_.tree_param(
mur, ev.incoming(), ev.outgoing()
)*ME_kin*MEt2_.virtual_corrections(
mur, ev.incoming(), ev.outgoing()
);
cur_ME = known_ME.emplace(mur, ME).first;
}
result.variations.emplace_back(cur_ME->second);
}
assert(result.variations.size() == ev.variations().size());
return result;
}
double EventReweighter::tree_matrix_element(Event const & ev) const{
assert(ev.variations().empty());
- assert(treat_.count(last_event_class_) > 0);
- if(treat_.find(last_event_class_)->second == EventTreatment::keep){
+ assert(treat_.count(last_event_type_) > 0);
+ if(treat_.find(last_event_type_)->second == EventTreatment::keep){
return fixed_order_scale_ME(ev).central;
}
return MEt2_.tree(
ev.central().mur,
ev.incoming(), ev.outgoing(),
false
);
}
EventReweighter::EventFactors
EventReweighter::fixed_order_scale_ME(Event const & ev) const{
const int alpha_s_power = std::count_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Sparticle const & p){ return is_parton(p); }
);
EventFactors result;
result.central = pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
for(auto const & var: ev.variations()){
result.variations.emplace_back(
pow(pdf_.Halphas(var.mur), alpha_s_power)
);
}
return result;
}
}
diff --git a/src/config.cc b/src/config.cc
index 2468110..d1be3b5 100644
--- a/src/config.cc
+++ b/src/config.cc
@@ -1,458 +1,458 @@
#include "RHEJ/config.hh"
#include "RHEJ/exceptions.hh"
#include <set>
#include <string>
#include <vector>
#include <iostream>
#include <stdexcept>
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",
"scales", "scale factors", "max scale ratio",
"log correction", "unweight", "event output", "analysis",
"RanLux init"
};
// 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] = "";
}
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;
}
}
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>());
}
}
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;
}
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 == "hepmc3") 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;
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
*/
if(name != "analysis"){
try{
assert_all_options_known(conf[name], supported[name]);
}
catch(unknown_option const & ex){
throw unknown_option{name + ": " + ex.what()};
}
}
}
}
}
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 RHEJ{
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;
}
// simple (as in non-composite) scale functions
/**
* An example for a simple scale function would be Ht,
* Ht/2 is then composite (take Ht and then divide by 2)
*/
std::unique_ptr<ScaleFunction> make_simple_ScaleFunction_ptr(
std::string const & scale_fun
){
using ret_type = std::unique_ptr<ScaleFunction>;
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
if(scale_fun == "input") return ret_type{new InputScales{}};
if(scale_fun == "Ht") return ret_type{new Ht{}};
if(scale_fun == "max jet pperp") return ret_type{new MaxJetPperp{}};
if(scale_fun == "jet invariant mass"){
return ret_type{new JetInvariantMass{}};
}
try{
const double scale = to_double(scale_fun);
return ret_type{new 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;
}
std::unique_ptr<ScaleFunction> make_ScaleFunction_ptr(
std::string const & scale_fun
){
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 make_simple_ScaleFunction_ptr(scale_fun);
}
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;
std::unique_ptr<ScaleFunction> fun;
if(scale_fun[delim] == '/'){
factor = 1/to_double(second);
fun = make_simple_ScaleFunction_ptr(first);
}
else{
assert(scale_fun[delim] == '*');
try{
factor = to_double(second);
fun = make_simple_ScaleFunction_ptr(first);
}
catch(std::invalid_argument const &){
factor = to_double(first);
fun = make_simple_ScaleFunction_ptr(second);
}
}
assert(fun != nullptr);
return std::unique_ptr<ScaleFunction>{
new Product{factor, std::move(fun)}
};
}
EventTreatMap get_event_treatment(
YAML::Node const & yaml
){
- using namespace event_class;
+ 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}
};
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"};
}
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");
set_from_yaml_if_defined(config.RanLux_init, yaml, "RanLux init");
set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
config.scale_gen = ScaleGenerator{to_ScaleConfig(yaml)};
config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
return config;
}
} // anonymous namespace
ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
ScaleConfig config;
std::vector<std::string> scales;
set_from_yaml(scales, yaml, "scales");
config.scales.reserve(scales.size());
std::transform(
begin(scales), end(scales), std::back_inserter(config.scales),
make_ScaleFunction_ptr
);
set_from_yaml_if_defined(config.scale_factors, yaml, "scale factors");
config.max_scale_ratio = std::numeric_limits<double>::infinity();
set_from_yaml_if_defined(config.max_scale_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;
}
}
Event ScaleGenerator::operator()(Event ev) const{
assert(ev.variations().empty());
assert(! settings_.scales.empty());
auto const & scales = settings_.scales;
auto const & scale_factors = settings_.scale_factors;
const double max_scale_ratio = settings_.max_scale_ratio;
// check if we are doing scale variation at all
if(scales.size() == 1 && scale_factors.empty()){
ev.central() = (*scales.front())(ev);
return ev;
}
for(auto && base_scale: scales){
const EventParameters cur_base = (*base_scale)(ev);
ev.variations().emplace_back(cur_base);
//multiplicative scale variation
for(double mur_factor: scale_factors){
const double mur = cur_base.mur*mur_factor;
for(double muf_factor: scale_factors){
if(muf_factor == 1. && mur_factor == 1.) continue;
const double muf = cur_base.muf*muf_factor;
if(mur/muf < 1/max_scale_ratio || mur/muf > max_scale_ratio) continue;
ev.variations().emplace_back(
EventParameters{mur, muf, cur_base.weight}
);
}
}
};
ev.central() = (*scales.front())(ev);
return ev;
}
}
diff --git a/src/event_classes.cc b/src/event_types.cc
similarity index 92%
rename from src/event_classes.cc
rename to src/event_types.cc
index a840fd1..1cf158f 100644
--- a/src/event_classes.cc
+++ b/src/event_types.cc
@@ -1,163 +1,163 @@
-#include "RHEJ/event_classes.hh"
+#include "RHEJ/event_types.hh"
#include <array>
#include <vector>
#include <algorithm>
#include "RHEJ/Event.hh"
namespace RHEJ{
- using EventClass = event_class::EventClass;
+ using EventType = event_type::EventType;
namespace{
// 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.
if(is_AWZH_boson(out[1])) return false;
if(!is_FKL(in, {next(begin(out)), 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 return event with all pz -> -pz
*/
Event mirrored(Event const & orig_ev){
UnclusteredEvent mirrored = orig_ev.unclustered();
for(auto & in: mirrored.incoming){
in.p.reset_momentum(in.p.px(), in.p.py(), -in.p.pz(), in.p.E());
}
std::swap(mirrored.incoming[0], mirrored.incoming[1]);
for(auto & out: mirrored.outgoing){
out.p.reset_momentum(out.p.px(), out.p.py(), -out.p.pz(), out.p.E());
}
std::reverse(begin(mirrored.outgoing), end(mirrored.outgoing));
return {mirrored, orig_ev.jet_def(), orig_ev.min_jet_pt()};
}
/**
* \brief Checks for a forward unordered gluon emission
* @param ev Event
* @returns Is the event a forward unordered emission
*
* See is_unordered_backward. Performs the same checks on the event
* after using mirrored on the event itself.
*/
bool is_unordered_forward(Event const & ev){
// check whether the mirrored event is FKL with unordered backward emission
return is_unordered_backward(mirrored(ev));
}
}
- EventClass classify(Event const & ev){
- if(! final_state_ok(ev.outgoing())) return EventClass::bad_final_state;
- if(! has_2_jets(ev)) return EventClass::no_2_jets;
- if(is_FKL(ev.incoming(), ev.outgoing())) return EventClass::FKL;
+ 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 EventClass::unordered_backward;
+ return EventType::unordered_backward;
}
if(is_unordered_forward(ev)){
- return EventClass::unordered_forward;
+ return EventType::unordered_forward;
}
- return EventClass::nonFKL;
+ return EventType::nonFKL;
}
}
diff --git a/src/main.cc b/src/main.cc
index 763cf65..78cdbdb 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -1,147 +1,147 @@
/**
* Name: main.cc
* Authors: Tuomas Hapola, Andreas Maier <andreas.maier@durham.ac.uk>
*
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include "yaml-cpp/yaml.h"
#include "LHEF/LHEF.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/get_analysis.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/EventReweighter.hh"
#include "RHEJ/config.hh"
#include "RHEJ/stream.hh"
int event_number(std::string const & record){
size_t start = record.rfind("Number of Events");
start = record.find_first_of("123456789", start);
if(start == std::string::npos) {
throw std::invalid_argument("no event number record found");
}
const size_t end = record.find_first_not_of("0123456789", start);
return std::stoi(record.substr(start, end - start));
}
RHEJ::Config load_config(char const * filename){
try{
return RHEJ::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<RHEJ::Analysis> get_analysis(
YAML::Node const & parameters
){
try{
return RHEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
int main(int argn, char** argv) {
using clock = std::chrono::system_clock;
if (argn < 3) {
std::cerr << "\n# Usage:\n./HEJ config_file input_file\n\n";
return EXIT_FAILURE;
}
const auto start_time = clock::now();
// read configuration
const RHEJ::Config config = load_config(argv[1]);
RHEJ::istream in{argv[2]};
std::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
LHEF::Reader reader{in};
RHEJ::CombinedEventWriter writer{config.output, reader.heprup};
double global_reweight = 1.;
int max_events = std::numeric_limits<int>::max();
if(argn > 3){
max_events = std::stoi(argv[3]);
const int input_events = event_number(reader.headerBlock);
global_reweight = input_events/static_cast<double>(max_events);
std::cout << "processing " << max_events
<< " out of " << input_events << " events\n";
}
RHEJ::EventReweighter rhej{
reader.heprup,
config.resummation_jets.def, config.resummation_jets.min_pt,
config.treat,
config.min_extparton_pt, config.max_ext_soft_pt_fraction,
config.log_correction,
config.Higgs_coupling
};
if(config.RanLux_init){
RHEJ::PhaseSpacePoint::reset_ranlux(*config.RanLux_init);
}
int nevent = 0;
- std::array<int, RHEJ::event_class::last_class + 1>
- nevent_class{0}, nfailed_class{0};
+ std::array<int, RHEJ::event_type::last_type + 1>
+ nevent_type{0}, nfailed_type{0};
// Loop over the events in the inputfile
while(reader.readEvent()){
// reweight events so that the total cross section is conserved
reader.hepeup.setWeight(0, global_reweight * reader.hepeup.weight());
if(nevent == max_events) break;
nevent++;
if (nevent % 10000 == 0)
std::cout << "event number " << nevent << std::endl;
// calculate rHEJ weight
RHEJ::Event FO_event{
RHEJ::UnclusteredEvent{reader.hepeup},
config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
};
auto resummed_events = rhej.reweight(
FO_event,
config.trials, config.scale_gen
);
- const auto ev_class = rhej.last_event_class();
- ++nevent_class[ev_class];
+ const auto ev_type = rhej.last_event_type();
+ ++nevent_type[ev_type];
- if(resummed_events.empty()) ++nfailed_class[ev_class];
+ if(resummed_events.empty()) ++nfailed_type[ev_type];
for(auto const & ev: resummed_events){
//TODO: move pass_cuts to after phase space point generation
if(analysis->pass_cuts(ev)){
analysis->fill(ev);
writer.write(ev);
}
}
} // main event loop
- using namespace RHEJ::event_class;
+ using namespace RHEJ::event_type;
std::cout<< "Events processed: " << nevent << '\n';
- for(size_t ev_class = first_class; ev_class <= last_class; ++ev_class){
- std::cout << '\t' << names[ev_class] << ": " << nevent_class[ev_class]
- << ", failed to reconstruct " << nfailed_class[ev_class]
+ for(size_t ev_type = first_type; ev_type <= last_type; ++ev_type){
+ std::cout << '\t' << names[ev_type] << ": " << nevent_type[ev_type]
+ << ", failed to reconstruct " << nfailed_type[ev_type]
<< '\n';
}
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n";
}
diff --git a/t/check_res.cc b/t/check_res.cc
index 87dba75..623636a 100644
--- a/t/check_res.cc
+++ b/t/check_res.cc
@@ -1,71 +1,71 @@
#include <iostream>
#include "LHEF/LHEF.h"
#include "RHEJ/stream.hh"
#include "RHEJ/EventReweighter.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_class;
+ using namespace RHEJ::event_type;
RHEJ::EventTreatMap treat{
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{nonFKL, 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::EventReweighter rhej{
reader.heprup,
jet_def, jetptmin,
treat,
extpartonptmin, max_ext_soft_pt_fraction,
log_corr,
RHEJ::HiggsCouplingSettings{}
};
double xsec = 0.;
while(reader.readEvent()){
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;
}
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/test_classify.cc b/t/test_classify.cc
index c888d61..dbea328 100644
--- a/t/test_classify.cc
+++ b/t/test_classify.cc
@@ -1,51 +1,51 @@
#include "LHEF/LHEF.h"
#include "RHEJ/stream.hh"
-#include "RHEJ/event_classes.hh"
+#include "RHEJ/event_types.hh"
#include "RHEJ/Event.hh"
namespace{
constexpr double min_jet_pt = 30.;
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
- using namespace RHEJ::event_class;
- static const std::vector<EventClass> results{
+ using namespace RHEJ::event_type;
+ static const std::vector<EventType> results{
unob,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,unob,FKL,FKL,FKL,unof,FKL,unob,FKL,
FKL,unob,unob,FKL,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unof,
FKL,FKL,unof,FKL,FKL,FKL,FKL,FKL,unof,FKL,FKL,FKL,unof,FKL,FKL,unob,unof,
FKL,unof,FKL,unob,FKL,FKL,unob,FKL,unob,unof,unob,unof,FKL,FKL,FKL,FKL,FKL,
FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,
FKL,FKL,FKL,unof,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,
FKL,FKL,FKL,FKL,unob,FKL,unob,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unof,unob,FKL
};
}
int main(int argn, char** argv) {
if(argn != 2){
std::cerr << "Usage: test_classify eventfile";
return EXIT_FAILURE;
}
RHEJ::istream in{argv[1]};
LHEF::Reader reader{in};
LHEF::Writer writer{std::cerr};
writer.heprup = reader.heprup;
for(auto const & expected: results){
reader.readEvent();
const RHEJ::Event ev{
RHEJ::UnclusteredEvent{reader.hepeup},
jet_def, min_jet_pt
};
- const auto ev_class = classify(ev);
- if(ev_class != expected){
- using RHEJ::event_class::names;
+ const auto ev_type = classify(ev);
+ if(ev_type != expected){
+ using RHEJ::event_type::names;
writer.hepeup = reader.hepeup;
std::cerr << "wrong classification of event:\n";
writer.writeEvent();
- std::cerr << "classified as " << names[ev_class]
+ std::cerr << "classified as " << names[ev_type]
<< ", is " << names[expected] << '\n';
return EXIT_FAILURE;
}
}
}
diff --git a/t/test_psp.cc b/t/test_psp.cc
index 4fa423c..28fe4dc 100644
--- a/t/test_psp.cc
+++ b/t/test_psp.cc
@@ -1,64 +1,64 @@
#include "LHEF/LHEF.h"
#include "RHEJ/stream.hh"
-#include "RHEJ/event_classes.hh"
+#include "RHEJ/event_types.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/PhaseSpacePoint.hh"
namespace{
constexpr int max_trials = 100;
constexpr double extpartonptmin = 45.;
constexpr double max_ext_soft_pt_fraction =
std::numeric_limits<double>::infinity();
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
constexpr double min_jet_pt = 50;
};
int main(int argn, char** argv) {
if(argn != 2){
std::cerr << "Usage: test_psp eventfile";
return EXIT_FAILURE;
}
RHEJ::istream in{argv[1]};
LHEF::Reader reader{in};
LHEF::Writer writer{std::cerr};
writer.heprup = reader.heprup;
while(reader.readEvent()){
const RHEJ::Event ev{
RHEJ::UnclusteredEvent{reader.hepeup},
jet_def, min_jet_pt
};
- const auto in_class = classify(ev);
+ const auto in_type = classify(ev);
for(int trial = 0; trial < max_trials; ++trial){
RHEJ::PhaseSpacePoint psp{
ev, jet_def, min_jet_pt,
extpartonptmin, max_ext_soft_pt_fraction
};
if(psp.weight() != 0){
RHEJ::UnclusteredEvent tmp_ev;
tmp_ev.incoming = psp.incoming();
tmp_ev.outgoing = psp.outgoing();
tmp_ev.central = {0,0,0};
RHEJ::Event out_ev{
std::move(tmp_ev),
jet_def, min_jet_pt
};
- const auto out_class = classify(out_ev);
- if(out_class != in_class){
- using RHEJ::event_class::names;
+ const auto out_type = classify(out_ev);
+ if(out_type != in_type){
+ using RHEJ::event_type::names;
std::cerr << "Wrong class of phase space point:\n"
- "original event of class " << names[in_class] << ":\n";
+ "original event of class " << names[in_type] << ":\n";
writer.hepeup = reader.hepeup;
writer.writeEvent();
- std::cerr << "Phase space point of class " << names[out_class] << ":\n";
+ std::cerr << "Phase space point of class " << names[out_type] << ":\n";
writer.hepeup = to_HEPEUP(out_ev, &writer.heprup);
writer.writeEvent();
return EXIT_FAILURE;
}
}
}
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 2:19 AM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4238779
Default Alt Text
(67 KB)

Event Timeline