No OneTemporary
Actions

Size

42 KB

Subscribers

None

View Options

	diff --git a/FixedOrderGen/configFO.yml b/FixedOrderGen/configFO.yml
	index 26ccc5c..0633304 100644
	--- a/FixedOrderGen/configFO.yml
	+++ b/FixedOrderGen/configFO.yml
	@@ -1,69 +1,68 @@
	# number of generated events
	events: 200

	jets:
	min pt: 20
	peak pt: 30
	algorithm: antikt
	R: 0.4
	max rapidity: 5

	beam:
	energy: 6500
	particles: [p, p]

	pdf: 230000

	process: p p => h 4j

	# fraction of events with two extremal emissions in one direction
	# that contain an subleading emission e.g. unordered emission
	subleading fraction: 0.05
	# Allow different subleading configurations
	# By default all implemented processes are allowed
	#
	-# subleading channels: unordered
	+# subleading channels: [unordered, qqx]

	scales: max jet pperp

	event output:
	- HEJFO.lhe
	-# - HEJ.lhe
	-# - HEJ_events.hepmc
	+# - HEJFO_events.hepmc

	particle properties:
	Higgs:
	mass: 125
	width: 0.004165
	decays: {into: [photon, photon], branching ratio: 0.0023568762400521404}

	random generator:
	name: mixmax
	# seed: 1

	# unweighting parameters
	# remove to obtain weighted events
	unweight:
	sample size: 200 # should be similar to "events:", but not more than ~10000
	max deviation: 0

	# to use a rivet analysis
	#
	# analysis:
	# rivet: MC_XS # rivet analysis name
	# output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added
	#
	# to use a custom analysis
	#
	# analysis:
	# plugin: /path/to/libmyanalysis.so
	# my analysis parameter: some value

	# parameters for Higgs-gluon couplings
	# this requires compilation with qcdloop
	#
	# Higgs coupling:
	# use impact factors: false
	# mt: 174
	# include bottom: true
	# mb: 4.7
	diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
	index b471662..93ac98c 100644
	--- a/FixedOrderGen/src/PhaseSpacePoint.cc
	+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
	@@ -1,665 +1,663 @@
	#include "PhaseSpacePoint.hh"

	#include <random>
	#include <algorithm>

	#include "HEJ/Constants.hh"
	#include "HEJ/kinematics.hh"
	#include "HEJ/utility.hh"
	#include "HEJ/exceptions.hh"

	#include "Process.hh"
	#include "Subleading.hh"

	#include <CLHEP/Random/Randomize.h>
	#include <CLHEP/Random/RanluxEngine.h>

	using namespace HEJ;

	namespace HEJFOG{

	static_assert(
	std::numeric_limits<double>::has_quiet_NaN,
	"no quiet NaN for double"
	);
	constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
	HEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
	HEJ::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),
	HEJ::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{
	bool can_swap_to_uno(
	HEJ::Particle const & p1, HEJ::Particle const & p2
	){
	return is_parton(p1)
	&& p1.type != pid::gluon
	&& p2.type == pid::gluon;
	}

	size_t count_gluons(std::vector<Particle>::const_iterator first,
	std::vector<Particle>::const_iterator last){
	return std::count_if(first, last, [](Particle const & p)
	{return p.type == pid::gluon;});
	}

	/** assumes FKL configurations between first and last,
	* else there can be a quark in a non-extreme position
	* e.g. uno configuration gqg would pass
	*/
	bool can_change_to_qqx(
	std::vector<Particle>::const_iterator first,
	std::vector<Particle>::const_iterator last){
	return 1 < count_gluons(first,last);
	}
	bool is_AWZ_proccess(Process const & proc){
	return proc.boson && is_AWZ_boson(*proc.boson);
	}

	bool is_up_type(Particle const & part){
	return HEJ::is_anyquark(part) && !(abs(part.type)%2);
	}
	bool is_down_type(Particle const & part){
	return HEJ::is_anyquark(part) && abs(part.type)%2;
	}
	/// true iff parton can couple to a W
	bool can_couple_to_W(Particle const & part, int const sign_W){
	return abs(part.type)<5
	&& ( (sign_W*part.type > 0 && is_up_type(part))
	\|\| (sign_W*part.type < 0 && is_down_type(part)) );
	}
	}

	void PhaseSpacePoint::maybe_turn_to_subl(
	double chance,
	unsigned int const channels,
	Process const & proc,
	HEJ::RNG & ran
	){
	if(proc.njets <= 2) return;
	assert(outgoing_.size() >= 2);

	// decide what kind of subleading process is allowed
	bool allow_uno = false;
	bool allow_strange = true;
	const size_t nout = outgoing_.size();
	- const bool can_be_uno_backward = channels&Subleading::uno?
	- can_swap_to_uno(outgoing_[0], outgoing_[1]):false;
	- const bool can_be_uno_forward = channels&Subleading::uno?
	- can_swap_to_uno(outgoing_[nout-1], outgoing_[nout-2]):false;
	+ const bool can_be_uno_backward = (channels&Subleading::uno)
	+ && can_swap_to_uno(outgoing_[0], outgoing_[1]);
	+ const bool can_be_uno_forward = (channels&Subleading::uno)
	+ && can_swap_to_uno(outgoing_[nout-1], outgoing_[nout-2]);
	allow_uno = can_be_uno_backward \|\| can_be_uno_forward;

	bool allow_qqx = false;
	if(is_AWZ_proccess(proc)) {
	- allow_qqx = channels&Subleading::qqx?
	- can_change_to_qqx(outgoing_.cbegin(), outgoing_.cend()):false;
	+ allow_qqx = (channels&Subleading::qqx)
	+ && can_change_to_qqx(outgoing_.cbegin(), outgoing_.cend());
	const int sign_W = *proc.boson>0?1:-1;
	if(std::none_of(outgoing_.cbegin(), outgoing_.cend(),
	[sign_W](Particle const & p){ return can_couple_to_W(p, sign_W);})) {
	- if(!(channels&Subleading::qqx))
	- throw std::logic_error(
	- "Invalid configuration, require qqx pair but channel forbidden");
	// enforce qqx if A/W/Z can't couple somewhere else
	+ assert(allow_qqx);
	allow_uno = false;
	chance = 1.;
	// strange not allowed for W
	if(abs(*proc.boson)== pid::Wp) allow_strange = false;
	}
	}

	if(!allow_uno && !allow_qqx) return;
	if(ran.flat() < chance){
	weight_ /= chance;
	if(allow_uno && !allow_qqx){
	turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
	} else if (!allow_uno && allow_qqx) {
	turn_to_qqx(allow_strange, ran);
	} else {
	assert( allow_uno && allow_qqx);
	if(ran.flat() < 0.5) turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
	else turn_to_qqx(allow_strange, ran);
	weight_ *= 2.;
	}
	} else weight_ /= 1 - chance;

	}

	void PhaseSpacePoint::turn_to_uno(
	const bool can_be_uno_backward, const bool can_be_uno_forward,
	HEJ::RNG & ran
	){
	if(!can_be_uno_backward && !can_be_uno_forward) return;
	const size_t nout = outgoing_.size();
	if(can_be_uno_backward && can_be_uno_forward){
	if(ran.flat() < 0.5){
	std::swap(outgoing_[0].type, outgoing_[1].type);
	} else {
	std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
	}
	weight_ *= 2.;
	} else if(can_be_uno_backward){
	std::swap(outgoing_[0].type, outgoing_[1].type);
	} else {
	assert(can_be_uno_forward);
	std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
	}
	}

	void PhaseSpacePoint::turn_to_qqx(const bool allow_stange, HEJ::RNG & ran){
	/// find first and last gluon in FKL chain
	auto first = std::find_if(outgoing_.begin(), outgoing_.end(),
	[](Particle const & p){return p.type == pid::gluon;});
	std::vector<Particle*> FKL_gluons;
	for(auto p = first; p<outgoing_.end(); ++p){
	if((p).type == pid::gluon) FKL_gluons.push_back(&p);
	else if(is_quark(p) \|\| is_antiquark(p)) break;
	}
	const size_t ng = FKL_gluons.size();
	if(ng < 2)
	throw std::logic_error("not enough gluons to create qqx");
	// select flavour of quark
	const double r1 = 2.*ran.flat()-1.;
	const double max_flavour = allow_stange?n_f:n_f-1;
	weight_ = max_flavour2;
	int flavour = pid::down;
	for (double sum = 1./max_flavour; sum < std::abs(r1); sum += 1./max_flavour)
	++flavour;
	flavour*=r1<0.?-1:1;
	// select gluon for switch
	const size_t idx = floor((ng-1) * ran.flat());
	weight_ *= (ng-1);
	FKL_gluons[idx]->type = ParticleID(flavour);
	FKL_gluons[idx+1]->type = ParticleID(-flavour);
	}

	template<class ParticleMomenta>
	fastjet::PseudoJet PhaseSpacePoint::gen_last_momentum(
	ParticleMomenta const & other_momenta,
	const double mass_square, const double y
	) const {
	std::array<double,2> pt{0.,0.};
	for (auto const & p: other_momenta) {
	pt[0]-= p.px();
	pt[1]-= p.py();
	}

	const double mperp = sqrt(pt[0]pt[0]+pt[1]pt[1]+mass_square);
	const double pz=mperp*sinh(y);
	const double E=mperp*cosh(y);

	return {pt[0], pt[1], pz, E};
	}

	PhaseSpacePoint::PhaseSpacePoint(
	Process const & proc,
	JetParameters const & jet_param,
	HEJ::PDF & pdf, double E_beam,
	double const subl_chance,
	unsigned int const subl_channels,
	ParticlesPropMap const & particles_properties,
	HEJ::RNG & ran
	)
	{
	assert(proc.njets >= 2);
	if(proc.boson
	&& particles_properties.find(*(proc.boson))
	== particles_properties.end())
	throw HEJ::missing_option("Boson "
	+std::to_string(*(proc.boson))+" can't be generated: missing properties");
	status_ = good;
	weight_ = 1;

	const int nout = proc.njets + (proc.boson?1:0);
	outgoing_.reserve(nout);

	// generate parton momenta
	const bool is_pure_jets = !proc.boson;
	auto partons = gen_LO_partons(
	proc.njets, is_pure_jets, jet_param, E_beam, ran
	);
	// pre fill flavour with gluons
	for(auto&& p_out: partons) {
	outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out)});
	}
	if(status_ != good) return;

	// create boson
	if(proc.boson){
	const auto & boson_prop = particles_properties.at(*proc.boson);
	auto boson(gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran));
	const auto pos = std::upper_bound(
	begin(outgoing_),end(outgoing_),boson,rapidity_less{}
	);
	outgoing_.insert(pos, std::move(boson));
	if(! boson_prop.decays.empty()){
	const size_t boson_idx = std::distance(begin(outgoing_), pos);
	decays_.emplace(
	boson_idx,
	decay_boson(outgoing_[boson_idx], boson_prop.decays, ran)
	);
	}
	}
	// normalisation of momentum-conserving delta function
	weight_ = pow(2M_PI, 4);

	reconstruct_incoming(proc, subl_channels, pdf,E_beam,jet_param.min_pt, ran);
	if(status_ != good) return;
	// set outgoing states
	most_backward_FKL(outgoing_).type = incoming_[0].type;
	most_forward_FKL(outgoing_).type = incoming_[1].type;

	maybe_turn_to_subl(subl_chance, subl_channels, proc, ran);

	if(proc.boson) couple_boson(*proc.boson, ran);
	}

	double PhaseSpacePoint::gen_hard_pt(
	int np , double ptmin, double ptmax, double y,
	HEJ::RNG & ran
	) {
	// heuristic parameters for pt sampling
	const double ptpar = ptmin + np/5.;
	const double arg_small_y = atan((ptmax - ptmin)/ptpar);
	const double y_cut = 3.;

	const double r1 = ran.flat();
	if(y < y_cut){
	const double pt = ptmin + ptpartan(r1arg_small_y);
	const double temp = cos(r1*arg_small_y);
	weight_ = ptptpararg_small_y/(temptemp);
	return pt;
	}

	const double ptpar2 = ptpar/(1 + 5*(y-y_cut));
	const double temp = 1. - std::exp((ptmin-ptmax)/ptpar2);
	const double pt = ptmin - ptpar2std::log(1-r1temp);
	weight_ = ptptpar2temp/(1-r1temp);
	return pt;
	}

	double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, HEJ::RNG & ran) {
	constexpr double ptpar = 4.;

	const double r = ran.flat();
	const double pt = max_pt + ptpar/np*std::log(r);
	weight_ = ptptpar/(np*r);
	return pt;
	}

	double PhaseSpacePoint::gen_parton_pt(
	int count, JetParameters const & jet_param, double max_pt, double y,
	HEJ::RNG & ran
	) {
	constexpr double p_small_pt = 0.02;

	if(! jet_param.peak_pt) {
	return gen_hard_pt(count, jet_param.min_pt, max_pt, y, ran);
	}
	const double r = ran.flat();
	if(r > p_small_pt) {
	weight_ /= 1. - p_small_pt;
	return gen_hard_pt(count, *jet_param.peak_pt, max_pt, y, ran);
	}
	weight_ /= p_small_pt;
	const double pt = gen_soft_pt(count, *jet_param.peak_pt, ran);
	if(pt < jet_param.min_pt) {
	weight_=0.0;
	status_ = not_enough_jets;
	return jet_param.min_pt;
	}
	return pt;
	}

	std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_LO_partons(
	int np, bool is_pure_jets,
	JetParameters const & jet_param,
	double max_pt,
	HEJ::RNG & ran
	){
	if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"};

	weight_ /= pow(16.*pow(M_PI,3),np);
	weight_ /= std::tgamma(np+1); //remove rapidity ordering
	std::vector<fastjet::PseudoJet> partons;
	partons.reserve(np);
	for(int i = 0; i < np; ++i){
	const double y = -jet_param.max_y + 2jet_param.max_yran.flat();
	weight_ = 2jet_param.max_y;

	const bool is_last_parton = i+1 == np;
	if(is_pure_jets && is_last_parton) {
	constexpr double parton_mass_sq = 0.;
	partons.emplace_back(gen_last_momentum(partons, parton_mass_sq, y));
	break;
	}

	const double phi = 2M_PIran.flat();
	weight_ = 2.0M_PI;

	const double pt = gen_parton_pt(np, jet_param, max_pt, y, ran);
	if(weight_ == 0.0) return {};

	partons.emplace_back(fastjet::PtYPhiM(pt, y, phi));
	assert(jet_param.min_pt <= partons[i].pt());
	assert(partons[i].pt() <= max_pt+1e-5);
	}
	// Need to check that at LO, the number of jets = number of partons;
	fastjet::ClusterSequence cs(partons, jet_param.def);
	auto cluster_jets=cs.inclusive_jets(jet_param.min_pt);
	if (cluster_jets.size()!=unsigned(np)){
	weight_=0.0;
	status_ = not_enough_jets;
	return {};
	}

	std::sort(begin(partons), end(partons), rapidity_less{});

	return partons;
	}

	Particle PhaseSpacePoint::gen_boson(
	HEJ::ParticleID bosonid, double mass, double width,
	HEJ::RNG & ran
	){
	// Usual phase space measure
	weight_ /= 16.*pow(M_PI, 3);

	// Generate a y Gaussian distributed around 0
	/// @TODO: magic number only for Higgs
	/// @TODO better sampling for W
	const double y = random_normal(1.6, ran);
	const double r1 = ran.flat();
	const double sH = mass*(
	mass + widthtan(M_PI/2.r1 + (r1-1.)*atan(mass/width))
	);

	auto p = gen_last_momentum(outgoing_, sH, y);

	return Particle{bosonid, std::move(p)};
	}

	Particle const & PhaseSpacePoint::most_backward_FKL(
	std::vector<Particle> const & partons
	) const{
	if(!HEJ::is_parton(partons[0])) return partons[1];
	return partons[0];
	}

	Particle const & PhaseSpacePoint::most_forward_FKL(
	std::vector<Particle> const & partons
	) const{
	const size_t last_idx = partons.size() - 1;
	if(!HEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
	return partons[last_idx];
	}

	Particle & PhaseSpacePoint::most_backward_FKL(
	std::vector<Particle> & partons
	) const{
	if(!HEJ::is_parton(partons[0])) return partons[1];
	return partons[0];
	}

	Particle & PhaseSpacePoint::most_forward_FKL(
	std::vector<Particle> & partons
	) const{
	const size_t last_idx = partons.size() - 1;
	if(!HEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
	return partons[last_idx];
	}

	namespace {
	/// partons are ordered: even = anti, 0 = gluon
	ParticleID index_to_pid(size_t i){
	if(!i) return pid::gluon;
	return static_cast<ParticleID>(i%2?(i+1)/2:-i/2);
	}

	/// partons are ordered: even = anti, 0 = gluon
	size_t pid_to_index(ParticleID id){
	if(id==pid::gluon) return 0;
	return id>0?id2-1:abs(id)2;
	}

	std::bitset<11> init_allowed(ParticleID const id){
	if(abs(id) == pid::proton)
	return ~0;
	std::bitset<11> out = 0;
	if(is_parton(id))
	out[pid_to_index(id)] = 1;
	return out;
	}

	/// decides which "index" (see index_to_pid) are allowed for process
	std::bitset<11> allowed_quarks(ParticleID const boson){
	std::bitset<11> allowed = ~0;
	if(abs(boson) == pid::Wp){
	// special case W:
	// Wp: anti-down or up-type quark, no b/t -> 0001100110(1) = 205
	// Wm: down or anti-up-type quark, no b/t -> 0010011001(1) = 307
	allowed = boson>0?205:307;
	}
	return allowed;
	}
	}

	/**
	* checks which partons are allowed as initial state:
	* 1. only allow what is given in the Runcard (p -> all)
	* 2. A/W/Z require something to couple to
	* a) no qqx => no incoming gluon
	* b) 2j => no incoming gluon
	* c) 3j => can couple OR is gluon => 2 gluons become qqx later
	*/
	std::array<std::bitset<11>,2> PhaseSpacePoint::filter_partons(
	Process const & proc, unsigned int const subl_channels, HEJ::RNG & ran
	){
	std::array<std::bitset<11>,2> allowed_partons{
	init_allowed(proc.incoming[0]),
	init_allowed(proc.incoming[1])
	};
	bool const allow_qqx = subl_channels&Subleading::qqx;
	// special case A/W/Z
	if(is_AWZ_proccess(proc) && ((proc.njets < 4) \|\| !allow_qqx)){
	// all possible incoming states
	auto allowed(allowed_quarks(*proc.boson));
	if(proc.njets == 2 \|\| !allow_qqx) allowed[0]=0;

	// possible states per leg
	std::array<std::bitset<11>,2> const maybe_partons{
	allowed_partons[0]&allowed, allowed_partons[1]&allowed};

	if(maybe_partons[0].any() && maybe_partons[1].any()){
	// two options to get allowed initial state => choose one at random
	const size_t idx = ran.flat() < 0.5;
	allowed_partons[idx] = maybe_partons[idx];
	+ // else choose the possible
	} else if(maybe_partons[0].any()) {
	- // only first possible => choose
	allowed_partons[0] = maybe_partons[0];
	} else if(maybe_partons[1].any()) {
	- // only second possible => choose
	allowed_partons[1] = maybe_partons[1];
	} else{
	throw std::invalid_argument{"Incoming state not allowed."};
	}
	}
	return allowed_partons;
	}

	void PhaseSpacePoint::reconstruct_incoming(
	Process const & proc, unsigned int const subl_channels,
	HEJ::PDF & pdf, double E_beam,
	double uf,
	HEJ::RNG & ran
	){
	std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
	// calculate xa, xb
	const double sqrts=2*E_beam;
	const double xa=(incoming_[0].p.e()-incoming_[0].p.pz())/sqrts;
	const double xb=(incoming_[1].p.e()+incoming_[1].p.pz())/sqrts;
	// abort if phase space point is outside of collider energy reach
	if (xa>1. \|\| xb>1.){
	weight_=0;
	status_ = too_much_energy;
	return;
	}
	// pick pdfs
	/** @TODO
	* ufa, ufb don't correspond to our final scale choice.
	* The HEJ scale generators currently expect a full event as input,
	* so fixing this is not completely trivial
	*/
	auto const & ids = proc.incoming;
	std::array<std::bitset<11>,2> allowed_partons(
	filter_partons(proc, subl_channels, ran));
	for(size_t i = 0; i < 2; ++i){
	if(ids[i] == pid::proton \|\| ids[i] == pid::p_bar){
	incoming_[i].type =
	generate_incoming_id(i, i?xb:xa, uf, pdf, allowed_partons[i], ran);
	} else {
	assert(allowed_partons[i][pid_to_index(ids[i])]);
	incoming_[i].type = ids[i];
	}
	}
	assert(momentum_conserved(1e-7));
	}

	HEJ::ParticleID PhaseSpacePoint::generate_incoming_id(
	size_t const beam_idx, double const x, double const uf,
	HEJ::PDF & pdf, std::bitset<11> allowed_partons, HEJ::RNG & ran
	){
	std::array<double,11> pdf_wt;
	pdf_wt[0] = allowed_partons[0]?fabs(pdf.pdfpt(beam_idx,x,uf,pid::gluon)):0.;
	double pdftot = pdf_wt[0];
	for(size_t i = 1; i < pdf_wt.size(); ++i){
	pdf_wt[i] = allowed_partons[i]?4./9.*fabs(pdf.pdfpt(beam_idx,x,uf,index_to_pid(i))):0;
	pdftot += pdf_wt[i];
	}
	const double r1 = pdftot * ran.flat();
	double sum = 0;
	for(size_t i=0; i < pdf_wt.size(); ++i){
	if (r1 < (sum+=pdf_wt[i])){
	weight_*= pdftot/pdf_wt[i];
	return index_to_pid(i);
	}
	}
	std::cerr << "Error in choosing incoming parton: "<<x<<" "<<uf<<" "
	<<sum<<" "<<pdftot<<" "<<r1<<std::endl;
	throw std::logic_error{"Failed to choose parton flavour"};
	}

	void PhaseSpacePoint::couple_boson(
	HEJ::ParticleID const boson, HEJ::RNG & ran
	){
	if(abs(boson) != pid::Wp) return; // only matters for W
	+ /// @TODO this could be use to sanity check gamma and Z

	// find all possible quarks
	const int sign_W = boson>0?1:-1;
	std::vector<Particle*> allowed_parts;
	for(auto & part: outgoing_){
	// Wp -> up OR anti-down, Wm -> anti-up OR down, no bottom
	if ( can_couple_to_W(part, sign_W) )
	allowed_parts.push_back(&part);
	}
	if(allowed_parts.size() == 0){
	throw std::logic_error{"Found no parton for coupling with boson"};
	}

	// select one and flip it
	size_t idx = 0;
	if(allowed_parts.size() > 1){
	/// @TODO more efficient sampling
	/// old code: probability[i] = exp(parton[i].y - W.y)
	idx = floor(ran.flat()*allowed_parts.size());
	weight_ *= allowed_parts.size();
	}
	allowed_parts[idx]->type =
	static_cast<ParticleID>( allowed_parts[idx]->type - sign_W );
	}

	double PhaseSpacePoint::random_normal(
	double stddev,
	HEJ::RNG & ran
	){
	const double r1 = ran.flat();
	const double r2 = ran.flat();
	const double lninvr1 = -log(r1);
	const double result = stddevsqrt(2.lninvr1)cos(2.M_PI*r2);
	weight_ = exp(resultresult/(2stddevstddev))sqrt(2.M_PI)*stddev;
	return result;
	}

	bool PhaseSpacePoint::momentum_conserved(double ep) const{
	fastjet::PseudoJet diff;
	for(auto const & in: incoming()) diff += in.p;
	for(auto const & out: outgoing()) diff -= out.p;
	return nearby_ep(diff, fastjet::PseudoJet{}, ep);
	}

	Decay PhaseSpacePoint::select_decay_channel(
	std::vector<Decay> const & decays,
	HEJ::RNG & ran
	){
	double br_total = 0.;
	for(auto const & decay: decays) br_total += decay.branching_ratio;
	// adjust weight
	// this is given by (channel branching ratio)/(chance to pick channel)
	// where (chance to pick channel) =
	// (channel branching ratio)/(total branching ratio)
	weight_ *= br_total;
	const double r1 = br_total*ran.flat();
	double br_sum = 0.;
	for(auto const & decay: decays){
	br_sum += decay.branching_ratio;
	if(r1 < br_sum) return decay;
	}
	throw std::logic_error{"unreachable"};
	}

	std::vector<Particle> PhaseSpacePoint::decay_boson(
	HEJ::Particle const & parent,
	std::vector<Decay> const & decays,
	HEJ::RNG & ran
	){
	const auto channel = select_decay_channel(decays, ran);
	if(channel.products.size() != 2){
	throw HEJ::not_implemented{
	"only decays into two particles are implemented"
	};
	}
	std::vector<Particle> decay_products(channel.products.size());
	for(size_t i = 0; i < channel.products.size(); ++i){
	decay_products[i].type = channel.products[i];
	}
	// choose polar and azimuth angle in parent rest frame
	const double E = parent.m()/2;
	const double theta = 2.M_PIran.flat();
	const double cos_phi = 2.*ran.flat()-1.;
	const double sin_phi = sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
	const double px = Ecos(theta)sin_phi;
	const double py = Esin(theta)sin_phi;
	const double pz = E*cos_phi;
	decay_products[0].p.reset(px, py, pz, E);
	decay_products[1].p.reset(-px, -py, -pz, E);
	for(auto & particle: decay_products) particle.p.boost(parent.p);
	return decay_products;
	}
	}
	diff --git a/FixedOrderGen/t/W_nj_classify.cc b/FixedOrderGen/t/W_nj_classify.cc
	index b3cb317..06ff18f 100644
	--- a/FixedOrderGen/t/W_nj_classify.cc
	+++ b/FixedOrderGen/t/W_nj_classify.cc
	@@ -1,177 +1,177 @@
	#ifdef NDEBUG
	#undef NDEBUG
	#endif

	#include <algorithm>

	#include "JetParameters.hh"
	#include "ParticleProperties.hh"
	#include "PhaseSpacePoint.hh"
	#include "Process.hh"
	#include "Subleading.hh"

	#include "HEJ/Event.hh"
	#include "HEJ/Mixmax.hh"
	#include "HEJ/PDF.hh"
	#include "HEJ/utility.hh"

	using namespace HEJFOG;
	using namespace HEJ;

	namespace {
	void print_psp(PhaseSpacePoint const & psp){
	std::cerr << "Process:\n"
	<< psp.incoming()[0].type << " + "<< psp.incoming()[1].type << " -> ";
	for(auto const & out: psp.outgoing()){
	std::cerr << out.type << " ";
	}
	std::cerr << "\n";
	}
	void bail_out(PhaseSpacePoint const & psp, std::string msg){
	print_psp(psp);
	throw std::logic_error{msg};
	}
	}

	int main(){
	- constexpr size_t n_psp_base = 19375;
	+ constexpr size_t n_psp_base = 10375;
	const JetParameters jet_para{
	fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30, 5, 30};
	PDF pdf(11000, pid::proton, pid::proton);
	constexpr double E_cms = 13000.;
	constexpr double subl_change = 0.8;
	const ParticlesPropMap boson_prop{
	{pid::Wp, {91.1876, 2.085, {Decay{ {pid::e_bar, pid::nu_e}, 1.}} }},
	{pid::Wm, {91.1876, 2.085, {Decay{ {pid::e, pid::nu_e_bar}, 1.}} }}
	};
	HEJ::Mixmax ran{};

	auto subl_channels = Subleading::all;
	std::vector<event_type::EventType> allowed_types{event_type::FKL,
	event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf};

	std::cout << "Wp3j" << std::endl;
	// Wp3j
	Process proc {{pid::proton,pid::proton}, 3, pid::Wp};
	size_t n_psp = n_psp_base;
	std::unordered_map<event_type::EventType, size_t> type_counter;
	for( size_t i = 0; i<n_psp; ++i){
	const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
	boson_prop, ran};
	if(psp.status()==good){
	const Event ev(to_UnclusteredEvent(psp), jet_para.def, jet_para.min_pt);
	++type_counter[ev.type()];
	if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
	== allowed_types.cend()) {
	bail_out(psp, "Found not allowed event of type "
	+std::string(event_type::names[ev.type()]));
	}
	} else { // bad process -> try again
	++n_psp;
	}
	}
	std::cout << "Wp+3j: Took " << n_psp << " to generate "
	<< n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
	std::cout << "States by classification:\n";
	for(auto const & entry: type_counter){
	const double fraction = static_cast<double>(entry.second)/n_psp_base;
	const int percent = std::round(100*fraction);
	std::cout << std::left << std::setw(25)
	<< (event_type::names[entry.first] + std::string(":"))
	<< entry.second << " (" << percent << "%)\n";

	}
	for(auto const & t: allowed_types){
	if(type_counter[t] < 0.05 * n_psp_base){
	std::cerr << "Less than 5% of the events are of type " << event_type::names[t] << std::endl;
	return EXIT_FAILURE;
	}
	}

	// Wm3j - only uno
	proc = Process{{pid::proton,pid::proton}, 3, pid::Wm};
	n_psp = n_psp_base;

	subl_channels = Subleading::uno;
	allowed_types = {event_type::FKL, event_type::unob, event_type::unof};
	type_counter.clear();

	for( size_t i = 0; i<n_psp; ++i){
	const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
	boson_prop, ran};
	if(psp.status()==good){
	const Event ev(to_UnclusteredEvent(psp), jet_para.def, jet_para.min_pt);
	++type_counter[ev.type()];
	if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
	== allowed_types.cend()) {
	bail_out(psp, "Found not allowed event of type "
	+std::string(event_type::names[ev.type()]));
	}
	} else { // bad process -> try again
	++n_psp;
	}
	}
	std::cout << "Wm+3j (only uno): Took " << n_psp << " to generate "
	<< n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
	std::cout << "States by classification:\n";
	for(auto const & entry: type_counter){
	const double fraction = static_cast<double>(entry.second)/n_psp_base;
	const int percent = std::round(100*fraction);
	std::cout << std::left << std::setw(25)
	<< (event_type::names[entry.first] + std::string(":"))
	<< entry.second << " (" << percent << "%)\n";

	}
	for(auto const & t: allowed_types){
	if(type_counter[t] < 0.05 * n_psp_base){
	std::cerr << "Less than 5% of the events are of type " << event_type::names[t] << std::endl;
	return EXIT_FAILURE;
	}
	}

	// Wm4j
	proc = Process{{pid::proton,pid::proton}, 4, pid::Wm};
	n_psp = n_psp_base;

	subl_channels = Subleading::all;
	allowed_types = {event_type::FKL,
	event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf,
	event_type::qqxmid};
	type_counter.clear();

	for( size_t i = 0; i<n_psp; ++i){
	const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
	boson_prop, ran};
	if(psp.status()==good){
	const Event ev(to_UnclusteredEvent(psp), jet_para.def, jet_para.min_pt);
	++type_counter[ev.type()];
	if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
	== allowed_types.cend()) {
	bail_out(psp, "Found not allowed event of type "
	+std::string(event_type::names[ev.type()]));
	}
	} else { // bad process -> try again
	++n_psp;
	}
	}
	std::cout << "Wm+4j: Took " << n_psp << " to generate "
	<< n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl;
	std::cout << "States by classification:\n";
	for(auto const & entry: type_counter){
	const double fraction = static_cast<double>(entry.second)/n_psp_base;
	const int percent = std::round(100*fraction);
	std::cout << std::left << std::setw(25)
	<< (event_type::names[entry.first] + std::string(":"))
	<< entry.second << " (" << percent << "%)\n";

	}
	for(auto const & t: allowed_types){
	if(type_counter[t] < 0.03 * n_psp_base){
	std::cerr << "Less than 3% of the events are of type " << event_type::names[t] << std::endl;
	return EXIT_FAILURE;
	}
	}

	std::cout << "All processes passed." << std::endl;
	return EXIT_SUCCESS;
	}
	diff --git a/doc/sphinx/HEJFOG.rst b/doc/sphinx/HEJFOG.rst
	index 42c8af6..02c1e89 100644
	--- a/doc/sphinx/HEJFOG.rst
	+++ b/doc/sphinx/HEJFOG.rst
	@@ -1,281 +1,291 @@
	The HEJ Fixed Order Generator
	=============================

	For high jet multiplicities event generation with standard fixed-order
	generators becomes increasingly cumbersome. For example, the
	leading-order production of a Higgs Boson with five or more jets is
	computationally prohibitively expensive.

	To this end, HEJ 2 provides the ``HEJFOG`` fixed-order generator
	that allows to generate events with high jet multiplicities. To
	facilitate the computation the limit of Multi-Regge Kinematics with
	large invariant masses between all outgoing particles is assumed in the
	matrix elements. The typical use of the ``HEJFOG`` is to supplement
	low-multiplicity events from standard generators with high-multiplicity
	events before using the HEJ 2 program to add high-energy
	resummation.

	Installation
	------------

	The ``HEJFOG`` comes bundled together with HEJ 2 and the
	installation is very similar. After downloading HEJ 2 and
	installing the prerequisites as described in :ref:`Installation` the
	``HEJFOG`` can be installed with::

	cmake /path/to/FixedOrderGen -DCMAKE_INSTALL_PREFIX=target/directory -DCMAKE_BUILD_TYPE=Release
	make install

	where :file:`/path/to/FixedOrderGen` refers to the :file:`FixedOrderGen`
	subdirectory in the HEJ 2 directory. If HEJ 2 was
	installed to a non-standard location, it may be necessary to specify the
	directory containing :file:`HEJ-config.cmake`. If the base installation
	directory is :file:`/path/to/HEJ`, :file:`HEJ-config.cmake` should be
	found in :file:`/path/to/HEJ/lib/cmake/HEJ` and the commands for
	installing the ``HEJFOG`` would read::

	cmake /path/to/FixedOrderGen -DHEJ_DIR=/path/to/HEJ/lib/cmake/HEJ -DCMAKE_INSTALL_PREFIX=target/directory -DCMAKE_BUILD_TYPE=Release
	make install

	The installation can be tested with::

	make test

	provided that the NNPDF 3.0 PDF set is installed.

	Running the fixed-order generator
	---------------------------------

	After installing the ``HEJFOG`` you can modify the provided
	configuration file :file:`configFO.yml` and run the generator with::

	HEJFOG configFO.yml

	The resulting event file, by default named :file:`HEJFO.lhe`, can then be
	fed into HEJ 2 like any event file generated from a standard
	fixed-order generator, see :ref:`Running HEJ 2`.

	Settings
	--------

	Similar to HEJ 2, the ``HEJFOG`` uses a `YAML
	<http://yaml.org/>`_ configuration file. The settings are

	.. _`process`:

	process
	The scattering process for which events are being generated. The
	format is

	:code:`in1 in2 => out1 out2 ...`

	The incoming particles, :code:`in1`, :code:`in2` can be

	- quarks: :code:`u`, :code:`d`, :code:`u_bar`, and so on
	- gluons: :code:`g`
	- protons :code:`p` or antiprotons :code:`p_bar`

	At most one of the outgoing particles can be a boson, the rest has to be
	partonic. At the moment only the Higgs boson :code:`h` is supported. All
	other outgoing particles are jets. Multiple jets can be grouped together, so
	:code:`p p => h j j` is the same as :code:`p p => h 2j`. There have to be at
	least two jets. Further decays of the boson can be added through the
	:ref:`particle properties<particle properties: particle: decays>`.

	.. _`events`:

	events
	Specifies the number of events to generate.

	.. _`jets`:

	jets
	Defines the properties of the generated jets.

	.. _`jets: min pt`:

	min pt
	Minimum jet transverse momentum in GeV.

	.. _`jets: peak pt`:

	peak pt
	Optional setting to specify the dominant jet transverse momentum
	in GeV. If the generated events are used as input for HEJ
	resummation, this should be set to the minimum transverse momentum
	of the resummation jets. The effect is that only a small fraction
	of jets will be generated with a transverse momentum below the
	value of this setting.

	.. _`jets: algorithm`:

	algorithm
	The algorithm used to define jets. Allowed settings are
	:code:`kt`, :code:`cambridge`, :code:`antikt`,
	:code:`cambridge for passive`. See the `FastJet
	<http://fastjet.fr/>`_ documentation for a description of these
	algorithms.

	.. _`jets: R`:

	R
	The R parameter used in the jet algorithm.

	.. _`jets: max rapidity`:

	max rapidity
	Maximum absolute value of the jet rapidity.

	.. _`beam`:

	beam
	Defines various properties of the collider beam.

	.. _`beam: energy`:

	energy
	The beam energy in GeV. For example, the 13
	TeV LHC corresponds to a value of 6500.

	.. _`beam: particles`:

	particles
	A list :code:`[p1, p2]` of two beam particles. The only supported
	entries are protons :code:`p` and antiprotons :code:`p_bar`.

	.. _`pdf`:

	pdf
	The `LHAPDF number <https://lhapdf.hepforge.org/pdfsets>`_ of the PDF set.
	For example, 230000 corresponds to an NNPDF 2.3 NLO PDF set.

	.. _`subleading fraction`:

	subleading fraction
	This setting is related to the fraction of events, that are not a FKL
	- configuration. Currently only unordered emissions are implemented, and only
	- for Higgs plus multijet processes.
	- Typically, this value should be between 0.01 and 0.1.
	+ configuration. All possible subleading process are listed in
	+ :ref:`subleading channels<subleading channels>`. Typically, this value
	+ should be between 0.01 and 0.1.

	.. _`subleading channels`:

	subleading channels
	- Optional parameters to select a specific unordered process, multiple
	- channels can be selected at once. Only matters if :code:`subleading fraction`
	- is non-zero. Currently three values are allowed:
	+ Optional parameters to select a specific subleading process, multiple
	+ channels can be selected at once. If multiple subleading configurations
	+ are possible one will be selected at random for each event, thus each
	+ final state will include at most one subleading process, e.g. either
	+ :code:`unordered` or :code:`qqx`. The following values are allowed:

	- :code:`all`: All channels allowed, default if nothing else set
	- :code:`none`: No subleading contribution, only FKL allowed
	Equivalent to :code:`subleading fraction: 0`
	- :code:`unordered`: Unordered emission allowed
	Unordered emission are any rapidity ordering, where exactly one gluon is
	emitted outside the rapidity ordering required in FKL events. More
	precisely, if at least one of the incoming particles is a quark or
	antiquark and there are more than two jets in the final state,
	:code:`subleading fraction` states the probability that the flavours
	of the outgoing particles are assigned in such a way that an unordered
	- configuration arises.
	+ configuration arises. Unordered emissions are currently not implemented
	+ for pure jet production.
	+ - :code:`qqx`: Production of quark-antiquark pair
	+ If allowed processes with at least three jets might include a
	+ quark-antiquark pair, i.e. one quark (in rapidity) next to one antiquark.
	+ The pair can be at any position in the FKL chain. Similar to above
	+ :code:`subleading fraction` gives the probability of turning two gluons
	+ into the quark-antiquark pair. Quark-antiquark pairs are currently only
	+ implemented for W boson production.

	.. _`unweight`:

	unweight
	This setting defines the parameters for the partial unweighting of
	events. You can disable unweighting by removing this entry from the
	configuration file.

	.. _`unweight: sample size`:

	sample size
	The number of weighted events used to calibrate the unweighting.
	A good default is to set this to the number of target
	`events`_. If the number of `events`_ is large this can
	lead to significant memory consumption and a lower value should be
	chosen. Contrarily, for large multiplicities the unweighting
	efficiency becomes worse and the sample size should be increased.

	.. _`unweight: max deviation`:

	max deviation
	Controls the range of events to which unweighting is applied. A
	larger value means that a larger fraction of events are unweighted.
	Typical values are between -1 and 1.

	.. _`particle properties`:

	particle properties
	Specifies various properties of the different particles (Higgs, W or Z).
	This is only relevant if the chosen `process`_ is the production of the
	corresponding particles with jets. E.g. for the `process`_
	:code:`p p => h 2j` the :code:`mass`, :code:`width` and (optionally)
	:code:`decays` of the :code:`Higgs` boson are required, while all other
	particle properties will be ignored.

	.. _`particle properties: particle`:

	Higgs, Wp, Wm or Z
	Name of the boson. Can be any of the once above.

	.. _`particle properties: particle: mass`:

	mass
	The mass of the particle in GeV.

	.. _`particle properties: particle: width`:

	width
	The total decay width of the particle in GeV.

	.. _`particle properties: particle: decays`:

	decays
	Optional setting specifying the decays of the particle. Only the decay
	into two particles is implemented. Each decay has the form

	:code:`{into: [p1,p2], branching ratio: r}`

	where :code:`p1` and :code:`p2` are the particle names of the
	decay product (e.g. :code:`photon`) and :code:`r` is the branching
	ratio.

	Decays of a Higgs boson are treated as the production and subsequent
	decay of an on-shell Higgs boson, so decays into e.g. Z bosons are not
	supported.

	.. _`scales`:

	scales
	Specifies the renormalisation and factorisation scales for the output
	events. For details, see the corresponding entry in the HEJ 2
	:ref:`HEJ 2 settings`. Note that this should usually be a
	single value, as the weights resulting from additional scale choices
	will simply be ignored by HEJ 2.

	.. _`event output`:

	event output
	Specifies the name of a single event output file or a list of such
	files. See the corresponding entry in the HEJ 2
	:ref:`HEJ 2 settings` for details.

	.. _`RanLux init`:

	.. _`random generator`:

	random generator
	Sets parameters for random number generation. See the HEJ 2
	:ref:`HEJ 2 settings` for details.

	.. _`analysis`:

	analysis
	Specifies the name and settings for a custom analysis library. This
	can be useful to specify cuts at the fixed-order level. See the
	corresponding entry in the HEJ 2 :ref:`HEJ 2 settings`
	for details.

	.. _`Higgs coupling`:

	Higgs coupling
	This collects a number of settings concerning the effective coupling
	of the Higgs boson to gluons. See the corresponding entry in the
	HEJ 2 :ref:`HEJ 2 settings` for details.

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 10:25 PM (1 d, 14 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4242639
Default Alt Text: (42 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions