Page MenuHomeHEPForge
Files F7879917

No OneTemporary
Actions

View Options

diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index a995969..5c7c8cf 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,268 +1,289 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/EventReweighter.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <functional>
#include <string>
#include <unordered_map>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/PseudoJet.hh"
#include "LHEF/LHEF.h"
#include "HEJ/Event.hh"
#include "HEJ/Fraction.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PhaseSpacePoint.hh"
#include "HEJ/exceptions.hh"
namespace HEJ {
EventReweighter::EventReweighter(
LHEF::HEPRUP const & heprup,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
std::shared_ptr<RNG> ran
):
EventReweighter{
Beam{
heprup.EBMUP.first,
{{
static_cast<ParticleID>(heprup.IDBMUP.first),
static_cast<ParticleID>(heprup.IDBMUP.second)
}}
},
heprup.PDFSUP.first,
std::move(scale_gen),
std::move(conf),
std::move(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 const & beam,
int pdf_id,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
std::shared_ptr<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_{std::move(ran)}
{
// legacy code: override new variable with old
if(param_.psp_config.max_ext_soft_pt_fraction){
param_.psp_config.soft_pt_regulator = *param_.psp_config.max_ext_soft_pt_fraction;
param_.psp_config.max_ext_soft_pt_fraction = {};
}
assert(ran_);
}
PDF const & EventReweighter::pdf() const{
return pdf_;
}
std::vector<Event> EventReweighter::reweight(
Event const & input_ev, std::size_t 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_(std::move(event));
return rescale(input_ev, std::move(res_events));
}
EventTreatment EventReweighter::treatment(EventType type) const {
return param_.treat.at(type);
}
std::vector<Event> EventReweighter::gen_res_events(
Event const & ev,
std::size_t phase_space_points
){
assert(ev.variations().empty());
status_.clear();
switch(treatment(ev.type())){
case EventTreatment::discard: {
status_.emplace_back(StatusCode::discard);
return {};
}
case EventTreatment::keep:
if(! jets_pass_resummation_cuts(ev)) {
status_.emplace_back(StatusCode::failed_resummation_cuts);
return {};
}
else {
status_.emplace_back(StatusCode::good);
return {ev};
}
default:;
}
const double Born_shat = shat(ev);
std::vector<Event> resummation_events;
status_.reserve(phase_space_points);
for(std::size_t psp_number = 0; psp_number < phase_space_points; ++psp_number){
PhaseSpacePoint psp{ev, param_.psp_config, *ran_};
status_.emplace_back(psp.status());
assert(psp.status() != StatusCode::unspecified);
if(psp.status() != StatusCode::good) continue;
assert(psp.weight() != 0.);
if(psp.incoming()[0].E() > E_beam_ || psp.incoming()[1].E() > E_beam_) {
status_.back() = StatusCode::too_much_energy;
continue;
}
resummation_events.emplace_back(
to_EventData( std::move(psp) ).cluster(
param_.jet_param().def, param_.jet_param().min_pt
)
);
auto & new_event = resummation_events.back();
assert( new_event.valid_hej_state(
param_.psp_config.soft_pt_regulator,
param_.psp_config.min_extparton_pt ) );
if( new_event.type() != ev.type() ) {
throw std::logic_error{
"Resummation Event does not match Born event: "
+ name(new_event.type())
+ " != "
+ name(ev.type())
};
}
new_event.generate_colours(*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);
+
+ // Option to only keep lowpt events
+ // Keep only events where there is a fixed order event with at least 1
+ // jet below the resummation jet pt but all resummation jets are above
+ // the resummation jet pt
+ bool passes_pt_cut = true;
+ bool jet_below_cut = false;
+
+ for (size_t jet_it = 0; jet_it < ev.jets().size(); ++jet_it){
+ jet_below_cut = jet_below_cut || ev.jets()[jet_it].pt()
+ < param_.jet_param().min_pt;
+ passes_pt_cut = passes_pt_cut && new_event.jets()[jet_it].pt()
+ > param_.jet_param().min_pt;
+ }
+
+ if (!passes_pt_cut || !jet_below_cut) {
+ new_event.central().weight *= 0;
+ status_.emplace_back(StatusCode::too_much_energy);
+ return {};
+ }
+
}
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.parameters() *= pdf*ME/(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 if(part.type == pid::Higgs) {
alpha_s_power += 2;
}
// nothing to do for other uncoloured particles
}
Weights result;
result.central = std::pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
for(auto const & var: ev.variations()){
result.variations.emplace_back(
std::pow(pdf_.Halphas(var.mur), alpha_s_power)
);
}
return result;
}
} // namespace HEJ

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 9:22 PM (23 h, 1 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806264
Default Alt Text
(9 KB)

Event Timeline