Page MenuHomeHEPForge
Files F8724529

No OneTemporary
Actions

View Options

diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc
index eb99958..cc35e45 100644
--- a/FixedOrderGen/src/EventGenerator.cc
+++ b/FixedOrderGen/src/EventGenerator.cc
@@ -1,80 +1,82 @@
#include "EventGenerator.hh"
#include "Process.hh"
#include "Beam.hh"
#include "JetParameters.hh"
#include "PhaseSpacePoint.hh"
#include "HEJ/Event.hh"
#include "HEJ/config.hh"
namespace HEJFOG{
EventGenerator::EventGenerator(
Process process,
Beam beam,
HEJ::ScaleGenerator scale_gen,
JetParameters jets,
int pdf_id,
double subl_change,
unsigned int subl_channels,
ParticlesPropMap particles_properties,
HEJ::HiggsCouplingSettings Higgs_coupling,
HEJ::RNG & ran
):
pdf_{pdf_id, beam.particles[0], beam.particles[1]},
ME_{
[this](double mu){ return pdf_.Halphas(mu); },
HEJ::MatrixElementConfig{
false,
std::move(Higgs_coupling)
}
},
scale_gen_{std::move(scale_gen)},
process_{std::move(process)},
jets_{std::move(jets)},
beam_{std::move(beam)},
subl_change_{subl_change},
subl_channels_{subl_channels},
particles_properties_{std::move(particles_properties)},
ran_{ran}
{
}
HEJ::Event EventGenerator::gen_event(){
HEJFOG::PhaseSpacePoint psp{
process_,
jets_,
pdf_, beam_.energy,
subl_change_, subl_channels_,
particles_properties_,
ran_
};
status_ = psp.status();
if(status_ != good) return {};
HEJ::Event ev = scale_gen_(
HEJ::Event{
to_UnclusteredEvent(std::move(psp)),
jets_.def, jets_.min_pt
}
);
+ ev.generate_colour_flow(ran_);
+
const double shat = HEJ::shat(ev);
const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy);
const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy);
// evaluate matrix element on this point
const auto ME_weights = ME_.tree(ev);
ev.central().weight *= ME_weights.central/(shat*shat);
ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type);
ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type);
for(size_t i = 0; i < ev.variations().size(); ++i){
auto & var = ev.variations(i);
var.weight *= ME_weights.variations[i]/(shat*shat);
var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type);
var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type);
}
return ev;
}
}
diff --git a/include/HEJ/Particle.hh b/include/HEJ/Particle.hh
index 5b287a1..67af9db 100644
--- a/include/HEJ/Particle.hh
+++ b/include/HEJ/Particle.hh
@@ -1,144 +1,143 @@
/**
* \file Particle.hh
* \brief Contains the particle struct
*
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <utility>
#include "fastjet/PseudoJet.hh"
#include "HEJ/optional.hh"
#include "HEJ/PDG_codes.hh"
namespace HEJ {
using Colour = std::pair<int,int>;
//! Class representing a particle
struct Particle {
//! particle type
ParticleID type;
//! particle momentum
fastjet::PseudoJet p;
//! (optional) colour & anti-colour
- //! @TODO resolve compiler warnings of "missing initializer for member"
optional<Colour> colour;
//! get rapidity
double rapidity() const{
return p.rapidity();
}
//! get transverse momentum
double perp() const{
return p.perp();
}
//! get momentum in x direction
double px() const{
return p.px();
}
//! get momentum in y direction
double py() const{
return p.py();
}
//! get momentum in z direction
double pz() const{
return p.pz();
}
//! get energy
double E() const{
return p.E();
}
//! get mass
double m() const{
return p.m();
}
};
//! Functor to compare rapidities
/**
* This can be used whenever a rapidity comparison function is needed,
* for example in many standard library functions.
*
* @see pz_less
*/
struct rapidity_less{
template<class FourVector>
bool operator()(FourVector const & p1, FourVector const & p2){
return p1.rapidity() < p2.rapidity();
}
};
//! Functor to compare momenta in z direction
/**
* This can be used whenever a pz comparison function is needed,
* for example in many standard library functions.
*
* @see rapidity_less
*/
struct pz_less{
template<class FourVector>
bool operator()(FourVector const & p1, FourVector const & p2){
return p1.pz() < p2.pz();
}
};
//! Convert a vector of Particles to a vector of particle momenta
inline
std::vector<fastjet::PseudoJet> to_PseudoJet(
std::vector<Particle> const & v
){
std::vector<fastjet::PseudoJet> result;
for(auto && sp: v) result.emplace_back(sp.p);
return result;
}
//! Check if a particle is a parton, i.e. quark, antiquark, or gluon
inline
bool is_parton(Particle const & p){
return is_parton(p.type);
}
//! Check if a particle is a quark
inline
bool is_quark(Particle const & p){
return is_quark(p.type);
}
//! Check if a particle is an anti-quark
inline
bool is_antiquark(Particle const & p){
return is_antiquark(p.type);
}
//! Check if a particle is a quark or anit-quark
inline
bool is_anyquark(Particle const & p){
return is_anyquark(p.type);
}
//! Check if a particle is a photon, W, Z, or Higgs boson
inline bool is_AWZH_boson(Particle const & particle){
return is_AWZH_boson(particle.type);
}
//! Extract all partons from a vector of particles
inline
std::vector<Particle> filter_partons(
std::vector<Particle> const & v
){
std::vector<Particle> result;
result.reserve(v.size());
std::copy_if(
begin(v), end(v), std::back_inserter(result),
[](Particle const & p){ return is_parton(p); }
);
return result;
}
}
diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index 88ac45d..a8e6b00 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,275 +1,276 @@
/**
* \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();
UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
UnclusteredEvent 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_UnclusteredEvent(std::move(psp)),
param_.jet_param.def, param_.jet_param.min_pt
);
auto & new_event = resummation_events.back();
+ new_event.generate_colour_flow(ran_);
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

File Metadata

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

Event Timeline