No OneTemporary
Actions

Size

40 KB

Subscribers

None

View Options

	diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
	index 17afb92..b339d20 100644
	--- a/FixedOrderGen/src/PhaseSpacePoint.cc
	+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
	@@ -1,945 +1,950 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "PhaseSpacePoint.hh"

	#include <algorithm>
	#include <cassert>
	#include <cmath>
	#include <cstdlib>
	#include <iostream>
	#include <iterator>
	#include <limits>
	#include <tuple>
	#include <type_traits>
	#include <utility>

	#include "fastjet/ClusterSequence.hh"

	#include "HEJ/Constants.hh"
	#include "HEJ/EWConstants.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/kinematics.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/PDF.hh"
	#include "HEJ/RNG.hh"
	#include "HEJ/utility.hh"

	#include "JetParameters.hh"
	#include "Process.hh"

	namespace {
	static_assert(
	std::numeric_limits<double>::has_quiet_NaN,
	"no quiet NaN for double"
	);
	constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
	} // namespace anonymous

	namespace HEJFOG {
	HEJ::Event::EventData to_EventData(PhaseSpacePoint psp){
	//! @TODO Same function already in HEJ
	HEJ::Event::EventData result;
	result.incoming = std::move(psp).incoming_;
	assert(result.incoming.size() == 2);
	result.outgoing = std::move(psp).outgoing_;
	// technically Event::EventData doesn't have to be sorted,
	// but PhaseSpacePoint should be anyway
	assert(
	std::is_sorted(
	begin(result.outgoing), end(result.outgoing),
	HEJ::rapidity_less{}
	)
	);
	assert(result.outgoing.size() >= 2);
	result.decays = std::move(psp).decays_;
	result.parameters.central = {NaN, NaN, psp.weight()};
	return result;
	}

	PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::begin_partons() const {
	return cbegin_partons();
	}
	PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::cbegin_partons() const {
	return {HEJ::is_parton, cbegin(outgoing()), cend(outgoing())};
	}

	PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::end_partons() const {
	return cend_partons();
	}
	PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::cend_partons() const {
	return {HEJ::is_parton, cend(outgoing()), cend(outgoing())};
	}

	PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::rbegin_partons() const {
	return crbegin_partons();
	}
	PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::crbegin_partons() const {
	return std::reverse_iterator<ConstPartonIterator>( cend_partons() );
	}

	PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::rend_partons() const {
	return crend_partons();
	}
	PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::crend_partons() const {
	return std::reverse_iterator<ConstPartonIterator>( cbegin_partons() );
	}

	PhaseSpacePoint::PartonIterator PhaseSpacePoint::begin_partons() {
	return {HEJ::is_parton, begin(outgoing_), end(outgoing_)};
	}

	PhaseSpacePoint::PartonIterator PhaseSpacePoint::end_partons() {
	return {HEJ::is_parton, end(outgoing_), end(outgoing_)};
	}

	PhaseSpacePoint::ReversePartonIterator PhaseSpacePoint::rbegin_partons() {
	return std::reverse_iterator<PartonIterator>( end_partons() );
	}

	PhaseSpacePoint::ReversePartonIterator PhaseSpacePoint::rend_partons() {
	return std::reverse_iterator<PartonIterator>( begin_partons() );
	}

	namespace {
	bool can_swap_to_uno(
	HEJ::Particle const & p1, HEJ::Particle const & p2
	){
	assert(is_parton(p1) && is_parton(p2));
	return p1.type != HEJ::pid::gluon
	&& p2.type == HEJ::pid::gluon;
	}

	size_t count_gluons(PhaseSpacePoint::ConstPartonIterator first,
	PhaseSpacePoint::ConstPartonIterator last
	){
	return std::count_if(first, last, [](HEJ::Particle const & p)
	{return p.type == HEJ::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
	*/
	Subleading possible_qqx(
	PhaseSpacePoint::ConstPartonIterator first,
	PhaseSpacePoint::ConstReversePartonIterator last
	){
	using namespace subleading;
	assert( std::distance( first,last.base() )>2 );
	Subleading channels = ALL;
	channels.reset(eqqx);
	channels.reset(cqqx);
	auto const ngluon = count_gluons(first,last.base());
	if(ngluon < 2) return channels;
	if(first->type==HEJ::pid::gluon \|\| last->type==HEJ::pid::gluon){
	channels.set(eqqx);
	}
	if(std::distance(first,last.base())>=4){
	channels.set(cqqx);
	}
	return channels;
	}

	bool is_AWZ_proccess(Process const & proc){
	return proc.boson && HEJ::is_AWZ_boson(*proc.boson);
	}

	bool is_up_type(HEJ::Particle const & part){
	return is_anyquark(part) && !(std::abs(part.type)%2);
	}
	bool is_down_type(HEJ::Particle const & part){
	return is_anyquark(part) && std::abs(part.type)%2;
	}
	bool can_couple_to_W(
	HEJ::Particle const & part, HEJ::pid::ParticleID const W_id
	){
	const int W_charge = W_id>0?1:-1;
	return std::abs(part.type)<5
	&& ( (W_charge*part.type > 0 && is_up_type(part))
	\|\| (W_charge*part.type < 0 && is_down_type(part)) );
	}
	}

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

	// decide what kind of subleading process is allowed
	bool can_be_uno_backward = can_swap_to_uno(
	cbegin_partons(), (++cbegin_partons()) );
	bool can_be_uno_forward = can_swap_to_uno(
	crbegin_partons(), (++crbegin_partons()) );
	// Special case: Higgs can not be outside of uno
	if(proc.boson && *proc.boson==pid::Higgs){
	if(outgoing_.begin()->type == pid::Higgs
	\|\| (++outgoing_.begin())->type==pid::Higgs){
	can_be_uno_backward = false;
	}
	if(outgoing_.rbegin()->type == pid::Higgs
	\|\| (++outgoing_.rbegin())->type==pid::Higgs){
	can_be_uno_forward = false;
	}
	}
	if(channels[subleading::uno]){
	channels.set(subleading::uno, can_be_uno_backward \|\| can_be_uno_forward);
	}
	channels &= possible_qqx(cbegin_partons(), crbegin_partons());

	bool allow_strange = true;
	if(is_AWZ_proccess(proc)) {
	if(std::none_of(cbegin_partons(), cend_partons(),
	[&proc](Particle const & p){
	return can_couple_to_W(p, *proc.boson);})) {
	// enforce qqx if A/W/Z can't couple somewhere else
	// this is ensured to work through filter_partons in reconstruct_incoming
	channels.reset(subleading::uno);
	assert(channels.any());
	chance = 1.;
	// strange not allowed for W
	if(std::abs(*proc.boson)== pid::Wp) allow_strange = false;
	}
	}

	std::size_t const nchannels = channels.count();
	// no subleading
	if(nchannels==0) return;
	if(ran.flat() >= chance){
	weight_ /= 1 - chance;
	return;
	}
	weight_ /= chance;

	// select channel
	double const step = 1./nchannels;
	weight_*=nchannels;
	unsigned selected = subleading::first;
	double rnd = nchannels>1?ran.flat():0.;
	for(; selected<=subleading::last; ++selected){
	assert(rnd>=0);
	if(channels[selected]){
	if(rnd<step) break;
	rnd-=step;
	}
	}
	switch(selected){
	case subleading::uno:
	return turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
	case subleading::cqqx:
	return turn_to_cqqx(allow_strange, ran);
	case subleading::eqqx:
	return turn_to_eqqx(allow_strange, ran);
	default:
	throw std::logic_error{"unreachable"};
	}
	}

	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;
	if(can_be_uno_backward && can_be_uno_forward){
	weight_ *= 2.;
	if(ran.flat() < 0.5){
	return std::swap(begin_partons()->type, (++begin_partons())->type);
	}
	return std::swap(rbegin_partons()->type, (++rbegin_partons())->type);
	}
	if(can_be_uno_backward){
	return std::swap(begin_partons()->type, (++begin_partons())->type);
	}
	assert(can_be_uno_forward);
	std::swap(rbegin_partons()->type, (++rbegin_partons())->type);
	}

	//! select flavour of quark
	HEJ::ParticleID PhaseSpacePoint::select_qqx_flavour(
	const bool allow_strange, HEJ::RNG & ran
	){
	const double r1 = 2.*ran.flat()-1.;
	const double max_flavour = allow_strange?HEJ::N_F:HEJ::N_F-1;
	weight_ = max_flavour2;
	double const flavour = HEJ::pid::down + std::floor(std::abs(r1)*max_flavour);
	return static_cast<HEJ::ParticleID>(flavour*(r1<0.?-1:1));
	}

	void PhaseSpacePoint::turn_to_cqqx(const bool allow_strange, HEJ::RNG & ran){
	// we assume all FKL partons to be gluons
	auto first = ++begin_partons();
	auto last = ++rbegin_partons();
	auto const ng = std::distance(first, last.base());
	if(ng < 2)
	throw std::logic_error("not enough gluons to create qqx");
	auto flavour = select_qqx_flavour(allow_strange, ran);
	// select gluon for switch
	if(ng!=2){
	const double steps = 1./(ng-1);
	weight_ /= steps;
	for(auto rnd = ran.flat(); rnd>steps; ++first){
	rnd-=steps;
	}
	}
	first->type = flavour;
	(++first)->type = anti(flavour);
	}

	void PhaseSpacePoint::turn_to_eqqx(const bool allow_strange, HEJ::RNG & ran){
	/// find first and last gluon in FKL chain
	auto first = begin_partons();
	const bool can_forward = !is_anyquark(*first);
	auto last = rbegin_partons();
	const bool can_backward = !is_anyquark(*last);
	if(std::distance(first, last.base()) < 2)
	throw std::logic_error("not enough gluons to create qqx");
	auto flavour = select_qqx_flavour(allow_strange, ran);
	// select gluon for switch
	if(can_forward && !can_backward){
	first->type = flavour;
	(++first)->type = anti(flavour);
	return;
	}
	if(!can_forward && can_backward){
	last->type = flavour;
	(++last)->type = anti(flavour);
	return;
	}
	assert(can_forward && can_backward);
	weight_*=2.;
	if(ran.flat()>0.5){
	first->type = flavour;
	(++first)->type = anti(flavour);
	return;
	}
	last->type = flavour;
	(++last)->type = anti(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 = std::sqrt(pt[0]pt[0]+pt[1]pt[1]+mass_square);
	const double pz=mperp*std::sinh(y);
	const double E=mperp*std::cosh(y);

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

	namespace {
	//! adds a particle to target (in correct rapidity ordering)
	//! @returns positon of insertion
	auto insert_particle(std::vector<HEJ::Particle> & target,
	HEJ::Particle && particle
	){
	const auto pos = std::upper_bound(
	begin(target),end(target),particle,HEJ::rapidity_less{}
	);
	target.insert(pos, std::move(particle));
	return pos;
	}
	}

	PhaseSpacePoint::PhaseSpacePoint(
	Process const & proc,
	JetParameters const & jet_param,
	HEJ::PDF & pdf, double E_beam,
	double const subl_chance,
	Subleading subl_channels,
	ParticlesDecayMap const & particle_decays,
	HEJ::EWConstants const & ew_parameters,
	HEJ::RNG & ran
	){
	assert(proc.njets >= 2);
	status_ = Status::good;
	weight_ = 1;

	// ensure that all setting are consistent
	if(subl_chance == 0.)
	subl_channels.reset();

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

	// generate parton momenta
	const bool is_pure_jets = (nout == proc.njets);
	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(HEJ::Particle{HEJ::pid::gluon, std::move(p_out), {}});
	}
	if(status_ != Status::good) return;

	if(proc.boson){ // decay boson
	auto const & boson_prop = ew_parameters.prop(*proc.boson) ;
	auto boson{ gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran) };
	const auto pos{insert_particle(outgoing_, std::move(boson))};
	const size_t boson_idx = std::distance(begin(outgoing_), pos);
	auto const & boson_decay = particle_decays.find(*proc.boson);
	if( !proc.boson_decay.empty() ){ // decay given in proc
	decays_.emplace(
	boson_idx,
	decay_boson(outgoing_[boson_idx], proc.boson_decay, ran)
	);
	} else if( boson_decay != particle_decays.end()
	&& !boson_decay->second.empty() ){ // decay given explicitly
	decays_.emplace(
	boson_idx,
	decay_boson(outgoing_[boson_idx], boson_decay->second, ran)
	);
	}
	}
	// normalisation of momentum-conserving delta function
	weight_ = std::pow(2M_PI, 4);

	/** @TODO
	* uf (jet_param.min_pt) doesn'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
	*/
	reconstruct_incoming(proc, subl_chance, subl_channels, pdf, E_beam, jet_param.min_pt, ran);
	if(status_ != Status::good) return;
	// set outgoing states
	begin_partons()->type = incoming_[0].type;
	rbegin_partons()->type = incoming_[1].type;

	maybe_turn_to_subl(subl_chance, subl_channels, proc, ran);

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

	// pt generation, see eq:pt_sampling in developer manual
	double PhaseSpacePoint::gen_hard_pt(
	const int np , const double ptmin, const double ptmax, const double /* y */,
	HEJ::RNG & ran
	){
	// heuristic parameter for pt sampling, see eq:pt_par in developer manual
	const double ptpar = ptmin + np/5.;

	const double arctan = std::atan((ptmax - ptmin)/ptpar);
	const double xpt = ran.flat();
	const double pt = ptmin + ptparstd::tan(xptarctan);
	const double cosine = std::cos(xpt*arctan);
	weight_ = ptptpararctan/(cosinecosine);
	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_ = 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_ /= std::pow(16.*std::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_ = Status::not_enough_jets;
	return {};
	}

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

	return partons;
	}

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

	// Generate a y Gaussian distributed around 0
	/// @TODO check magic numbers for different boson Higgs
	/// @TODO better sampling for W
	const double stddev_y = 1.6;
	const double y = random_normal(stddev_y, ran);
	const double r1 = ran.flat();
	const double s_boson = mass*(
	mass + widthstd::tan(M_PI/2.r1 + (r1-1.)*std::atan(mass/width))
	);
	// off-shell s_boson sampling, compensates for Breit-Wigner
	/// @TODO use a flag instead
	if(std::abs(bosonid) == HEJ::pid::Wp){
	weight_/=M_PIM_PI16.; //Corrects B-W factors, see git issue 132
	weight_= masswidth( M_PI+2.std::atan(mass/width) )
	/ ( 1. + std::cos( M_PIr1 + 2.(r1-1.)*std::atan(mass/width) ) );
	}

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

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

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

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

	using part_mask = std::bitset<11>; //!< Selection mask for partons

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

	/// decides which "index" (see index_to_pid) are allowed for process
	part_mask allowed_quarks(HEJ::ParticleID const boson){
	if(std::abs(boson) != HEJ::pid::Wp){
	return ~1; // not a gluon
	}
	// special case W:
	// Wp: anti-down or up-type quark, no b/t
	// Wm: down or anti-up-type quark, no b/t
	return boson>0?0b00011001100
	:0b00100110010;
	}
	}

	- std::array<std::bitset<11>,2> PhaseSpacePoint::incoming_AWZ(
	+ std::array<part_mask,2> PhaseSpacePoint::incoming_AWZ(
	Process const & proc, Subleading const subl_channels,
	- std::array<std::bitset<11>,2> allowed_partons,
	+ std::array<part_mask,2> allowed_partons,
	HEJ::RNG & ran
	){
	assert(proc.boson);
	auto couple_parton = allowed_quarks(*proc.boson);
	- if(
	- // coupling possible through input
	- allowed_partons[0] != (couple_parton&allowed_partons[0])
	- && allowed_partons[1] != (couple_parton&allowed_partons[1])
	- // qqx not possible
	- && (proc.njets < 4
	- \|\| (!subl_channels[subleading::eqqx]
	- && !subl_channels[subleading::cqqx])
	- )
	- )
	- {
	- // eqqx only possible if one incoming is a gluon
	- if(proc.njets >= 3 && subl_channels[subleading::eqqx]){
	- couple_parton.set(pid_to_index(HEJ::ParticleID::gluon));
	+ if( // coupling possible through input
	+ allowed_partons[0] == (couple_parton&allowed_partons[0])
	+ \|\| allowed_partons[1] == (couple_parton&allowed_partons[1])
	+ // cqqx possible
	+ \|\| (proc.njets >= 4 && subl_channels[subleading::cqqx])
	+ ){
	+ return allowed_partons;
	+ }
	+ // eqqx only possible if one incoming is a gluon
	+ if(proc.njets >= 3 && subl_channels[subleading::eqqx]){
	+ couple_parton.set(pid_to_index(HEJ::ParticleID::gluon));
	+ }
	+ // only first can couple
	+ if( (allowed_partons[0]&couple_parton).any()
	+ &&(allowed_partons[1]&couple_parton).none()
	+ ){
	+ return {allowed_partons[0]&couple_parton, allowed_partons[1]};
	+ }
	+ // only second can couple
	+ if( (allowed_partons[0]&couple_parton).none()
	+ && (allowed_partons[1]&couple_parton).any()
	+ ){
	+ return {allowed_partons[0], allowed_partons[1]&couple_parton};
	+ }
	+ // both can couple
	+ if( (allowed_partons[0]&couple_parton).any()
	+ && (allowed_partons[1]&couple_parton).any()
	+ ){
	+ double rnd = ran.flat();
	+ weight_*=3.;
	+ if(rnd<1./3.){
	+ return {
	+ allowed_partons[0] & couple_parton,
	+ allowed_partons[1] & ~couple_parton
	+ };
	}
	- if( (allowed_partons[0]&couple_parton).any()
	- &&(allowed_partons[1]&couple_parton).none()
	- ){
	- allowed_partons[0]&=couple_parton;
	- } else if( (allowed_partons[0]&couple_parton).none()
	- && (allowed_partons[1]&couple_parton).any()
	- ){
	- allowed_partons[1]&=couple_parton;
	- } else if( (allowed_partons[0]&couple_parton).any()
	- && (allowed_partons[1]&couple_parton).any()
	- ){
	- double rnd = ran.flat();
	- weight_*=3.;
	- if(rnd<1./3.){
	- allowed_partons[0] &= couple_parton;
	- allowed_partons[1] &= ~couple_parton;
	- } else if(rnd<2./3.){
	- allowed_partons[0] &= ~couple_parton;
	- allowed_partons[1] &= couple_parton;
	- } else {
	- allowed_partons[0] &= couple_parton;
	- allowed_partons[1] &= couple_parton;
	- }
	- } else {
	- throw std::invalid_argument{"Incoming state not allowed."};
	+ if(rnd<2./3.){
	+ return {
	+ allowed_partons[0] & ~couple_parton,
	+ allowed_partons[1] & couple_parton
	+ };
	}
	+ return {
	+ allowed_partons[0] & couple_parton,
	+ allowed_partons[1] & couple_parton
	+ };
	}
	- return allowed_partons;
	+ throw std::invalid_argument{"Incoming state not allowed."};
	}

	- std::array<std::bitset<11>,2> PhaseSpacePoint::incoming_eqqx(
	- std::array<std::bitset<11>,2> allowed_partons, HEJ::RNG & ran
	+ std::array<part_mask,2> PhaseSpacePoint::incoming_eqqx(
	+ std::array<part_mask,2> allowed_partons, HEJ::RNG & ran
	){
	auto const gluon_idx = pid_to_index(HEJ::pid::gluon);
	auto & first_beam = allowed_partons[0];
	auto & second_beam = allowed_partons[1];
	if(first_beam[gluon_idx] && !second_beam[gluon_idx]){
	first_beam.reset();
	first_beam.set(gluon_idx);
	return allowed_partons;
	}
	if(!first_beam[gluon_idx] && second_beam[gluon_idx]) {
	second_beam.reset();
	second_beam.set(gluon_idx);
	return allowed_partons;
	}
	if(first_beam[gluon_idx] && second_beam[gluon_idx]) {
	// both beams can be gluons
	// if one can't be a quark everything is good
	auto first_quarks = first_beam;
	first_quarks.reset(gluon_idx);
	auto second_quarks = second_beam;
	second_quarks.reset(gluon_idx);
	if(first_quarks.none() \|\| second_quarks.none()){
	return allowed_partons;
	}
	// else choose one to be a gluon
	double rnd = ran.flat();
	weight_*=3.;
	if(rnd<1./3.){
	allowed_partons[0].reset();
	allowed_partons[0].set(gluon_idx);
	allowed_partons[1].reset(gluon_idx);
	} else if(rnd<2./3.){
	allowed_partons[1].reset();
	allowed_partons[1].set(gluon_idx);
	allowed_partons[0].reset(gluon_idx);
	} else {
	allowed_partons[0].reset();
	allowed_partons[0].set(gluon_idx);
	allowed_partons[1].reset();
	allowed_partons[1].set(gluon_idx);
	}
	return allowed_partons;
	}
	throw std::invalid_argument{
	"Incoming state not allowed for pure extremal qqx."};
	}

	- std::array<std::bitset<11>,2> PhaseSpacePoint::incoming_uno(
	- std::array<std::bitset<11>,2> allowed_partons, HEJ::RNG & ran
	+ std::array<part_mask,2> PhaseSpacePoint::incoming_uno(
	+ std::array<part_mask,2> allowed_partons, HEJ::RNG & ran
	){
	auto const gluon_idx = pid_to_index(HEJ::pid::gluon);
	auto & first_beam = allowed_partons[0];
	auto first_quarks = first_beam;
	first_quarks.reset(gluon_idx);
	auto & second_beam = allowed_partons[1];
	auto second_quarks = second_beam;
	second_quarks.reset(gluon_idx);
	if(first_quarks.any() && second_quarks.none()){
	first_beam.reset(gluon_idx);
	return allowed_partons;
	}
	if(first_quarks.none() && second_quarks.any()) {
	second_beam.reset(gluon_idx);
	return allowed_partons;
	}
	if(first_quarks.any() && second_quarks.any()) {
	// both beams can be quarks
	// if one can't be gluon everything is good
	if(!first_beam[gluon_idx] \|\| !second_beam[gluon_idx]){
	return allowed_partons;
	}
	// else choose one to be a quark
	double rnd = ran.flat();
	weight_*=3.;
	if(rnd<1./3.){
	allowed_partons[0].reset(gluon_idx);
	allowed_partons[1].reset();
	allowed_partons[1].set(gluon_idx);
	} else if(rnd<2./3.){
	allowed_partons[1].reset(gluon_idx);
	allowed_partons[0].reset();
	allowed_partons[0].set(gluon_idx);
	} else {
	allowed_partons[0].reset(gluon_idx);
	allowed_partons[1].reset(gluon_idx);
	}
	return allowed_partons;
	}
	throw std::invalid_argument{
	"Incoming state not allowed for pure unordered."};
	}

	/**
	* @brief Returns list of all allowed initial states partons
	*
	* 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
	* 3. pure eqqx requires at least one gluon
	* 4. pure uno requires at least one quark
	*/
	std::array<part_mask,2> PhaseSpacePoint::allowed_incoming(
	Process const & proc,
	double const subl_chance, Subleading const subl_channels,
	HEJ::RNG & ran
	){
	// all possible incoming states
	std::array<part_mask,2> allowed_partons{
	init_allowed(proc.incoming[0]),
	init_allowed(proc.incoming[1])
	};
	// special case A/W/Z
	if(proc.boson && is_AWZ_proccess(proc)){
	allowed_partons = incoming_AWZ(proc, subl_channels, allowed_partons, ran);
	}
	// special case: pure subleading
	if(subl_chance!=1.){
	return allowed_partons;
	}
	auto other_channels = subl_channels;
	// pure eqqx
	other_channels.reset(subleading::eqqx);
	if(other_channels.none()){
	return incoming_eqqx(allowed_partons, ran);
	}
	other_channels = subl_channels;
	// pure uno
	other_channels.reset(subleading::uno);
	if(other_channels.none()){
	return incoming_uno(allowed_partons, ran);
	}
	return allowed_partons;
	}


	void PhaseSpacePoint::reconstruct_incoming(
	Process const & proc,
	double const subl_chance, Subleading 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].E()-incoming_[0].pz())/sqrts;
	const double xb=(incoming_[1].E()+incoming_[1].pz())/sqrts;
	// abort if phase space point is outside of collider energy reach
	if (xa>1. \|\| xb>1.){
	weight_=0;
	status_ = Status::too_much_energy;
	return;
	}
	auto const & ids = proc.incoming;
	std::array<part_mask,2> allowed_partons
	= allowed_incoming(proc, subl_chance, subl_channels, ran);
	for(size_t i = 0; i < 2; ++i){
	if(ids[i] == HEJ::pid::proton \|\| ids[i] == HEJ::pid::p_bar){
	// pick ids according to pdfs
	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, part_mask allowed_partons, HEJ::RNG & ran
	){
	std::array<double,11> pdf_wt;
	pdf_wt[0] = allowed_partons[0]?
	std::fabs(pdf.pdfpt(beam_idx,x,uf,HEJ::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.*std::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(std::abs(boson) != HEJ::pid::Wp) return; // only matters for W
	/// @TODO this could be use to sanity check gamma and Z

	// find all possible quarks
	std::vector<PartonIterator> allowed_parts;
	for(auto part_it=begin_partons(); part_it!=end_partons(); ++part_it){
	// Wp -> up OR anti-down, Wm -> anti-up OR down, no bottom
	if ( can_couple_to_W(*part_it, boson) )
	allowed_parts.push_back(part_it);
	}
	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 = std::floor(ran.flat()*allowed_parts.size());
	weight_ *= allowed_parts.size();
	}
	const int W_charge = boson>0?1:-1;
	allowed_parts[idx]->type =
	static_cast<HEJ::ParticleID>( allowed_parts[idx]->type - W_charge );
	}

	double PhaseSpacePoint::random_normal( double stddev, HEJ::RNG & ran ){
	const double r1 = ran.flat();
	const double r2 = ran.flat();
	const double lninvr1 = -std::log(r1);
	const double result = stddevstd::sqrt(2.lninvr1)std::cos(2.M_PI*r2);
	weight_ = exp(resultresult/(2stddevstddev))std::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 HEJ::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;
	if(decays.size()==1) return decays.front();
	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<HEJ::Particle> PhaseSpacePoint::decay_boson(
	HEJ::Particle const & parent,
	std::vector<Decay> const & decays,
	HEJ::RNG & ran
	){
	const auto channel = select_decay_channel(decays, ran);
	return decay_boson(parent, channel.products, ran);
	}

	std::vector<HEJ::Particle> PhaseSpacePoint::decay_boson(
	HEJ::Particle const & parent,
	std::vector<HEJ::ParticleID> const & decays,
	HEJ::RNG & ran
	){
	if(decays.size() != 2){
	throw HEJ::not_implemented{
	"only decays into two particles are implemented"
	};
	}
	std::vector<HEJ::Particle> decay_products(decays.size());
	for(size_t i = 0; i < decays.size(); ++i){
	decay_products[i].type = decays[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.; // Jacobian Factors for W in line 418
	const double sin_phi = std::sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
	const double px = Estd::cos(theta)sin_phi;
	const double py = Estd::sin(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/CMakeLists.txt b/FixedOrderGen/t/CMakeLists.txt
	index 0cc08e5..ef2d76f 100644
	--- a/FixedOrderGen/t/CMakeLists.txt
	+++ b/FixedOrderGen/t/CMakeLists.txt
	@@ -1,200 +1,205 @@
	set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}")
	set(runcard_dir ${tst_dir}/runcards)

	function(add_long_test)
	set(options NONE)
	set(oneValueArgs NAME)
	set(multiValueArgs COMMAND)
	cmake_parse_arguments(TEST "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN})
	add_test(
	NAME ${TEST_NAME}_short
	COMMAND ${TEST_COMMAND} short
	)
	if(TEST_ALL)
	add_test(
	NAME ${TEST_NAME}
	COMMAND ${TEST_COMMAND}
	)
	endif()
	endfunction()

	foreach(tst W_reconstruct_enu W_2j_classify W_nj_classify)
	add_executable(test_${tst} ${tst_dir}/${tst}.cc)
	target_link_libraries(test_${tst} hejfog_lib)
	add_test(NAME ${tst}
	COMMAND test_${tst}
	WORKING_DIRECTORY ${runcard_dir}
	)
	endforeach()
	# this only tests if the runcard actually works, not if the result is correct
	add_test(
	NAME main_2j
	COMMAND HEJFOG ${runcard_dir}/2j.yml
	)
	add_test(
	NAME main_h2j
	COMMAND HEJFOG ${runcard_dir}/h2j.yml
	)
	add_test(
	NAME main_h2j_decay
	COMMAND HEJFOG ${runcard_dir}/h2j_decay.yml
	)
	add_test(
	NAME peakpt
	COMMAND HEJFOG ${runcard_dir}/2j_peak.yml
	)
	+# special case where coupling not trivial
	+add_test(
	+ NAME main_W4j_uno+eqqx
	+ COMMAND HEJFOG ${runcard_dir}/Wp4j_uno+eqqx.yml
	+)

	# check that uno emission doesn't change FKL xs
	add_executable(FKL_uno FKL_uno.cc)
	target_link_libraries(FKL_uno hejfog_lib)
	add_test(
	NAME FKL_uno
	# calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	COMMAND FKL_uno ${runcard_dir}/h3j_uno.yml 0.0243548 0.000119862
	)
	# xs tests
	add_executable(xs_gen xs_gen.cc)
	target_link_libraries(xs_gen hejfog_lib)
	## Higgs
	add_test(
	NAME xs_hqQ
	# calculated with Sherpa see #132
	COMMAND xs_gen ${runcard_dir}/hqQ.yml 1.612e-02 1.26303e-04
	)
	add_test(
	NAME xs_h2j
	# calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	COMMAND xs_gen ${runcard_dir}/h2j.yml 2.04928 0.00956022
	)
	add_test(
	NAME xs_h3j
	# calculated with HEJ revision bd4388fe55cbc3f5a7b6139096456c551294aa31
	COMMAND xs_gen ${runcard_dir}/h3j.yml 1.07807 0.0132409
	)
	add_long_test(
	NAME xs_h3j_uno
	# calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	COMMAND xs_gen ${runcard_dir}/h3j_uno.yml 0.00347538 3.85875e-05
	)
	add_test(
	NAME xs_h5j
	# calculated with HEJ revision dbde2ffbb3b383ae6709b2424d8f0f9d5658270b
	COMMAND xs_gen ${runcard_dir}/h5j.yml 0.0112504 0.000199633
	)
	add_test(
	NAME xs_h2j_decay
	# calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	COMMAND xs_gen ${runcard_dir}/h2j_decay.yml 0.00466994 2.20995e-05
	)
	## pure jets
	add_test(
	NAME xs_qQ
	# calculated with Sherpa see #132
	COMMAND xs_gen ${runcard_dir}/qQ.yml 7.354e+05 1.905e+03
	)
	add_test(
	NAME xs_2j
	# calculated with "combined" HEJ svn r3480
	COMMAND xs_gen ${runcard_dir}/2j.yml 86.42031848e06 590570
	)
	add_test(
	NAME xs_3j_uno
	# calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1
	COMMAND xs_gen ${runcard_dir}/3j_uno.yml 900 14.3131
	)
	add_test(
	NAME xs_3j_qqx
	# calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1
	COMMAND xs_gen ${runcard_dir}/3j_qqx.yml 62040 1005
	)
	add_test(
	NAME xs_4j_qqx
	# calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1
	COMMAND xs_gen ${runcard_dir}/4j_qqx.yml 28936 550
	)
	add_test(
	NAME xs_4j
	# calculated with HEJ revision 13207b5f67a5f40a2141aa7ee515b022bd4efb65
	COMMAND xs_gen ${runcard_dir}/4j.yml 915072 15402.4
	)
	## W
	add_test(
	NAME xs_WqQ
	# calculated with Sherpa see #132
	COMMAND xs_gen ${runcard_dir}/WpqQ.yml 3.086e+00 4.511e-02
	)
	add_test(
	NAME xs_W2j
	# calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1
	COMMAND xs_gen ${runcard_dir}/Wm2j.yml 4.177443e+02 7.446928e+00
	)
	add_long_test(
	NAME xs_W3j_eqqx
	# calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1
	COMMAND xs_gen ${runcard_dir}/Wp3j_qqx.yml 8.710751e+00 1.245725e-01
	)
	add_long_test(
	NAME xs_W3j_uno
	# calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1
	COMMAND xs_gen ${runcard_dir}/Wp3j_uno.yml 3.000955e-01 5.831799e-03
	)
	add_long_test(
	NAME xs_W4j_qqx
	# calculated with HEJ revision be065dc9a21e5965ce57583f6c0a3d953664b82b
	COMMAND xs_gen ${runcard_dir}/Wp4j_qqx.yml 9.274718e-02 4.875742e-03
	)

	# Test that sum of partons == proton
	add_executable(PSP_channel PSP_channel.cc)
	target_link_libraries(PSP_channel hejfog_lib)
	# pure jets
	add_test(
	NAME channel_2j
	COMMAND PSP_channel ${runcard_dir}/2j.yml
	)
	add_test(
	NAME channel_3j_qqx
	COMMAND PSP_channel ${runcard_dir}/3j_qqx.yml
	)
	add_test(
	NAME channel_3j_uno
	COMMAND PSP_channel ${runcard_dir}/3j_uno.yml
	)
	add_long_test(
	NAME channel_4j_qqx
	COMMAND PSP_channel ${runcard_dir}/4j_qqx.yml
	)
	# W+jets
	# pure jets
	add_test(
	NAME channel_W2j
	COMMAND PSP_channel ${runcard_dir}/Wm2j.yml
	)
	add_test(
	NAME channel_W3j_uno
	COMMAND PSP_channel ${runcard_dir}/Wp3j_uno.yml
	)
	add_test(
	NAME channel_W3j_eqqx
	COMMAND PSP_channel ${runcard_dir}/Wp3j_qqx.yml
	)
	add_long_test(
	NAME channel_W4j_qqx
	COMMAND PSP_channel ${runcard_dir}/Wp4j_qqx.yml
	)

	# Test that each subleading channel is generated consistently
	add_executable(PSP_subl PSP_subl.cc)
	target_link_libraries(PSP_subl hejfog_lib)
	add_long_test(
	NAME subl_5j
	COMMAND PSP_subl ${runcard_dir}/5j.yml
	)
	add_long_test(
	NAME subl_h5j
	COMMAND PSP_subl ${runcard_dir}/h5j.yml
	)
	add_long_test(
	NAME subl_W5j
	COMMAND PSP_subl ${runcard_dir}/Wm5j.yml
	)
	diff --git a/FixedOrderGen/t/runcards/Wp4j_uno+eqqx.yml b/FixedOrderGen/t/runcards/Wp4j_uno+eqqx.yml
	new file mode 100644
	index 0000000..fc223db
	--- /dev/null
	+++ b/FixedOrderGen/t/runcards/Wp4j_uno+eqqx.yml
	@@ -0,0 +1,42 @@
	+events: 10000
	+
	+jets:
	+ min pt: 40
	+ R: 0.4
	+ algorithm: antikt
	+ max rapidity: 6
	+
	+beam:
	+ energy: 6000
	+ particles: [p, p]
	+
	+pdf: 11000
	+
	+# code should enforce coupling at proton leg
	+process: d p => W+ 4j
	+
	+subleading fraction: 1.
	+subleading channels:
	+ - eqqx
	+ - uno
	+
	+scales: 125
	+
	+particle properties:
	+ Higgs:
	+ mass: 125
	+ width: 0
	+ W:
	+ mass: 80.385
	+ width: 2.085
	+ Z:
	+ mass: 91.187
	+ width: 2.495
	+
	+decays:
	+ Wp: {into: [e+, nu_e], branching ratio: 1}
	+
	+random generator:
	+ name: mixmax
	+
	+vev: 246.2196508

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 10:15 AM (1 d, 16 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5111071
Default Alt Text: (40 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions