Page MenuHomeHEPForge
Files F8725643

No OneTemporary
Actions

View Options

diff --git a/include/RHEJ/EventReweighter.hh b/include/RHEJ/EventReweighter.hh
index 4de8cf2..6985f6c 100644
--- a/include/RHEJ/EventReweighter.hh
+++ b/include/RHEJ/EventReweighter.hh
@@ -1,192 +1,192 @@
/** \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/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;
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 */
);
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 */
);
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 */
);
- //! Constructor for last_event_class
+ //! The EventClass for the Last Event
EventClass last_event_class() 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 */
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 bf6a3ba..dd63eab 100644
--- a/include/RHEJ/config.hh
+++ b/include/RHEJ/config.hh
@@ -1,106 +1,112 @@
/** \file config.hh
* \brief The file which handles the configuration file parameters
*
- * Contains the JetParameters Struct and ScaleConfig Struct. Also
+ * 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 "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 */
+ 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.
+ * 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
+ *
+ * If called upon with EventTreatMap.at(EventClass) it will return the treatment which has
+ * been specified for that EventClass.
+ */
using EventTreatMap = std::map<event_class::EventClass, 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.
*/
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 */
};
Config load_config(std::string const & config_file);
}
diff --git a/src/Event.cc b/src/Event.cc
index 40f9a13..35a046d 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,137 +1,145 @@
#include "RHEJ/Event.hh"
#include "RHEJ/debug.hh"
namespace RHEJ{
namespace{
constexpr int status_in = -1;
constexpr int status_out = 1;
}
UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup):
central(EventParameters{
hepeup.scales.mur, hepeup.scales.muf, hepeup.weight()
})
{
size_t in_idx = 0;
for (int i = 0; i < hepeup.NUP; ++i) {
// Only consider final state particles
if (abs(hepeup.ISTUP[i]) != status_out) continue;
// Save particle information
Sparticle particle{
static_cast<ParticleID>(hepeup.IDUP[i]),
fastjet::PseudoJet{
hepeup.PUP[i][0], hepeup.PUP[i][1],
hepeup.PUP[i][2], hepeup.PUP[i][3]
}
};
if (hepeup.ISTUP[i]>0) outgoing.emplace_back(std::move(particle));
else{
if(in_idx > incoming.size()) {
throw std::invalid_argument{
"Event has too many incoming particles"
};
}
incoming[in_idx++] = std::move(particle);
}
}
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{});
}
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}
{}
std::vector<fastjet::PseudoJet> Event::jets() const{
return cs_.inclusive_jets(min_jet_pt_);
}
+ /**
+ * \brief Returns the invarient mass of the event
+ * @param ev Event
+ * @returns s hat
+ *
+ * Makes use of the FastJet PseudoJet function m2().
+ * Applies this function to the sum of the incoming partons.
+ */
double shat(Event const & ev){
return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
}
namespace{
// colour flow according to Les Houches standard
// TODO: stub
std::vector<std::pair<int, int>> colour_flow(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> const & outgoing
){
std::vector<std::pair<int, int>> result(
incoming.size() + outgoing.size()
);
for(auto & col: result){
col = std::make_pair(-1, -1);
}
return result;
}
}
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){
LHEF::HEPEUP result;
result.heprup = heprup;
result.weights = {{event.central().weight, nullptr}};
for(auto const & var: event.variations()){
result.weights.emplace_back(var.weight, nullptr);
}
const size_t num_particles =
event.incoming().size() + event.outgoing().size();
result.NUP = num_particles;
// the following entries are pretty much meaningless
result.IDPRUP = 1;
result.AQEDUP = 1./128.;
// end meaningless part
result.XWGTUP = event.central().weight;
result.SCALUP = event.central().muf;
result.scales.muf = event.central().muf;
result.scales.mur = event.central().mur;
result.scales.SCALUP = event.central().muf;
result.pdfinfo.p1 = event.incoming().front().type;
result.pdfinfo.p2 = event.incoming().back().type;
result.pdfinfo.scale = event.central().muf;
for(Sparticle const & in: event.incoming()){
result.IDUP.emplace_back(in.type);
result.ISTUP.emplace_back(status_in);
result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()});
result.MOTHUP.emplace_back(0, 0);
}
for(Sparticle const & out: event.outgoing()){
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()});
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)
);
}
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/EventReweighter.cc b/src/EventReweighter.cc
index a02567d..2e645ed 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,316 +1,317 @@
#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;
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
):
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
}
{
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
):
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_{jet_def, jetptmin, log_corr},
treat_{std::move(treat)}
{}
PDF const & EventReweighter::pdf() const{
return pdf_;
}
EventClass EventReweighter::last_event_class() const{
return last_event_class_;
}
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));
}
- // main generation/reweighting function
- // generate phase space points and divide out Born factors
+ /**
+ * \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_)){
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){
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_(
pdf_.Halphas(ev.central().mur), 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 alpha_s = pdf_.Halphas(mur);
const double ME = MEt2_.tree_param(
alpha_s, mur, ev.incoming(), ev.outgoing()
)*ME_kin*MEt2_.virtual_corrections(
alpha_s, 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){
return fixed_order_scale_ME(ev).central;
}
return MEt2_.tree(
pdf_.Halphas(ev.central().mur), 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/event_classes.cc b/src/event_classes.cc
index fd5920f..a840fd1 100644
--- a/src/event_classes.cc
+++ b/src/event_classes.cc
@@ -1,140 +1,163 @@
#include "RHEJ/event_classes.hh"
#include <array>
#include <vector>
#include <algorithm>
#include "RHEJ/Event.hh"
namespace RHEJ{
using EventClass = event_class::EventClass;
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;
}
- // return event with all pz -> -pz
+ /**
+ * \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;
if(is_unordered_backward(ev)){
return EventClass::unordered_backward;
}
if(is_unordered_forward(ev)){
return EventClass::unordered_forward;
}
return EventClass::nonFKL;
}
}

File Metadata

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

Event Timeline