Page MenuHomeHEPForge
Files F8724499

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/EventReweighter.hh b/include/HEJ/EventReweighter.hh
index b6b1c2d..b183a91 100644
--- a/include/HEJ/EventReweighter.hh
+++ b/include/HEJ/EventReweighter.hh
@@ -1,190 +1,190 @@
/** \file
* \brief Declares the EventReweighter class
*
* EventReweighter is the main class used within HEJ 2. It reweights the
* resummation events.
*
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <functional>
#include <utility>
#include <vector>
#include "HEJ/config.hh"
#include "HEJ/event_types.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/PDF.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/RNG.hh"
#include "HEJ/ScaleFunction.hh"
+#include "HEJ/Weights.hh"
namespace LHEF {
class HEPRUP;
}
namespace HEJ{
class Event;
- class Weights;
//! Beam parameters
/**
* Currently, only symmetric beams are supported,
* so there is a single beam energy.
*/
struct Beam{
double E; /**< Beam energy */
std::array<ParticleID, 2> type; /**< Beam particles */
};
//! Main class for reweighting events in HEJ.
class EventReweighter{
using EventType = event_type::EventType;
public:
EventReweighter(
Beam beam, /**< Beam Energy */
int pdf_id, /**< PDF ID */
ScaleGenerator scale_gen, /**< Scale settings */
EventReweighterConfig conf, /**< Configuration parameters */
HEJ::RNG & ran /**< Random number generator */
);
EventReweighter(
LHEF::HEPRUP const & heprup, /**< LHEF event header */
ScaleGenerator scale_gen, /**< Scale settings */
EventReweighterConfig conf, /**< Configuration parameters */
HEJ::RNG & ran /**< Random number generator */
);
//! Get the used pdf
PDF const & pdf() const;
//! Generate resummation events for a given fixed-order event
/**
* @param ev Fixed-order event corresponding
* to the resummation events
* @param num_events Number of trial resummation configurations.
* @returns A vector of resummation events.
*
* The result vector depends on the type of the input event and the
* treatment of different types as specified in the constructor:
*
* \ref reweight The result vector contains between
* 0 and num_events resummation events.
*
* \ref keep If the input event passes the resummation jet cuts
* the result vector contains one event. Otherwise it is empty.
*
* \ref discard The result vector is empty
*/
std::vector<Event> reweight(
Event const & ev,
int num_events
);
private:
template<typename... T>
PDF const & pdf(T&& ...);
/** \internal
* \brief main generation/reweighting function:
* generate phase space points and divide out Born factors
*/
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;
/** \internal
* \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;
/** \internal
* \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.
*/
Weights pdf_factors(Event const & ev) const;
/** \internal
* \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.
*/
Weights matrix_elements(Event const & ev) const;
/** \internal
* \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.
*/
Weights fixed_order_scale_ME(Event const & ev) const;
/** \internal
* \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;
//! \internal General parameters
EventReweighterConfig param_;
//! \internal Beam energy
double E_beam_;
//! \internal PDF
PDF pdf_;
//! \internal Object to calculate the square of the matrix element
MatrixElement MEt2_;
//! \internal Object to calculate event renormalisation and factorisation scales
ScaleGenerator scale_gen_;
/** \internal random number generator
*
* \note We use a reference_wrapper so that EventReweighter objects can
* still be copied (which would be impossible with a reference).
*/
std::reference_wrapper<HEJ::RNG> ran_;
};
template<typename... T>
PDF const & EventReweighter::pdf(T&&... t){
return pdf_ = PDF{std::forward<T>(t)...};
}
}
diff --git a/include/HEJ/Weights.hh b/include/HEJ/Weights.hh
index fd34e92..8de138e 100644
--- a/include/HEJ/Weights.hh
+++ b/include/HEJ/Weights.hh
@@ -1,49 +1,97 @@
/** \file
- * \brief Container for event weights
+ * \brief Container for event Weights
*
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
+#include "HEJ/exceptions.hh"
+
#include <vector>
namespace HEJ{
- //! Collection of weights assigned to a single event
+ //! Collection of parameters, e.g. Weights, assigned to a single event
/**
- * A number of member functions of the MatrixElement class return Weights
+ * A number of member functions of the MatrixElement class return Parameters
* objects containing the squares of the matrix elements for the various
* scale choices.
*/
- struct Weights {
- double central;
- std::vector<double> variations;
-
- Weights& operator*=(Weights const & other);
- Weights& operator*=(double factor);
- Weights& operator/=(Weights const & other);
- Weights& operator/=(double factor);
+ template<class T>
+ struct Parameters {
+ T central;
+ std::vector<T> variations;
+
+ template<class T_ext>
+ Parameters<T>& operator*=(Parameters<T_ext> const & other);
+ Parameters<T>& operator*=(double factor);
+ template<class T_ext>
+ Parameters<T>& operator/=(Parameters<T_ext> const & other);
+ Parameters<T>& operator/=(double factor);
};
- inline
- Weights operator*(Weights a, Weights const & b) {
+ template<class T1, class T2> inline
+ Parameters<T1> operator*(Parameters<T1> a, Parameters<T2> const & b) {
return a*=b;
}
- inline
- Weights operator*(Weights a, double b) {
+ template<class T> inline
+ Parameters<T> operator*(Parameters<T> a, double b) {
return a*=b;
}
- inline
- Weights operator*(double b, Weights a) {
+ template<class T> inline
+ Parameters<T> operator*(double b, Parameters<T> a) {
return a*=b;
}
- inline
- Weights operator/(Weights a, Weights const & b) {
+ template<class T1, class T2> inline
+ Parameters<T1> operator/(Parameters<T1> a, Parameters<T2> const & b) {
return a/=b;
}
- inline
- Weights operator/(Weights a, double b) {
+ template<class T> inline
+ Parameters<T> operator/(Parameters<T> a, double b) {
return a/=b;
}
+
+ //! Shortcut for Weights, e.g. used by the MatrixElement
+ using Weights = Parameters<double>;
+
+ template<class T>
+ template<class T_ext>
+ Parameters<T>& Parameters<T>::operator*=(Parameters<T_ext> const & other) {
+ if(other.variations.size() != variations.size()) {
+ throw std::invalid_argument{"Wrong number of Parameters"};
+ }
+ central *= other.central;
+ for(std::size_t i = 0; i < variations.size(); ++i) {
+ variations[i] *= other.variations[i];
+ }
+ return *this;
+ };
+
+ template<class T>
+ Parameters<T>& Parameters<T>::operator*=(double factor) {
+ central *= factor;
+ for(auto & wt: variations) wt *= factor;
+ return *this;
+ };
+
+ template<class T>
+ template<class T_ext>
+ Parameters<T>& Parameters<T>::operator/=(Parameters<T_ext> const & other) {
+ if(other.variations.size() != variations.size()) {
+ throw std::invalid_argument{"Wrong number of Parameters"};
+ }
+ central /= other.central;
+ for(std::size_t i = 0; i < variations.size(); ++i) {
+ variations[i] /= other.variations[i];
+ }
+ return *this;
+ };
+
+ template<class T>
+ Parameters<T>& Parameters<T>::operator/=(double factor) {
+ central /= factor;
+ for(auto & wt: variations) wt /= factor;
+ return *this;
+ };
}
diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index c888370..a23923e 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,275 +1,274 @@
/**
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/EventReweighter.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <string>
#include <unordered_map>
#include "fastjet/ClusterSequence.hh"
#include "LHEF/LHEF.h"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/PhaseSpacePoint.hh"
-#include "HEJ/Weights.hh"
namespace HEJ{
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();
EventData to_EventData(PhaseSpacePoint const & psp){
EventData result;
result.incoming() = psp.incoming();
std::sort(
begin(result.incoming()), end(result.incoming()),
[](Particle o1, Particle 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.decays() = psp.decays();
result.central().mur = NaN;
result.central().muf = NaN;
result.central().weight = psp.weight();
return result;
}
} // namespace anonymous
EventReweighter::EventReweighter(
LHEF::HEPRUP const & heprup,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
HEJ::RNG & ran
):
EventReweighter{
HEJ::Beam{
heprup.EBMUP.first,
{{
static_cast<HEJ::ParticleID>(heprup.IDBMUP.first),
static_cast<HEJ::ParticleID>(heprup.IDBMUP.second)
}}
},
heprup.PDFSUP.first,
std::move(scale_gen),
std::move(conf),
ran
}
{
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,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
HEJ::RNG & ran
):
param_{std::move(conf)},
E_beam_{beam.E},
pdf_{pdf_id, beam.type.front(), beam.type.back()},
MEt2_{
[this](double mu){ return pdf_.Halphas(mu); },
param_.ME_config
},
scale_gen_(std::move(scale_gen)),
ran_{ran}
{}
PDF const & EventReweighter::pdf() const{
return pdf_;
}
std::vector<Event> EventReweighter::reweight(
Event const & input_ev, int num_events
){
auto res_events = gen_res_events(input_ev, num_events);
if(res_events.empty()) return {};
for(auto & event: res_events) event = scale_gen_(event);
return rescale(input_ev, std::move(res_events));
}
std::vector<Event> EventReweighter::gen_res_events(
Event const & ev,
int phase_space_points
){
assert(ev.variations().empty());
switch(param_.treat.at(ev.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, param_.psp_config, ran_};
if(psp.weight() == 0.) continue;
if(psp.incoming()[0].E() > E_beam_ || psp.incoming()[1].E() > E_beam_) continue;
resummation_events.emplace_back(
to_EventData(std::move(psp)),
param_.jet_param.def, param_.jet_param.min_pt
);
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, param_.jet_param.def};
return cs.inclusive_jets(param_.jet_param.min_pt).size() == ev.jets().size();
}
Weights
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_;
Weights 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 & ev_param: ev.variations()){
const double muf = ev_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;
}
Weights
EventReweighter::matrix_elements(Event const & ev) const{
assert(param_.treat.count(ev.type()) > 0);
if(param_.treat.find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev);
}
return MEt2_(ev);
}
double EventReweighter::tree_matrix_element(Event const & ev) const{
assert(ev.variations().empty());
assert(param_.treat.count(ev.type()) > 0);
if(param_.treat.find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev).central;
}
return MEt2_.tree(ev).central;
}
Weights
EventReweighter::fixed_order_scale_ME(Event const & ev) const{
int alpha_s_power = 0;
for(auto const & part: ev.outgoing()){
if(is_parton(part))
++alpha_s_power;
else {
switch(part.type){
case pid::Higgs: {
alpha_s_power += 2;
break;
}
// TODO
case pid::Wp:
case pid::Wm:
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
}
Weights 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;
}
} // namespace HEJ
diff --git a/src/Weights.cc b/src/Weights.cc
deleted file mode 100644
index de2d561..0000000
--- a/src/Weights.cc
+++ /dev/null
@@ -1,46 +0,0 @@
-/**
- * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
- * \date 2019
- * \copyright GPLv2 or later
- */
-#include "HEJ/Weights.hh"
-
-#include <stdexcept>
-
-namespace HEJ {
-
- Weights& Weights::operator*=(Weights const & other) {
- if(other.variations.size() != variations.size()) {
- throw std::invalid_argument{"Wrong number of weights"};
- }
- central *= other.central;
- for(std::size_t i = 0; i < variations.size(); ++i) {
- variations[i] *= other.variations[i];
- }
- return *this;
- }
-
- Weights& Weights::operator*=(double factor) {
- central *= factor;
- for(auto & wt: variations) wt *= factor;
- return *this;
- }
-
- Weights& Weights::operator/=(Weights const & other) {
- if(other.variations.size() != variations.size()) {
- throw std::invalid_argument{"Wrong number of weights"};
- }
- central /= other.central;
- for(std::size_t i = 0; i < variations.size(); ++i) {
- variations[i] /= other.variations[i];
- }
- return *this;
- }
-
- Weights& Weights::operator/=(double factor) {
- central /= factor;
- for(auto & wt: variations) wt /= factor;
- return *this;
- }
-
-}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 11:34 PM (1 d, 9 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242905
Default Alt Text
(19 KB)

Event Timeline