EventReweighter.cc
No OneTemporary
Actions

Size

7 KB

Referenced Files

None

Subscribers

None

EventReweighter.cc
View Options

	/**
	* \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"};
	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().weightBorn_shat*Born_shat/
	(phase_space_pointsresum_shatresum_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.parameters() = pdfME/(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-c
Expires: Tue, Sep 30, 5:50 AM (1 d, 18 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6537521
Default Alt Text: EventReweighter.cc (7 KB)

EventReweighter.ccNo OneTemporaryActions

EventReweighter.ccView Options

File Metadata

Event Timeline

EventReweighter.cc
No OneTemporary
Actions

EventReweighter.cc
View Options