ScaleFunction.cc
View Options

	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/ScaleFunction.hh"

	#include <cassert>

	#include "HEJ/Event.hh"
	#include "HEJ/exceptions.hh"

	namespace HEJ{

	double H_T(Event const & ev){
	double result = 0.;
	for(size_t i = 0; i < ev.outgoing().size(); ++i){
	auto const decay_products = ev.decays().find(i);
	if(decay_products == end(ev.decays())){
	result += ev.outgoing()[i].perp();
	}
	else{
	for(auto const & particle: decay_products->second){
	result += particle.perp();
	}
	}
	}
	return result;
	}

	double max_jet_pt(Event const & ev) {
	return sorted_by_pt(ev.jets()).front().pt();
	}

	double jet_invariant_mass(Event const & ev) {
	fastjet::PseudoJet sum;
	for(const auto & jet: ev.jets()) sum+=jet;
	return sum.m();
	}

	double m_j1j2(Event const & ev) {
	const auto jets = sorted_by_pt(ev.jets());
	assert(jets.size() >= 2);
	return (jets[0] + jets[1]).m();
	}

	ScaleFunction operator*(double factor, ScaleFunction base_scale) {
	base_scale.name_.insert(0, std::to_string(factor) + '*');
	auto new_fun =
	[factor,fun{std::move(base_scale.fun_)}](HEJ::Event const & ev) {
	return factor*fun(ev);
	};
	base_scale.fun_ = std::move(new_fun);
	return base_scale;
	}

	ScaleFunction operator*(ScaleFunction factor, ScaleFunction base_scale) {
	base_scale.name_.insert(0, factor.name_ + '*');
	auto new_fun =
	[fun1{std::move(factor.fun_)}, fun2{std::move(base_scale.fun_)}, name{base_scale.name_}]
	(HEJ::Event const & ev) {
	return fun1(ev)*fun2(ev);
	};
	base_scale.fun_ = std::move(new_fun);
	return base_scale;
	}

	ScaleFunction operator/(ScaleFunction base_scale, double denom) {
	base_scale.name_.append('/' + std::to_string(denom));
	auto new_fun =
	[denom,fun{std::move(base_scale.fun_)}](HEJ::Event const & ev) {
	return fun(ev)/denom;
	};
	base_scale.fun_ = std::move(new_fun);
	return base_scale;
	}

	ScaleFunction operator/(ScaleFunction base_scale, ScaleFunction denom) {
	base_scale.name_.append('/' + denom.name_);
	auto new_fun =
	[fun2{std::move(denom.fun_)}, fun1{std::move(base_scale.fun_)}]
	(HEJ::Event const & ev) {
	return fun1(ev)/fun2(ev);
	};
	base_scale.fun_ = std::move(new_fun);
	return base_scale;
	}

	// TODO: significant logic duplication with operator()
	void ScaleGenerator::gen_descriptions() {
	if(scales_.empty()) {
	throw std::logic_error{"Need at least one scale"};
	}
	descriptions_.emplace_back(
	std::make_shared<ParameterDescription>(scales_.front().name(), 1., 1.)
	);

	for(auto & base_scale: scales_){
	const auto base_name = base_scale.name();
	descriptions_.emplace_back(
	std::make_shared<ParameterDescription>(base_name, 1., 1.)
	);
	//multiplicative scale variation
	for(double mur_factor: scale_factors_){
	for(double muf_factor: scale_factors_){
	if(muf_factor == 1. && mur_factor == 1.) continue;
	if(
	mur_factor/muf_factor < 1/max_scale_ratio_
	\|\| mur_factor/muf_factor > max_scale_ratio_
	) continue;
	descriptions_.emplace_back(
	std::make_shared<ParameterDescription>(
	base_name, mur_factor, muf_factor
	)
	);
	}
	}
	}
	}

	Event ScaleGenerator::operator()(Event ev) const {
	if(! ev.variations().empty()) {
	throw std::invalid_argument{"Input event already has scale variation"};
	}
	assert(!scales_.empty());
	assert(!descriptions_.empty());

	size_t descr_idx = 0;
	const double mu_central = (scales_.front())(ev);
	ev.central().mur = mu_central;
	ev.central().muf = mu_central;
	assert(descr_idx < descriptions_.size());
	assert(descriptions_[descr_idx]);
	ev.central().description = descriptions_[descr_idx++];

	// check if we are doing scale variation at all
	if(scales_.size() == 1 && scale_factors_.empty()) return ev;

	for(auto & base_scale: scales_){
	const double mu_base = base_scale(ev);
	assert(descr_idx < descriptions_.size());
	assert(descriptions_[descr_idx]);
	ev.variations().emplace_back(
	EventParameters{
	mu_base, mu_base, ev.central().weight, descriptions_[descr_idx++]
	}
	);
	//multiplicative scale variation
	for(double mur_factor: scale_factors_){
	const double mur = mu_base*mur_factor;
	for(double muf_factor: scale_factors_){
	if(muf_factor == 1. && mur_factor == 1.) continue;
	const double muf = mu_base*muf_factor;
	if(
	mur/muf < 1/max_scale_ratio_
	\|\| mur/muf > max_scale_ratio_
	) continue;
	assert(descr_idx < descriptions_.size());
	assert(descriptions_[descr_idx]);
	ev.variations().emplace_back(
	EventParameters{
	mur, muf, ev.central().weight, descriptions_[descr_idx++]
	}
	);
	}
	}
	};
	return ev;
	}

	}

File Metadata

Mime Type: text/x-c
Expires: Tue, Sep 30, 6:10 AM (19 h, 18 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6459787
Default Alt Text: ScaleFunction.cc (5 KB)

ScaleFunction.cc
No OneTemporary
Actions

ScaleFunction.cc
View Options

File Metadata

Event Timeline

ScaleFunction.ccNo OneTemporaryActions

ScaleFunction.ccView Options

File Metadata

Event Timeline

ScaleFunction.cc
No OneTemporary
Actions

ScaleFunction.cc
View Options