No OneTemporary
Actions

Size

105 KB

Subscribers

None

View Options

	diff --git a/src/Hjets.cc b/src/Hjets.cc
	index 8d901f2..804896b 100644
	--- a/src/Hjets.cc
	+++ b/src/Hjets.cc
	@@ -1,483 +1,481 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "HEJ/jets.hh"
	#include "HEJ/Hjets.hh"

	#include <array>
	#include <cassert>
	#include <cmath>
	#include <complex>
	#include <limits>

	#include "HEJ/ConfigFlags.hh"
	#include "HEJ/Constants.hh"
	#include "HEJ/LorentzVector.hh"
	#include "HEJ/utility.hh"

	// generated headers
	#include "HEJ/currents/j_h_j.hh"
	#include "HEJ/currents/jgh_j.hh"
	#include "HEJ/currents/juno_h_j.hh"
	#include "HEJ/currents/juno_jgh.hh"

	#ifdef HEJ_BUILD_WITH_QCDLOOP

	#include "qcdloop/qcdloop.h"

	#else

	#include "HEJ/exceptions.hh"

	#endif

	namespace HEJ {
	namespace currents {

	namespace {
	using COM = std::complex<double>;

	constexpr double infinity = std::numeric_limits<double>::infinity(); // NOLINT

	// Loop integrals
	#ifdef HEJ_BUILD_WITH_QCDLOOP
	const COM LOOPRWFACTOR = (COM(0.,1.)M_PIM_PI)/std::pow((2.*M_PI),4);

	COM B0DD(HLV const & q, double mq)
	{
	static std::vector<std::complex<double>> result(3);
	static auto ql_B0 = [](){
	ql::Bubble<std::complex<double>,double,double> ql_B0;
	ql_B0.setCacheSize(100);
	return ql_B0;
	}();
	static std::vector<double> masses(2);
	static std::vector<double> momenta(1);
	for(auto & m: masses) m = mq*mq;
	momenta.front() = q.m2();
	ql_B0.integral(result, 1, masses, momenta);
	return result[0];
	}
	COM C0DD(HLV const & q1, HLV const & q2, double mq)
	{
	static std::vector<std::complex<double>> result(3);
	static auto ql_C0 = [](){
	ql::Triangle<std::complex<double>,double,double> ql_C0;
	ql_C0.setCacheSize(100);
	return ql_C0;
	}();
	static std::vector<double> masses(3);
	static std::vector<double> momenta(3);
	for(auto & m: masses) m = mq*mq;
	momenta[0] = q1.m2();
	momenta[1] = q2.m2();
	momenta[2] = (q1+q2).m2();
	ql_C0.integral(result, 1, masses, momenta);
	return result[0];
	}

	// Kallen lambda functions, see eq:lambda in developer manual
	double lambda(const double s1, const double s2, const double s3) {
	return s1s1 + s2s2 + s3s3 - 2s1s2 - 2s1s3 - 2s2*s3;
	}

	// eq: T_1 in developer manual
	COM T1(HLV const & q1, HLV const & q2, const double m) {
	const double q12 = q1.m2();
	const double q22 = q2.m2();
	const HLV ph = q1 - q2;
	const double ph2 = ph.m2();

	const double lam = lambda(q12, q22, ph2);

	assert(m > 0.);

	const double m2 = m*m;

	return
	- C0DD(q1, -q2, m)(2.m2 + 1./2.(q12 + q22 - ph2) + 2.q12q22ph2/lam)
	- (B0DD(q2, m) - B0DD(ph, m))(q22 - q12 - ph2)q22/lam
	- (B0DD(q1, m) - B0DD(ph, m))(q12 - q22 - ph2)q12/lam
	- 1.;
	}

	// eq: T_2 in developer manual
	COM T2(HLV const & q1, HLV const & q2, const double m) {
	const double q12 = q1.m2();
	const double q22 = q2.m2();
	const HLV ph = q1 - q2;
	const double ph2 = ph.m2();

	const double lam = lambda(q12, q22, ph2);

	assert(m > 0.);

	const double m2 = m*m;

	return
	C0DD(q1, -q2, m)*(
	4.m2/lam(ph2 - q12 - q22) - 1. - 4.q12q22/lam*(
	1 + 3.ph2(q12 + q22 - ph2)/lam
	)
	)
	- (B0DD(q2, m) - B0DD(ph, m))(1. + 6.q12/lam(q22 - q12 + ph2))2.*q22/lam
	- (B0DD(q1, m) - B0DD(ph, m))(1. + 6.q22/lam(q12 - q22 + ph2))2.*q12/lam
	- 2.*(q12 + q22 - ph2)/lam;
	}

	#else // no QCDloop

	COM T1(HLV const & /q1/, HLV const & /q2/, double /mt/){
	throw std::logic_error{"T1 called without QCDloop support"};
	}

	COM T2(HLV const & /q1/, HLV const & /q2/, double /mt/){
	throw std::logic_error{"T2 called without QCDloop support"};
	}

	#endif

	// prefactors of g^{\mu \nu} and q_2^\mu q_1^\nu in Higgs boson emission vertex
	// see eq:VH in developer manual, but without global factor \alpha_s
	std::array<COM, 2> TT(
	HLV const & qH1, HLV const & qH2,
	const double mt, const bool inc_bottom,
	const double mb, const double vev
	) {
	if(mt == infinity) {
	std::array<COM, 2> res = {qH1.dot(qH2), 1.};
	for(auto & tt: res) tt /= (3.M_PIvev);
	return res;
	}
	std::array<COM, 2> res = {T1(qH1, qH2, mt), T2(qH1, qH2, mt)};
	for(auto & tt: res) tt = mtmt;
	if(inc_bottom) {
	res[0] += mbmbT1(qH1, qH2, mb);
	res[1] += mbmbT2(qH1, qH2, mb);
	}
	for(auto & tt: res) tt /= M_PI*vev;
	return res;
	}

	/**
	* @brief Higgs+Jets FKL Contributions, function to handle all incoming types.
	* @param p1out Outgoing Particle 1
	* @param p1in Incoming Particle 1
	* @param p2out Outgoing Particle 2
	* @param p2in Incoming Particle 2
	* @param qH1 t-channel momenta into higgs vertex
	* @param qH2 t-channel momenta out of higgs vertex
	* @param mt top mass (inf or value)
	* @param inc_bottom whether to include bottom mass effects (true) or not (false)
	* @param mb bottom mass (value)
	* @param vev Higgs vacuum expectation value
	*
	* Calculates j^\mu H j_\mu. FKL with higgs vertex somewhere in the FKL chain.
	* Handles all possible incoming states.
	*/

	// }

	} // namespace

	double ME_H_qQ(
	HLV const & p1out, HLV const & p1in,
	HLV const & p2out, HLV const & p2in,
	HLV const & qH1, HLV const & qH2,
	const double mt, const bool inc_bottom, const double mb, const double vev

	){
	using helicity::plus;
	using helicity::minus;

	const auto qqH1 = split_into_lightlike(qH1);
	const HLV qH11 = qqH1.first;
	const HLV qH12 = -qqH1.second;
	const auto qqH2 = split_into_lightlike(qH2);
	const HLV qH21 = qqH2.first;
	const HLV qH22 = -qqH2.second;
	// since qH1.m2(), qH2.m2() < 0 the following assertions are always true
	assert(qH11.e() > 0);
	assert(qH12.e() > 0);
	assert(qH21.e() > 0);
	assert(qH22.e() > 0);

	const auto T_pref = TT(qH1, qH2, mt, inc_bottom, mb, vev);

	const COM amp_mm = HEJ::j_h_j<minus, minus>(
	p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1]
	);
	const COM amp_mp = HEJ::j_h_j<minus, plus>(
	p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1]
	);
	const COM amp_pm = HEJ::j_h_j<plus, minus>(
	p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1]
	);
	const COM amp_pp = HEJ::j_h_j<plus, plus>(
	p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1]
	);

	// square of amplitudes, averaged over helicities
	- const double amp2 = (norm(amp_mm) + norm(amp_mp) + norm(amp_pm) + norm(amp_pp));
	-
	- return amp2/((p1in-p1out).m2()(p2in-p2out).m2()qH1.m2()*qH2.m2());
	+ return (norm(amp_mm) + norm(amp_mp) + norm(amp_pm) + norm(amp_pp));
	}

	//@}

	double ME_jgH_j(
	HLV const & ph, HLV const & pa,
	HLV const & pn, HLV const & pb,
	const double mt, const bool inc_bottom, const double mb, const double vev
	){
	using helicity::plus;
	using helicity::minus;

	const auto pH = split_into_lightlike(ph);
	const HLV ph1 = pH.first;
	const HLV ph2 = pH.second;
	// since pH.m2() > 0 the following assertions are always true
	assert(ph1.e() > 0);
	assert(ph2.e() > 0);

	const auto T_pref = TT(pa, pa-ph, mt, inc_bottom, mb, vev);

	// only distinguish between same and different helicities,
	// the other two combinations just add a factor of 2
	const COM amp_mm = HEJ::jgh_j<minus, minus>(
	pa, pb, pn, ph1, ph2, T_pref[0], T_pref[1]
	);
	const COM amp_mp = HEJ::jgh_j<minus, plus>(
	pa, pb, pn, ph1, ph2, T_pref[0], T_pref[1]
	);
	constexpr double hel_factor = 2.;

	// sum over squares of helicity amplitudes
	return hel_factor*(norm(amp_mm) + norm(amp_mp));
	}

	namespace {

	template<Helicity h1, Helicity h2, Helicity hg>
	double amp_juno_jgh(
	HLV const & pg, HLV const & p1, HLV const & pa,
	HLV const & ph1, HLV const & ph2, HLV const & pb,
	std::array<COM, 2> const & T_pref
	) {
	// TODO: code duplication with Wjets and pure jets
	const COM u1 = U1_jgh<h1, h2, hg>(pg, p1, pa, ph1, ph2, pb, T_pref[0], T_pref[1]);
	const COM u2 = U2_jgh<h1, h2, hg>(pg, p1, pa, ph1, ph2, pb, T_pref[0], T_pref[1]);
	const COM l = L_jgh<h1, h2, hg>(pg, p1, pa, ph1, ph2, pb, T_pref[0], T_pref[1]);
	return C_Fstd::norm(u1+u2) - C_Astd::real((u1-l)*std::conj(u2+l));
	}

	} // namespace

	double ME_juno_jgH(
	HLV const & pg,
	HLV const & p1, HLV const & pa,
	HLV const & ph, HLV const & pb,
	const double mt, const bool inc_bottom, const double mb, const double vev
	) {
	using Helicity::plus;
	using Helicity::minus;

	const auto T_pref = TT(pb, pb-ph, mt, inc_bottom, mb, vev);

	const auto pH = split_into_lightlike(ph);
	const HLV ph1 = pH.first;
	const HLV ph2 = pH.second;
	// since pH.m2() > 0 the following assertions are always true
	assert(ph1.e() > 0);
	assert(ph2.e() > 0);

	// only 4 out of the 8 helicity amplitudes are independent
	// we still compute all of them for better numerical stability (mirror check)
	MultiArray<double, 2, 2, 2> amp;

	// NOLINTNEXTLINE
	#define ASSIGN_HEL(RES, J, H1, H2, HG) \
	RES[H1][H2][HG] = J<H1, H2, HG>( \
	pg, p1, pa, ph1, ph2, pb, T_pref \
	)

	ASSIGN_HEL(amp, amp_juno_jgh, minus, minus, minus);
	ASSIGN_HEL(amp, amp_juno_jgh, minus, minus, plus);
	ASSIGN_HEL(amp, amp_juno_jgh, minus, plus, minus);
	ASSIGN_HEL(amp, amp_juno_jgh, minus, plus, plus);
	ASSIGN_HEL(amp, amp_juno_jgh, plus, minus, minus);
	ASSIGN_HEL(amp, amp_juno_jgh, plus, minus, plus);
	ASSIGN_HEL(amp, amp_juno_jgh, plus, plus, minus);
	ASSIGN_HEL(amp, amp_juno_jgh, plus, plus, plus);

	#undef ASSIGN_HEL

	double ampsq = 0.;
	for(Helicity h1: {minus, plus}) {
	for(Helicity h2: {minus, plus}) {
	for(Helicity hg: {minus, plus}) {
	ampsq += amp[h1][h2][hg];
	}
	}
	}

	return ampsq;
	}


	namespace {

	template<Helicity h1, Helicity h2, Helicity hg>
	double amp_juno_h_j(
	HLV const & pa, HLV const & pb,
	HLV const & pg, HLV const & p1, HLV const & p2,
	HLV const & qH11, HLV const & qH12, HLV const & qH21, HLV const & qH22,
	std::array<COM, 2> const & T_pref
	) {
	// TODO: code duplication with Wjets and pure jets
	const COM u1 = U1_h_j<h1, h2, hg>(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]);
	const COM u2 = U2_h_j<h1, h2, hg>(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]);
	const COM l = L_h_j<h1, h2, hg>(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]);
	return 2.C_Fstd::real((l-u1)*std::conj(l+u2))
	+ 2.C_FC_F/3.*std::norm(u1+u2)
	;
	}


	//@{
	/**
	* @brief Higgs+Jets Unordered Contributions, function to handle all incoming types.
	* @param pg Unordered Gluon momenta
	* @param p1out Outgoing Particle 1
	* @param p1in Incoming Particle 1
	* @param p2out Outgoing Particle 2
	* @param p2in Incoming Particle 2
	* @param qH1 t-channel momenta into higgs vertex
	* @param qH2 t-channel momenta out of higgs vertex
	* @param mt top mass (inf or value)
	* @param inc_bottom whether to include bottom mass effects (true) or not (false)
	* @param mb bottom mass (value)
	*
	* Calculates j_{uno}^\mu H j_\mu. Unordered with higgs vertex
	* somewhere in the FKL chain. Handles all possible incoming states.
	*/
	double juno_h_j(
	HLV const & pg, HLV const & p1out, HLV const & p1in,
	HLV const & p2out, HLV const & p2in,
	HLV const & qH1, HLV const & qH2,
	const double mt, const bool incBot, const double mb, const double vev
	){
	using helicity::plus;
	using helicity::minus;

	const auto qqH1 = split_into_lightlike(qH1);
	const HLV qH11 = qqH1.first;
	const HLV qH12 = -qqH1.second;
	const auto qqH2 = split_into_lightlike(qH2);
	const HLV qH21 = qqH2.first;
	const HLV qH22 = -qqH2.second;
	// since qH1.m2(), qH2.m2() < 0 the following assertions are always true
	assert(qH11.e() > 0);
	assert(qH12.e() > 0);
	assert(qH21.e() > 0);
	assert(qH22.e() > 0);

	const auto T_pref = TT(qH1, qH2, mt, incBot, mb, vev);

	// only 4 out of the 8 helicity amplitudes are independent
	// we still compute all of them for better numerical stability (mirror check)
	MultiArray<double, 2, 2, 2> amp{};

	// NOLINTNEXTLINE
	#define ASSIGN_HEL(RES, J, H1, H2, HG) \
	RES[H1][H2][HG] = J<H1, H2, HG>( \
	p1in, p2in, pg, p1out, p2out, qH11, qH12, qH21, qH22, T_pref \
	)

	ASSIGN_HEL(amp, amp_juno_h_j, minus, minus, minus);
	ASSIGN_HEL(amp, amp_juno_h_j, minus, minus, plus);
	ASSIGN_HEL(amp, amp_juno_h_j, minus, plus, minus);
	ASSIGN_HEL(amp, amp_juno_h_j, minus, plus, plus);
	ASSIGN_HEL(amp, amp_juno_h_j, plus, minus, minus);
	ASSIGN_HEL(amp, amp_juno_h_j, plus, minus, plus);
	ASSIGN_HEL(amp, amp_juno_h_j, plus, plus, minus);
	ASSIGN_HEL(amp, amp_juno_h_j, plus, plus, plus);

	#undef ASSIGN_HEL

	const HLV q1 = p1in-p1out; // Top End
	const HLV q2 = p2out-p2in; // Bottom End
	const HLV qg = p1in-p1out-pg; // Extra bit post-gluon

	double ampsq = 0.;
	for(Helicity h1: {minus, plus}) {
	for(Helicity h2: {minus, plus}) {
	for(Helicity hg: {minus, plus}) {
	ampsq += amp[h1][h2][hg];
	}
	}
	}

	ampsq /= 16.qg.m2()qH1.m2()qH2.m2()q2.m2();

	// Factor of (Cf/Ca) for each quark to match ME_H_qQ.
	ampsq=C_FC_F/C_A/C_A;

	return ampsq;
	}
	} // namespace

	double ME_H_unob_qQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
	HLV const & p2out, HLV const & p2in, HLV const & qH1,
	HLV const & qH2, double mt, bool include_bottom, double mb,
	double vev
	){
	return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev);
	}

	double ME_H_unob_qbarQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
	HLV const & p2out, HLV const & p2in, HLV const & qH1,
	HLV const & qH2, double mt, bool include_bottom, double mb,
	double vev
	){
	return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev);
	}

	double ME_H_unob_qQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
	HLV const & p2out, HLV const & p2in, HLV const & qH1,
	HLV const & qH2, double mt, bool include_bottom, double mb,
	double vev
	){
	return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev);
	}

	double ME_H_unob_qbarQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
	HLV const & p2out, HLV const & p2in, HLV const & qH1,
	HLV const & qH2, double mt, bool include_bottom, double mb,
	double vev
	){
	return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev);
	}

	double ME_H_unob_gQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
	HLV const & p2out, HLV const & p2in, HLV const & qH1,
	HLV const & qH2, double mt, bool include_bottom, double mb,
	double vev
	){
	return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev)
	*K_g(p2out,p2in)/C_F;
	}

	double ME_H_unob_gQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
	HLV const & p2out, HLV const & p2in, HLV const & qH1,
	HLV const & qH2, double mt, bool include_bottom, double mb,
	double vev
	){
	return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev)
	*K_g(p2out,p2in)/C_F;
	}
	//@}
	} // namespace currents
	} // namespace HEJ
	diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
	index b8827c4..dc21bbf 100644
	--- a/src/MatrixElement.cc
	+++ b/src/MatrixElement.cc
	@@ -1,2439 +1,2444 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2022
	* \copyright GPLv2 or later
	*/
	#include "HEJ/MatrixElement.hh"

	#include <algorithm>
	#include <cassert>
	#include <cmath>
	#include <cstddef>
	#include <cstdlib>
	#include <iterator>
	#include <limits>
	#include <unordered_map>
	#include <utility>

	#include "fastjet/PseudoJet.hh"

	#include "HEJ/ConfigFlags.hh"
	#include "HEJ/Constants.hh"
	#include "HEJ/EWConstants.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/HiggsCouplingSettings.hh"
	#include "HEJ/Hjets.hh"
	#include "HEJ/LorentzVector.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/WWjets.hh"
	#include "HEJ/Wjets.hh"
	#include "HEJ/Zjets.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/jets.hh"
	#include "HEJ/utility.hh"

	namespace HEJ {
	namespace {

	// Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
	double K_g(double p1minus, double paminus) {
	return 1./2.(p1minus/paminus + paminus/p1minus)(C_A - 1./C_A) + 1./C_A;
	}
	double K_g(
	CLHEP::HepLorentzVector const & pout,
	CLHEP::HepLorentzVector const & pin
	) {
	if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
	return K_g(pout.minus(), pin.minus());
	}

	// Colour acceleration multiplier for quarks
	// see eq:K_q in developer manual
	constexpr double K_q = C_F;

	// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
	double K(
	ParticleID type,
	CLHEP::HepLorentzVector const & pout,
	CLHEP::HepLorentzVector const & pin
	){
	if(type == pid::gluon) return K_g(pout, pin);
	return K_q;
	}
	} // namespace

	double MatrixElement::omega0(
	double alpha_s, double mur,
	fastjet::PseudoJet const & q_j
	) const {
	const double lambda = param_.regulator_lambda;
	const double result = - alpha_sN_C/M_PIstd::log(q_j.perp2()/(lambda*lambda));
	if(! param_.log_correction) return result;
	return (
	1. + alpha_s/(4.M_PI)BETA0std::log(murmur/(q_j.perp()*lambda))
	)*result;
	}

	Weights MatrixElement::operator()(Event const & event) const {
	std::vector <double> tree_kin_part=tree_kin(event);
	std::vector <Weights> virtual_part=virtual_corrections(event);
	if(tree_kin_part.size() != virtual_part.size()) {
	throw std::logic_error("tree and virtuals have different sizes");
	}
	Weights sum = Weights{0., std::vector<double>(event.variations().size(), 0.)};
	for(size_t i=0; i<tree_kin_part.size(); ++i) {
	sum += tree_kin_part.at(i)*virtual_part.at(i);
	}
	return tree_param(event)*sum;
	}

	Weights MatrixElement::tree(Event const & event) const {
	std::vector <double> tree_kin_part=tree_kin(event);
	double sum = 0.;
	for(double i : tree_kin_part) {
	sum += i;
	}
	return tree_param(event)*sum;
	}

	Weights MatrixElement::tree_param(Event const & event) const {
	if(! is_resummable(event.type())) {
	return Weights{0., std::vector<double>(event.variations().size(), 0.)};
	}
	Weights result;
	// only compute once for each renormalisation scale
	std::unordered_map<double, double> known;
	result.central = tree_param(event, event.central().mur);
	known.emplace(event.central().mur, result.central);
	for(auto const & var: event.variations()) {
	const auto ME_it = known.find(var.mur);
	if(ME_it == end(known)) {
	const double wt = tree_param(event, var.mur);
	result.variations.emplace_back(wt);
	known.emplace(var.mur, wt);
	}
	else {
	result.variations.emplace_back(ME_it->second);
	}
	}
	return result;
	}

	std::vector<Weights> MatrixElement::virtual_corrections(Event const & event) const {
	if(!event.valid_hej_state(param_.soft_pt_regulator)) {
	return {Weights{0., std::vector<double>(event.variations().size(), 0.)}};
	}
	// only compute once for each renormalisation scale
	std::unordered_map<double, std::vector<double> > known_vec;
	std::vector<double> central_vec=virtual_corrections(event, event.central().mur);
	known_vec.emplace(event.central().mur, central_vec);
	for(auto const & var: event.variations()) {
	const auto ME_it = known_vec.find(var.mur);
	if(ME_it == end(known_vec)) {
	known_vec.emplace(var.mur, virtual_corrections(event, var.mur));
	}
	}
	// At this stage known_vec contains one vector of virtual corrections for each mur value
	// Now put this into a vector of Weights
	std::vector<Weights> result_vec;
	for(size_t i=0; i<central_vec.size(); ++i) {
	Weights result;
	result.central = central_vec.at(i);
	for(auto const & var: event.variations()) {
	const auto ME_it = known_vec.find(var.mur);
	result.variations.emplace_back(ME_it->second.at(i));
	}
	result_vec.emplace_back(result);
	}
	return result_vec;
	}

	template<class InputIterator>
	std::vector <double> MatrixElement::virtual_corrections_interference(
	InputIterator begin_parton, InputIterator end_parton,
	fastjet::PseudoJet const & q0_t,
	fastjet::PseudoJet const & q0_b,
	const double mur
	) const{
	const double alpha_s = alpha_s_(mur);

	auto qi_t = q0_t;
	auto qi_b = q0_b;

	double sum_top = 0.;
	double sum_bot = 0.;
	double sum_mix = 0.;
	for(auto parton_it = begin_parton; parton_it != end_parton; ++parton_it){
	Particle parton = *parton_it;
	Particle parton_next = *(parton_it+1);
	const double dy = parton_next.rapidity() - parton.rapidity();
	const double tmp_top = omega0(alpha_s, mur, qi_t)*dy;
	const double tmp_bot = omega0(alpha_s, mur, qi_b)*dy;
	sum_top += tmp_top;
	sum_bot += tmp_bot;
	sum_mix += (tmp_top + tmp_bot) / 2.;
	qi_t -= parton_next.p;
	qi_b -= parton_next.p;
	}

	if (param_.nlo.enabled){
	return {(sum_top), (sum_bot), (sum_mix)};
	}

	return {exp(sum_top), exp(sum_bot), exp(sum_mix)};


	}

	double MatrixElement::virtual_corrections_W(
	Event const & event,
	const double mur,
	Particle const & WBoson
	) const{
	auto const & in = event.incoming();
	const auto partons = filter_partons(event.outgoing());
	fastjet::PseudoJet const & pa = in.front().p;
	#ifndef NDEBUG
	fastjet::PseudoJet const & pb = in.back().p;
	double const norm = (in.front().p + in.back().p).E();
	#endif

	assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
	assert(partons.size() >= 2);
	assert(pa.pz() < pb.pz());

	fastjet::PseudoJet q = pa - partons[0].p;
	std::size_t first_idx = 0;
	std::size_t last_idx = partons.size() - 1;

	#ifndef NDEBUG
	bool wc = true;
	#endif
	bool wqq = false;

	// With extremal qqbar or unordered gluon outside the extremal
	// partons then it is not part of the FKL ladder and does not
	// contribute to the virtual corrections. W emitted from the
	// most backward leg must be taken into account in t-channel

	if (event.type() == event_type::unob) {
	q -= partons[1].p;
	++first_idx;
	if (in[0].type != partons[1].type ){
	q -= WBoson.p;
	#ifndef NDEBUG
	wc=false;
	#endif
	}
	}

	else if (event.type() == event_type::qqbar_exb) {
	q -= partons[1].p;
	++first_idx;
	if (std::abs(partons[0].type) != std::abs(partons[1].type)){
	q -= WBoson.p;
	#ifndef NDEBUG
	wc=false;
	#endif
	}
	}
	else {
	if(event.type() == event_type::unof
	\|\| event.type() == event_type::qqbar_exf){
	--last_idx;
	}
	if (in[0].type != partons[0].type ){
	q -= WBoson.p;
	#ifndef NDEBUG
	wc=false;
	#endif
	}
	}


	std::size_t first_idx_qqbar = last_idx;
	std::size_t last_idx_qqbar = last_idx;

	//if qqbarMid event, virtual correction do not occur between qqbar pair.
	if(event.type() == event_type::qqbar_mid){
	const auto backquark = std::find_if(
	begin(partons) + 1, end(partons) - 1 ,
	[](Particle const & s){ return (s.type != pid::gluon); }
	);
	if(backquark == end(partons) \|\| (backquark+1)->type==pid::gluon) return 0;
	if(std::abs(backquark->type) != std::abs((backquark+1)->type)) {
	wqq=true;
	#ifndef NDEBUG
	wc=false;
	#endif
	}
	last_idx = std::distance(begin(partons), backquark);
	first_idx_qqbar = last_idx+1;
	}
	double exponent = 0;
	const double alpha_s = alpha_s_(mur);
	for(std::size_t j = first_idx; j < last_idx; ++j){
	exponent += omega0(alpha_s, mur, q)*(
	partons[j+1].rapidity() - partons[j].rapidity()
	);
	q -=partons[j+1].p;
	} // End Loop one

	if (last_idx != first_idx_qqbar) q -= partons[last_idx+1].p;
	if (wqq) q -= WBoson.p;

	for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){
	exponent += omega0(alpha_s, mur, q)*(
	partons[j+1].rapidity() - partons[j].rapidity()
	);
	q -= partons[j+1].p;
	}

	#ifndef NDEBUG
	if (wc) q -= WBoson.p;
	assert(
	nearby(q, -1*pb, norm)
	\|\| is_AWZH_boson(partons.back().type)
	\|\| event.type() == event_type::unof
	\|\| event.type() == event_type::qqbar_exf
	);
	#endif
	if (param_.nlo.enabled){
	double nlo_virtual =1.;
	//Only apply virtual corrections to a nlo order event.
	if (partons.size() == param_.nlo.nj) nlo_virtual+=exponent;
	return nlo_virtual;
	}
	return std::exp(exponent);
	}



	std::vector <double> MatrixElement::virtual_corrections_WW(
	Event const & event,
	const double mur
	) const{
	auto const & in = event.incoming();
	const auto partons = filter_partons(event.outgoing());
	fastjet::PseudoJet const & pa = in.front().p;
	#ifndef NDEBUG
	fastjet::PseudoJet const & pb = in.back().p;
	#endif

	assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
	assert(partons.size() >= 2);
	assert(pa.pz() < pb.pz());
	assert(event.decays().size() == 2);

	std::vector<fastjet::PseudoJet> plbar;
	std::vector<fastjet::PseudoJet> pl;

	for(auto const & decay_pair : event.decays()) {
	auto const decay = decay_pair.second;
	if(decay.at(0).type < 0) {
	plbar.emplace_back(decay.at(0).p);
	pl .emplace_back(decay.at(1).p);
	}
	else {
	pl .emplace_back(decay.at(0).p);
	plbar.emplace_back(decay.at(1).p);
	}
	}

	fastjet::PseudoJet q_t = pa - partons[0].p - pl[0] - plbar[0];
	fastjet::PseudoJet q_b = pa - partons[0].p - pl[1] - plbar[1];

	auto const begin_parton = cbegin(partons);
	auto const end_parton = cend(partons) - 1;

	if (param_.nlo.enabled){
	std::vector<double> virt_corrections_nlo {1.0,1.0,1.0}; //set virtual corrections to 1.

	// Only apply virtual corrections to a nlo order event.
	if (partons.size() == param_.nlo.nj) {
	std::vector <double> virt_corrections_nlo_interference =
	virtual_corrections_interference(
	begin_parton, end_parton, q_t, q_b, mur
	);
	assert(
	virt_corrections_nlo_interference.size()
	== virt_corrections_nlo.size()
	);

	for(std::size_t i = 0; i < virt_corrections_nlo.size(); ++i){
	virt_corrections_nlo[i] += virt_corrections_nlo_interference[i];
	}
	}
	return virt_corrections_nlo;
	}

	return virtual_corrections_interference(begin_parton, end_parton, q_t, q_b, mur);

	}

	std::vector <double> MatrixElement::virtual_corrections_Z_qq(
	Event const & event,
	const double mur,
	Particle const & ZBoson
	) const{
	auto const & in = event.incoming();
	const auto partons = filter_partons(event.outgoing());
	fastjet::PseudoJet const & pa = in.front().p;
	#ifndef NDEBUG
	fastjet::PseudoJet const & pb = in.back().p;
	#endif

	assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
	assert(partons.size() >= 2);
	assert(pa.pz() < pb.pz());

	fastjet::PseudoJet q_t = pa - partons[0].p - ZBoson.p;
	fastjet::PseudoJet q_b = pa - partons[0].p;

	auto begin_parton = cbegin(partons);
	auto end_parton = cend(partons) - 1;

	// Unordered gluon does not contribute to the virtual corrections
	if (event.type() == event_type::unob) {
	// Gluon is partons[0] and is already subtracted
	// partons[1] is the backward quark
	q_t -= partons[1].p;
	q_b -= partons[1].p;
	++begin_parton;
	} else if (event.type() == event_type::unof) {
	// End sum at forward quark
	--end_parton;
	}

	if (param_.nlo.enabled){

	//set virtual corrections to 1.
	std::vector<double> virt_corrections_nlo {1.0,1.0,1.0};

	// Only apply virtual corrections to a nlo order event.
	if (partons.size() == param_.nlo.nj) {
	std::vector <double> virt_corrections_nlo_interference =
	virtual_corrections_interference(
	begin_parton, end_parton, q_t, q_b, mur
	);
	assert(
	virt_corrections_nlo.size()
	== virt_corrections_nlo_interference.size()
	);
	for(std::size_t i = 0; i < virt_corrections_nlo.size(); ++i){
	virt_corrections_nlo[i] += virt_corrections_nlo_interference[i];
	}
	}
	return virt_corrections_nlo;
	}

	return virtual_corrections_interference(begin_parton, end_parton, q_t, q_b, mur);
	}

	double MatrixElement::virtual_corrections_Z_qg(
	Event const & event,
	const double mur,
	Particle const & ZBoson,
	const bool is_gq_event
	) const{
	auto const & in = event.incoming();
	const auto partons = filter_partons(event.outgoing());
	fastjet::PseudoJet const & pa = in.front().p;
	#ifndef NDEBUG
	fastjet::PseudoJet const & pb = in.back().p;
	#endif

	assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
	assert(partons.size() >= 2);
	assert(pa.pz() < pb.pz());

	// If this is a gq event, don't subtract the Z momentum from first q
	fastjet::PseudoJet q = (is_gq_event ? pa - partons[0].p : pa - partons[0].p - ZBoson.p);
	size_t first_idx = 0;
	size_t last_idx = partons.size() - 1;

	// Unordered gluon does not contribute to the virtual corrections
	if (event.type() == event_type::unob) {
	// Gluon is partons[0] and is already subtracted
	// partons[1] is the backward quark
	q -= partons[1].p;
	++first_idx;
	} else if (event.type() == event_type::unof) {
	// End sum at forward quark
	--last_idx;
	}

	double sum=0.;
	const double alpha_s = alpha_s_(mur);
	for(size_t j = first_idx; j < last_idx; ++j){
	sum += omega0(alpha_s, mur, q)*(partons[j+1].rapidity()
	- partons[j].rapidity());
	q -= partons[j+1].p;
	}

	if (param_.nlo.enabled){
	double nlo_virtual =1.;
	//Only apply virtual corrections to a nlo order event.
	if (partons.size() == param_.nlo.nj) nlo_virtual+=sum;
	return nlo_virtual;
	}
	return exp(sum);
	}

	std::vector<double> MatrixElement::virtual_corrections(
	Event const & event,
	const double mur
	) const{
	auto const & in = event.incoming();
	auto const & out = event.outgoing();
	fastjet::PseudoJet const & pa = in.front().p;
	#ifndef NDEBUG
	fastjet::PseudoJet const & pb = in.back().p;
	double const norm = (in.front().p + in.back().p).E();
	#endif

	std::vector<Particle> bosons = filter_AWZH_bosons(out);

	if(event.jets().size() != param_.nlo.nj && param_.nlo.enabled) {
	throw std::logic_error{
	"Input event has number of jets different to stated NLO "
	"input in config file: " + std::to_string(event.jets().size())
	+ " vs " +std::to_string(param_.nlo.nj) + "\n"
	};
	}
	if(bosons.size() > 2) {
	throw not_implemented("Emission of >2 bosons is unsupported");
	}

	if(bosons.size() == 2) {
	if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp) {
	return virtual_corrections_WW(event, mur);
	}
	else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm) {
	return virtual_corrections_WW(event, mur);
	}
	throw not_implemented("Emission of bosons of unsupported type");
	}

	if(bosons.size() == 1) {
	const auto AWZH_boson = bosons[0];

	if(std::abs(AWZH_boson.type) == pid::Wp){
	return {virtual_corrections_W(event, mur, AWZH_boson)};
	}

	if(AWZH_boson.type == pid::Z_photon_mix){
	if(is_gluon(in.back().type)){
	// This is a qg event
	return {virtual_corrections_Z_qg(event, mur, AWZH_boson, false)};
	}
	if(is_gluon(in.front().type)){
	// This is a gq event
	return {virtual_corrections_Z_qg(event, mur, AWZH_boson, true)};
	}
	// This is a qq event
	return virtual_corrections_Z_qq(event, mur, AWZH_boson);
	}
	}


	assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
	assert(out.size() >= 2);
	assert(pa.pz() < pb.pz());

	fastjet::PseudoJet q = pa - out[0].p;
	std::size_t first_idx = 0;
	std::size_t last_idx = out.size() - 1;

	// if there is a Higgs boson _not_ emitted off an incoming gluon,
	// extremal qqbar or unordered gluon outside the extremal partons
	// then it is not part of the FKL ladder
	// and does not contribute to the virtual corrections
	if((out.front().type == pid::Higgs && in.front().type != pid::gluon)
	\|\| event.type() == event_type::unob
	\|\| event.type() == event_type::qqbar_exb){
	q -= out[1].p;
	++first_idx;
	}
	if((out.back().type == pid::Higgs && in.back().type != pid::gluon)
	\|\| event.type() == event_type::unof
	\|\| event.type() == event_type::qqbar_exf){
	--last_idx;
	}

	std::size_t first_idx_qqbar = last_idx;
	std::size_t last_idx_qqbar = last_idx;

	//if central qqbar event, virtual correction do not occur between q-qbar.
	if(event.type() == event_type::qqbar_mid){
	const auto backquark = std::find_if(
	begin(out) + 1, end(out) - 1 ,
	[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
	);
	if(backquark == end(out) \|\| (backquark+1)->type==pid::gluon) return {0.};
	last_idx = std::distance(begin(out), backquark);
	first_idx_qqbar = last_idx+1;
	}
	double exponent = 0;
	const double alpha_s = alpha_s_(mur);
	for(std::size_t j = first_idx; j < last_idx; ++j){
	exponent += omega0(alpha_s, mur, q)*(
	out[j+1].rapidity() - out[j].rapidity()
	);
	q -= out[j+1].p;
	}

	if (last_idx != first_idx_qqbar) q -= out[last_idx+1].p;

	for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){
	exponent += omega0(alpha_s, mur, q)*(
	out[j+1].rapidity() - out[j].rapidity()
	);
	q -= out[j+1].p;
	}
	assert(
	nearby(q, -1*pb, norm)
	\|\| (out.back().type == pid::Higgs && in.back().type != pid::gluon)
	\|\| event.type() == event_type::unof
	\|\| event.type() == event_type::qqbar_exf
	);

	if (param_.nlo.enabled){
	const auto partons = filter_partons(event.outgoing());
	double nlo_virtual =1.;
	//Only apply virtual corrections to a nlo order event.
	if (partons.size() == param_.nlo.nj) nlo_virtual+=exponent;
	return {nlo_virtual};
	}
	return {std::exp(exponent)};
	}

	namespace {

	//! Lipatov vertex for partons emitted into extremal jets
	CLHEP::HepLorentzVector CLipatov(
	CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
	CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
	) {
	const CLHEP::HepLorentzVector p5 = qav-qbv;
	const CLHEP::HepLorentzVector CL = -(qav+qbv)
	+ p1(qav.m2()/p5.dot(p1) + 2.p5.dot(p2)/p1.dot(p2))
	- p2(qbv.m2()/p5.dot(p2) + 2.p5.dot(p1)/p1.dot(p2));
	return CL;
	}

	double C2Lipatov(
	CLHEP::HepLorentzVector const & qav,
	CLHEP::HepLorentzVector const & qbv,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & p2
	){
	const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, p1, p2);
	return -CL.dot(CL);
	}

	//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
	double C2Lipatovots(
	CLHEP::HepLorentzVector const & qav,
	CLHEP::HepLorentzVector const & qbv,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & p2,
	const double lambda
	) {
	const double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));

	const double kperp=(qav-qbv).perp();
	if (kperp>lambda)
	return Cls;

	return Cls-4./(kperp*kperp);
	}

	double C2Lipatov_Mix(
	CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
	CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
	CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
	) {
	const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, p1, p2);
	const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, p1, p2);
	return -CL_t.dot(CL_b);
	}

	double C2Lipatovots_Mix(
	CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
	CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
	CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
	const double lambda
	) {
	const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, p1, p2)
	/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
	const double kperp = (qav_t - qbv_t).perp();
	if (kperp > lambda){
	return Cls;
	}
	return Cls - 4.0 / (kperp * kperp);

	}

	CLHEP::HepLorentzVector CLipatov(
	CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
	CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
	CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
	){
	const CLHEP::HepLorentzVector p5 = qav-qbv;
	const CLHEP::HepLorentzVector CL = -(qav+qbv)
	+ qav.m2()(1./p5.dot(pip)pip + 1./p5.dot(pop)*pop)/2.
	- qbv.m2()(1./p5.dot(pim)pim + 1./p5.dot(pom)*pom)/2.
	+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
	+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
	- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
	- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
	return CL;
	}

	//! Lipatov vertex
	double C2Lipatov( // B
	CLHEP::HepLorentzVector const & qav,
	CLHEP::HepLorentzVector const & qbv,
	CLHEP::HepLorentzVector const & pim,
	CLHEP::HepLorentzVector const & pip,
	CLHEP::HepLorentzVector const & pom,
	CLHEP::HepLorentzVector const & pop
	){
	const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, pim, pip, pom, pop);
	return -CL.dot(CL);
	}

	//! Lipatov vertex with soft subtraction
	double C2Lipatovots(
	CLHEP::HepLorentzVector const & qav,
	CLHEP::HepLorentzVector const & qbv,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & p2,
	const double lambda
	) {
	const double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));

	const double kperp=(qav-qbv).perp();
	if (kperp>lambda)
	return Cls;
	return Cls-4./(kperp*kperp);
	}

	double C2Lipatov_Mix(
	CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
	CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
	CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
	CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
	) {
	const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, pim, pip, pom, pop);
	const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, pim, pip, pom, pop);
	return -CL_t.dot(CL_b);
	}

	double C2Lipatovots_Mix(
	CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
	CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
	CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
	const double lambda
	) {
	const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, pa, pb, p1, p2)
	/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
	const double kperp = (qav_t - qbv_t).perp();
	if (kperp > lambda) {
	return Cls;
	}
	return Cls - 4.0 / (kperp * kperp);

	}

	/** Matrix element squared for tree-level current-current scattering
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pg Unordered gluon momentum
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @returns ME Squared for Tree-Level Current-Current Scattering
	*
	* @note The unof contribution can be calculated by reversing the argument ordering.
	*/
	double ME_uno_current(
	[[maybe_unused]] ParticleID aptype,
	ParticleID bptype,
	CLHEP::HepLorentzVector const & pg,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa
	){
	assert(aptype!=pid::gluon); // aptype cannot be gluon

	const double t1 = (pa - p1 - pg).m2();
	const double t2 = (pb - pn).m2();

	return K_q * K(bptype, pn, pb)currents::ME_unob_qq(pg, p1, pa, pn, pb) / (t1 t2);
	}

	/** Matrix element squared for tree-level current-current scattering
	* @param bptype Particle b PDG ID
	* @param pgin Incoming gluon momentum
	* @param p1 More backward (anti-)quark from splitting Momentum
	* @param p2 Less backward (anti-)quark from splitting Momentum
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @returns ME Squared for Tree-Level Current-Current Scattering
	*
	* @note The forward qqbar contribution can be calculated by reversing the
	* argument ordering.
	*/
	double ME_qqbar_current(
	ParticleID bptype,
	CLHEP::HepLorentzVector const & pgin,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & p2,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb
	){
	const double t1 = (pgin - p1 - p2).m2();
	const double t2 = (pn - pb).m2();

	return K(bptype, pn, pb)currents::ME_qqbar_qg(pb, pgin, pn, p2, p1) / (t1 t2);
	}

	/* \brief Matrix element squared for central qqbar tree-level current-current
	* scattering
	*
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param nabove Number of gluons emitted before central qqbarpair
	* @param nbelow Number of gluons emitted after central qqbarpair
	* @param pa Initial state a Momentum
	* @param pb Initial state b Momentum
	* @param pq Final state qbar Momentum
	* @param pqbar Final state q Momentum
	* @param partons Vector of all outgoing partons
	* @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering
	*/
	double ME_qqbar_mid_current(
	ParticleID aptype, ParticleID bptype, int nabove,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & pq,
	CLHEP::HepLorentzVector const & pqbar,
	std::vector<CLHEP::HepLorentzVector> const & partons
	){
	using namespace currents;
	// CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards)
	const bool swap_qqbar=pqbar.rapidity() < pq.rapidity();

	CLHEP::HepLorentzVector const & p1 = partons.front();
	CLHEP::HepLorentzVector const & pn = partons.back();

	const double t1 = (p1 - pa).m2();
	const double t2 = (pb - pn).m2();

	return K(aptype, p1, pa)
	*K(bptype, pn, pb)
	*ME_Cenqqbar_qq(
	pa, pb, partons,
	swap_qqbar, nabove
	) / (t1 * t2);
	}


	/** Matrix element squared for tree-level current-current scattering
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @returns ME Squared for Tree-Level Current-Current Scattering
	*/
	double ME_current(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa
	){
	const double t1 = (p1 - pa).m2();
	const double t2 = (pb - pn).m2();

	return K(aptype, p1, pa)
	* K(bptype, pn, pb)
	* currents::ME_qq(p1, pa, pn, pb)/(t1 * t2);
	}

	/** Matrix element squared for tree-level current-current scattering With W+Jets
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
	* @returns ME Squared for Tree-Level Current-Current Scattering
	*/
	double ME_W_current(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & plbar,
	CLHEP::HepLorentzVector const & pl,
	bool const wc, ParticleProperties const & Wprop
	){
	using namespace currents;
	// We know it cannot be gg incoming.
	assert(!(aptype==pid::gluon && bptype==pid::gluon));
	if (aptype==pid::gluon && bptype!=pid::gluon) {
	if (is_quark(bptype))
	return ME_W_qg(pn,plbar,pl,pb,p1,pa,Wprop);
	return ME_W_qbarg(pn,plbar,pl,pb,p1,pa,Wprop);
	}
	if (bptype==pid::gluon && aptype!=pid::gluon) {
	if (is_quark(aptype))
	return ME_W_qg(p1,plbar,pl,pa,pn,pb,Wprop);
	return ME_W_qbarg(p1,plbar,pl,pa,pn,pb,Wprop);
	}
	// they are both quark
	if (wc){ // emission off b, (first argument pbout)
	if (is_quark(bptype)) {
	if (is_quark(aptype))
	return ME_W_qQ(pn,plbar,pl,pb,p1,pa,Wprop);
	return ME_W_qQbar(pn,plbar,pl,pb,p1,pa,Wprop);
	}
	if (is_quark(aptype))
	return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa,Wprop);
	return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa,Wprop);
	}
	// emission off a, (first argument paout)
	if (is_quark(aptype)) {
	if (is_quark(bptype))
	return ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop);
	return ME_W_qQbar(p1,plbar,pl,pa,pn,pb,Wprop);
	}
	// a is anti-quark
	if (is_quark(bptype))
	return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb,Wprop);
	return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb,Wprop);
	}

	/** Matrix element squared for backwards uno tree-level current-current
	* scattering With W+Jets
	*
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param pg Unordered gluon momentum
	* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
	* @returns ME Squared for unob Tree-Level Current-Current Scattering
	*
	* @note The unof contribution can be calculated by reversing the argument ordering.
	*/
	double ME_W_uno_current(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pg,
	CLHEP::HepLorentzVector const & plbar,
	CLHEP::HepLorentzVector const & pl,
	bool const wc, ParticleProperties const & Wprop
	){
	using namespace currents;
	// we know they are not both gluons
	assert(bptype != pid::gluon \|\| aptype != pid::gluon);
	if (bptype == pid::gluon && aptype != pid::gluon) {
	// b gluon => W emission off a
	if (is_quark(aptype))
	return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	}
	// they are both quark
	if (wc) {// emission off b, i.e. b is first current
	if (is_quark(bptype)){
	if (is_quark(aptype))
	return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	}
	if (is_quark(aptype))
	return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	}
	// wc == false, emission off a, i.e. a is first current
	if (is_quark(aptype)) {
	if (is_quark(bptype)) //qq
	return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	//qqbar
	return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	}
	// a is anti-quark
	if (is_quark(bptype)) //qbarq
	return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	//qbarqbar
	return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	}

	/** \brief Matrix element squared for backward qqbar tree-level current-current
	* scattering With W+Jets
	*
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pa Initial state a Momentum
	* @param pb Initial state b Momentum
	* @param pq Final state q Momentum
	* @param pqbar Final state qbar Momentum
	* @param pn Final state n Momentum
	* @param plbar Final state anti-lepton momentum
	* @param pl Final state lepton momentum
	* @param swap_qqbar Ordering of qqbar pair. False: pqbar extremal.
	* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
	* @returns ME Squared for qqbarb Tree-Level Current-Current Scattering
	*
	* @note calculate forwards qqbar contribution by reversing argument ordering.
	*/
	double ME_W_qqbar_current(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & pq,
	CLHEP::HepLorentzVector const & pqbar,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & plbar,
	CLHEP::HepLorentzVector const & pl,
	bool const swap_qqbar, bool const wc,
	ParticleProperties const & Wprop
	){
	using namespace currents;
	// With qqbar we could have 2 incoming gluons and W Emission
	if (aptype==pid::gluon && bptype==pid::gluon) {
	//a gluon, b gluon gg->qqbarWg
	// This will be a wqqbar emission as there is no other possible W Emission
	// Site.
	if (swap_qqbar)
	return ME_WExqqbar_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb, Wprop);
	return ME_WExqqbar_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop);
	}
	assert(aptype==pid::gluon && bptype!=pid::gluon );
	//a gluon => W emission off b leg or qqbar
	if (!wc){ // W Emitted from backwards qqbar
	if (swap_qqbar)
	return ME_WExqqbar_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop);
	return ME_WExqqbar_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb, Wprop);
	}
	// W Must be emitted from forwards leg.
	return ME_W_Exqqbar_QQq(pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop);
	}

	/* \brief Matrix element squared for central qqbar tree-level current-current
	* scattering With W+Jets
	*
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param nabove Number of gluons emitted before central qqbarpair
	* @param nbelow Number of gluons emitted after central qqbarpair
	* @param pa Initial state a Momentum
	* @param pb Initial state b Momentum\
	* @param pq Final state qbar Momentum
	* @param pqbar Final state q Momentum
	* @param partons Vector of all outgoing partons
	* @param plbar Final state anti-lepton momentum
	* @param pl Final state lepton momentum
	* @param wqq Boolean. True siginfies W boson is emitted from Central qqbar
	* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
	* @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering
	*/
	double ME_W_qqbar_mid_current(
	ParticleID aptype, ParticleID bptype,
	int nabove, int nbelow,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & pq,
	CLHEP::HepLorentzVector const & pqbar,
	std::vector<CLHEP::HepLorentzVector> const & partons,
	CLHEP::HepLorentzVector const & plbar,
	CLHEP::HepLorentzVector const & pl,
	bool const wqq, bool const wc,
	ParticleProperties const & Wprop
	){
	using namespace currents;
	// CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards)
	const bool swap_qqbar=pqbar.rapidity() < pq.rapidity();
	double wt=1.;

	if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/C_F;
	if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F;

	if(wqq)
	return wt*ME_WCenqqbar_qq(pa, pb, pl, plbar, partons,
	is_antiquark(aptype),is_antiquark(bptype),
	swap_qqbar, nabove, Wprop);
	return wt*ME_W_Cenqqbar_qq(pa, pb, pl, plbar, partons,
	is_antiquark(aptype), is_antiquark(bptype),
	swap_qqbar, nabove, nbelow, wc, Wprop);
	}

	/** Matrix element squared for tree-level current-current scattering With Z+Jets
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param plbar Final state positron momentum
	* @param pl Final state electron momentum
	* @param Zprop Z properties
	* @param stw2 Value of sin(theta_w)^2
	* @param ctw Value of cos(theta_w)
	* @returns ME Squared for Tree-Level Current-Current Scattering
	*/
	std::vector<double> ME_Z_current(
	const ParticleID aptype, const ParticleID bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & plbar,
	CLHEP::HepLorentzVector const & pl,
	ParticleProperties const & Zprop,
	const double stw2, const double ctw
	){
	using namespace currents;
	// we know they are not both gluons
	assert(!is_gluon(aptype) \|\| !is_gluon(bptype));
	if(is_anyquark(aptype) && is_gluon(bptype)){
	// This is a qg event
	return { ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw) };
	}
	if(is_gluon(aptype) && is_anyquark(bptype)){
	// This is a gq event
	return { ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw) };
	}
	assert(is_anyquark(aptype) && is_anyquark(bptype));
	// This is a qq event
	return ME_Z_qQ(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
	}

	/** Matrix element squared for backwards uno tree-level current-current
	* scattering With Z+Jets
	*
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param pg Unordered gluon momentum
	* @param plbar Final state positron momentum
	* @param pl Final state electron momentum
	* @param Zprop Z properties
	* @param stw2 Value of sin(theta_w)^2
	* @param ctw Value of cos(theta_w)
	* @returns ME Squared for unob Tree-Level Current-Current Scattering
	*
	* @note The unof contribution can be calculated by reversing the argument ordering.
	*/

	std::vector<double> ME_Z_uno_current(
	const ParticleID aptype, const ParticleID bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pg,
	CLHEP::HepLorentzVector const & plbar,
	CLHEP::HepLorentzVector const & pl,
	ParticleProperties const & Zprop,
	const double stw2, const double ctw
	){
	using namespace currents;
	// we know they are not both gluons
	assert(!is_gluon(aptype) \|\| !is_gluon(bptype));
	if (is_anyquark(aptype) && is_gluon(bptype)) {
	// This is a qg event
	return { ME_Zuno_qg(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw) };
	}
	if (is_gluon(aptype) && is_anyquark(bptype)) {
	// This is a gq event
	return { ME_Zuno_qg(pb,pa,pg,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw) };
	}
	assert(is_anyquark(aptype) && is_anyquark(bptype));
	// This is a qq event
	return ME_Zuno_qQ(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
	}

	/** \brief Matrix element squared for tree-level current-current scattering with Higgs
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param qH t-channel momentum before Higgs
	* @param qHp1 t-channel momentum after Higgs
	* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
	*/
	double ME_Higgs_current(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
	CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
	double mt, bool include_bottom, double mb, double vev
	){
	+ const double t1 = (pa - p1).m2();
	+ const double tn = (pb - pn).m2();
	+
	return
	K(aptype, p1, pa)
	*K(bptype, pn, pb)
	- *currents::ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev);
	+ *currents::ME_H_qQ(
	+ pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev
	+ ) / (t1 * tn * qH.m2() * qHp1.m2());
	}

	/** \brief Current matrix element squared with Higgs and unordered backward emission
	* @param aptype Particle A PDG ID
	* @param bptype Particle B PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param pg Unordered back Particle Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param qH t-channel momentum before Higgs
	* @param qHp1 t-channel momentum after Higgs
	* @returns ME Squared with Higgs and unordered backward emission
	*
	* @note This function assumes unordered gluon backwards from pa-p1 current.
	* For unof, reverse call order
	*/
	double ME_Higgs_current_uno(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pg,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
	CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
	double mt, bool include_bottom, double mb, double vev
	){
	using namespace currents;
	if (bptype==pid::gluon && aptype!=pid::gluon) {
	if (is_quark(aptype))
	return ME_H_unob_gQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	return ME_H_unob_gQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	}
	// they are both quark
	if (is_quark(aptype)) {
	if (is_quark(bptype))
	return ME_H_unob_qQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	return ME_H_unob_qbarQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	}
	if (is_quark(bptype))
	return ME_H_unob_qQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	return ME_H_unob_qbarQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	}

	/** Matrix element squared for tree-level scattering with WW+Jets
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param pl1bar Particle l1bar Momentum
	* @param pl1 Particle l1 Momentum
	* @param pl2bar Particle l2bar Momentum
	* @param pl2 Particle l2 Momentum
	* @returns ME Squared for Tree-Level Current-Current Scattering
	*/
	std::vector <double> ME_WW_current(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pl1bar,
	CLHEP::HepLorentzVector const & pl1,
	CLHEP::HepLorentzVector const & pl2bar,
	CLHEP::HepLorentzVector const & pl2,
	ParticleProperties const & Wprop
	){
	using namespace currents;

	if (aptype > 0 && bptype > 0)
	return ME_WW_qQ(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
	if (aptype < 0 && bptype > 0)
	return ME_WW_qbarQ(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
	if (aptype > 0 && bptype < 0)
	return ME_WW_qQbar(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
	if (aptype < 0 && bptype < 0)
	return ME_WW_qbarQbar(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);

	throw std::logic_error("unreachable");
	}

	void validate(MatrixElementConfig const & config) {
	#ifndef HEJ_BUILD_WITH_QCDLOOP
	if(!config.Higgs_coupling.use_impact_factors) {
	throw std::invalid_argument{
	"Invalid Higgs coupling settings.\n"
	"HEJ without QCDloop support can only use impact factors.\n"
	"Set use_impact_factors to true or recompile HEJ.\n"
	};
	}
	#endif
	if(config.Higgs_coupling.use_impact_factors
	&& config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) {
	throw std::invalid_argument{
	"Conflicting settings: "
	"impact factors may only be used in the infinite top mass limit"
	};
	}
	}
	} // namespace

	MatrixElement::MatrixElement(
	std::function<double (double)> alpha_s,
	MatrixElementConfig conf
	):
	alpha_s_{std::move(alpha_s)},
	param_{std::move(conf)}
	{
	validate(param_);
	}

	std::vector<double> MatrixElement::tree_kin(
	Event const & ev
	) const {
	if(!ev.valid_hej_state(param_.soft_pt_regulator)) return {0.};

	if(!ev.valid_incoming()){
	throw std::invalid_argument{
	"Invalid momentum for one or more incoming particles. "
	"Incoming momenta must have vanishing mass and transverse momentum."
	};
	}

	std::vector<Particle> bosons = filter_AWZH_bosons(ev.outgoing());

	if(bosons.empty()) {
	return {tree_kin_jets(ev)};
	}

	if(bosons.size() == 1) {
	switch(bosons[0].type){
	case pid::Higgs:
	return {tree_kin_Higgs(ev)};
	case pid::Wp:
	case pid::Wm:
	return {tree_kin_W(ev)};
	case pid::Z_photon_mix:
	return tree_kin_Z(ev);
	// TODO
	case pid::photon:
	case pid::Z:
	default:
	throw not_implemented("Emission of boson of unsupported type");
	}
	}

	if(bosons.size() == 2) {
	if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp){
	return tree_kin_WW(ev);
	}
	else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm){
	return tree_kin_WW(ev);
	}

	throw not_implemented("Emission of bosons of unsupported type");
	}

	throw not_implemented("Emission of >2 bosons is unsupported");
	}

	namespace {
	constexpr int EXTREMAL_JET_IDX = 1;
	constexpr int NO_EXTREMAL_JET_IDX = 0;

	bool treat_as_extremal(Particle const & parton){
	return parton.p.user_index() == EXTREMAL_JET_IDX;
	}

	template<class InputIterator>
	double FKL_ladder_weight(
	InputIterator begin_gluon, InputIterator end_gluon,
	CLHEP::HepLorentzVector const & q0,
	CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
	double lambda
	){
	double wt = 1;
	auto qi = q0;
	for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
	assert(gluon_it->type == pid::gluon);
	const auto g = to_HepLorentzVector(*gluon_it);
	const auto qip1 = qi - g;

	if(treat_as_extremal(*gluon_it)){
	wt = C2Lipatovots(qip1, qi, pa, pb, lambda)C_A;
	} else{
	wt = C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)C_A;
	}

	qi = qip1;
	}
	return wt;
	}

	template<class InputIterator>
	std::vector <double> FKL_ladder_weight_mix(
	InputIterator begin_gluon, InputIterator end_gluon,
	CLHEP::HepLorentzVector const & q0_t, CLHEP::HepLorentzVector const & q0_b,
	CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
	const double lambda
	){
	double wt_top = 1;
	double wt_bot = 1;
	double wt_mix = 1;
	auto qi_t = q0_t;
	auto qi_b = q0_b;
	for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
	assert(gluon_it->type == pid::gluon);
	const auto g = to_HepLorentzVector(*gluon_it);
	const auto qip1_t = qi_t - g;
	const auto qip1_b = qi_b - g;
	if(treat_as_extremal(*gluon_it)){
	wt_top = C2Lipatovots(qip1_t, qi_t, pa, pb, lambda)C_A;
	wt_bot = C2Lipatovots(qip1_b, qi_b, pa, pb, lambda)C_A;
	wt_mix = C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, lambda)C_A;
	} else{
	wt_top = C2Lipatovots(qip1_t, qi_t, pa, pb, p1, pn, lambda)C_A;
	wt_bot = C2Lipatovots(qip1_b, qi_b, pa, pb, p1, pn, lambda)C_A;
	wt_mix = C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, p1, pn, lambda)C_A;
	}
	qi_t = qip1_t;
	qi_b = qip1_b;
	}
	return {wt_top, wt_bot, wt_mix};
	}

	std::vector<Particle> tag_extremal_jet_partons( Event const & ev ){
	auto out_partons = filter_partons(ev.outgoing());
	if(out_partons.size() == ev.jets().size()){
	// no additional emissions in extremal jets, don't need to tag anything
	for(auto & parton: out_partons){
	parton.p.set_user_index(NO_EXTREMAL_JET_IDX);
	}
	return out_partons;
	}
	auto const & jets = ev.jets();
	std::vector<fastjet::PseudoJet> extremal_jets;
	if(! is_backward_g_to_h(ev)) {
	auto most_backward = begin(jets);
	// skip jets caused by unordered emission or qqbar
	if(ev.type() == event_type::unob \|\| ev.type() == event_type::qqbar_exb){
	assert(jets.size() >= 2);
	++most_backward;
	}
	extremal_jets.emplace_back(*most_backward);
	}
	if(! is_forward_g_to_h(ev)) {
	auto most_forward = end(jets) - 1;
	if(ev.type() == event_type::unof \|\| ev.type() == event_type::qqbar_exf){
	assert(jets.size() >= 2);
	--most_forward;
	}
	extremal_jets.emplace_back(*most_forward);
	}
	const auto extremal_jet_indices = ev.particle_jet_indices(
	extremal_jets
	);
	assert(extremal_jet_indices.size() == out_partons.size());
	for(std::size_t i = 0; i < out_partons.size(); ++i){
	assert(is_parton(out_partons[i]));
	const int idx = (extremal_jet_indices[i]>=0)?
	EXTREMAL_JET_IDX:
	NO_EXTREMAL_JET_IDX;
	out_partons[i].p.set_user_index(idx);
	}
	return out_partons;
	}

	double tree_kin_jets_qqbar_mid(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
	std::vector<Particle> const & partons,
	double lambda
	){
	CLHEP::HepLorentzVector pq;
	CLHEP::HepLorentzVector pqbar;
	const auto backmidquark = std::find_if(
	begin(partons)+1, end(partons)-1,
	[](Particle const & s){ return s.type != pid::gluon; }
	);

	assert(backmidquark!=end(partons)-1);

	if (is_quark(backmidquark->type)){
	pq = to_HepLorentzVector(*backmidquark);
	pqbar = to_HepLorentzVector(*(backmidquark+1));
	}
	else {
	pqbar = to_HepLorentzVector(*backmidquark);
	pq = to_HepLorentzVector(*(backmidquark+1));
	}

	auto p1 = to_HepLorentzVector(partons[0]);
	auto pn = to_HepLorentzVector(partons[partons.size() - 1]);

	auto q0 = pa - p1;
	// t-channel momentum after qqbar
	auto qqbart = q0;

	const auto begin_ladder = cbegin(partons) + 1;
	const auto end_ladder_1 = (backmidquark);
	const auto begin_ladder_2 = (backmidquark+2);
	const auto end_ladder = cend(partons) - 1;
	for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
	qqbart -= to_HepLorentzVector(*parton_it);
	}

	const int nabove = std::distance(begin_ladder, backmidquark);

	std::vector<CLHEP::HepLorentzVector> partonsHLV;
	partonsHLV.reserve(partons.size());
	for (std::size_t i = 0; i != partons.size(); ++i) {
	partonsHLV.push_back(to_HepLorentzVector(partons[i]));
	}

	const double current_factor = ME_qqbar_mid_current(
	aptype, bptype, nabove, pa, pb,
	pq, pqbar, partonsHLV
	);

	const double ladder_factor = FKL_ladder_weight(
	begin_ladder, end_ladder_1,
	q0, pa, pb, p1, pn,
	lambda
	)*FKL_ladder_weight(
	begin_ladder_2, end_ladder,
	qqbart, pa, pb, p1, pn,
	lambda
	);
	return current_factor*ladder_factor;
	}


	template<class InIter, class partIter>
	double tree_kin_jets_qqbar(InIter BeginIn, InIter EndIn, partIter BeginPart,
	partIter EndPart, double lambda){
	const auto pgin = to_HepLorentzVector(*BeginIn);
	const auto pb = to_HepLorentzVector(*(EndIn-1));
	const auto p1 = to_HepLorentzVector(*BeginPart);
	const auto p2 = to_HepLorentzVector(*(BeginPart+1));
	const auto pn = to_HepLorentzVector(*(EndPart-1));

	assert((BeginIn)->type==pid::gluon); // Incoming a must be gluon.
	const double current_factor = ME_qqbar_current(
	(EndIn-1)->type, pgin, p1, p2, pn, pb
	)/(4.(N_CN_C - 1.));
	const double ladder_factor = FKL_ladder_weight(
	(BeginPart+2), (EndPart-1),
	pgin-p1-p2, pgin, pb, p1, pn, lambda
	);

	return current_factor*ladder_factor;
	}

	template<class InIter, class partIter>
	double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
	partIter EndPart, double lambda
	){

	const auto pa = to_HepLorentzVector(*BeginIn);
	const auto pb = to_HepLorentzVector(*(EndIn-1));

	const auto pg = to_HepLorentzVector(*BeginPart);
	const auto p1 = to_HepLorentzVector(*(BeginPart+1));
	const auto pn = to_HepLorentzVector(*(EndPart-1));

	const double current_factor = ME_uno_current(
	(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa
	)/(4.(N_CN_C - 1.));
	const double ladder_factor = FKL_ladder_weight(
	(BeginPart+2), (EndPart-1),
	pa-p1-pg, pa, pb, p1, pn, lambda
	);

	return current_factor*ladder_factor;
	}
	} // namespace

	double MatrixElement::tree_kin_jets(Event const & ev) const {
	auto const & incoming = ev.incoming();
	const auto partons = tag_extremal_jet_partons(ev);

	if (ev.type()==event_type::FKL){
	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	const auto p1 = to_HepLorentzVector(partons.front());
	const auto pn = to_HepLorentzVector(partons.back());
	return ME_current(
	incoming[0].type, incoming[1].type,
	pn, pb, p1, pa
	)/(4.(N_CN_C - 1.))*FKL_ladder_weight(
	begin(partons) + 1, end(partons) - 1,
	pa - p1, pa, pb, p1, pn,
	param_.regulator_lambda
	);
	}
	if (ev.type()==event_type::unordered_backward){
	return tree_kin_jets_uno(incoming.begin(), incoming.end(),
	partons.begin(), partons.end(),
	param_.regulator_lambda);
	}
	if (ev.type()==event_type::unordered_forward){
	return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
	partons.rbegin(), partons.rend(),
	param_.regulator_lambda);
	}
	if (ev.type()==event_type::extremal_qqbar_backward){
	return tree_kin_jets_qqbar(incoming.begin(), incoming.end(),
	partons.begin(), partons.end(),
	param_.regulator_lambda);
	}
	if (ev.type()==event_type::extremal_qqbar_forward){
	return tree_kin_jets_qqbar(incoming.rbegin(), incoming.rend(),
	partons.rbegin(), partons.rend(),
	param_.regulator_lambda);
	}
	if (ev.type()==event_type::central_qqbar){
	return tree_kin_jets_qqbar_mid(incoming[0].type, incoming[1].type,
	to_HepLorentzVector(incoming[0]),
	to_HepLorentzVector(incoming[1]),
	partons, param_.regulator_lambda);
	}
	throw std::logic_error("Cannot reweight non-resummable processes in Pure Jets");
	}

	namespace {
	double tree_kin_W_FKL(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
	std::vector<Particle> const & partons,
	CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
	double lambda, ParticleProperties const & Wprop
	){
	auto p1 = to_HepLorentzVector(partons[0]);
	auto pn = to_HepLorentzVector(partons[partons.size() - 1]);

	const auto begin_ladder = cbegin(partons) + 1;
	const auto end_ladder = cend(partons) - 1;

	bool wc = aptype==partons[0].type; //leg b emits w
	auto q0 = pa - p1;
	if(!wc)
	q0 -= pl + plbar;

	const double current_factor = ME_W_current(
	aptype, bptype, pn, pb,
	p1, pa, plbar, pl, wc, Wprop
	);

	const double ladder_factor = FKL_ladder_weight(
	begin_ladder, end_ladder,
	q0, pa, pb, p1, pn,
	lambda
	);
	return current_factor*ladder_factor;
	}


	template<class InIter, class partIter>
	double tree_kin_W_uno(InIter BeginIn, partIter BeginPart,
	partIter EndPart,
	const CLHEP::HepLorentzVector & plbar,
	const CLHEP::HepLorentzVector & pl,
	double lambda, ParticleProperties const & Wprop
	){
	const auto pa = to_HepLorentzVector(*BeginIn);
	const auto pb = to_HepLorentzVector(*(BeginIn+1));

	const auto pg = to_HepLorentzVector(*BeginPart);
	const auto p1 = to_HepLorentzVector(*(BeginPart+1));
	const auto pn = to_HepLorentzVector(*(EndPart-1));

	bool wc = (BeginIn)->type==(BeginPart+1)->type; //leg b emits w
	auto q0 = pa - p1 - pg;
	if(!wc)
	q0 -= pl + plbar;

	const double current_factor = ME_W_uno_current(
	(BeginIn)->type, (BeginIn+1)->type, pn, pb,
	p1, pa, pg, plbar, pl, wc, Wprop
	);

	const double ladder_factor = FKL_ladder_weight(
	BeginPart+2, EndPart-1,
	q0, pa, pb, p1, pn,
	lambda
	);
	return current_factor*ladder_factor;
	}

	template<class InIter, class partIter>
	double tree_kin_W_qqbar(InIter BeginIn, partIter BeginPart,
	partIter EndPart,
	const CLHEP::HepLorentzVector & plbar,
	const CLHEP::HepLorentzVector & pl,
	double lambda, ParticleProperties const & Wprop
	){
	const bool swap_qqbar=is_quark(*BeginPart);
	const auto pa = to_HepLorentzVector(*BeginIn);
	const auto pb = to_HepLorentzVector(*(BeginIn+1));
	const auto pq = to_HepLorentzVector(*(BeginPart+(swap_qqbar?0:1)));
	const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_qqbar?1:0)));
	const auto p1 = to_HepLorentzVector(*(BeginPart));
	const auto pn = to_HepLorentzVector(*(EndPart-1));

	const bool wc = (BeginIn+1)->type!=(EndPart-1)->type; //leg b emits w
	auto q0 = pa - pq - pqbar;
	if(!wc)
	q0 -= pl + plbar;

	const double current_factor = ME_W_qqbar_current(
	(BeginIn)->type, (BeginIn+1)->type, pa, pb,
	pq, pqbar, pn, plbar, pl, swap_qqbar, wc, Wprop
	);

	const double ladder_factor = FKL_ladder_weight(
	BeginPart+2, EndPart-1,
	q0, pa, pb, p1, pn,
	lambda
	);
	return current_factor*ladder_factor;
	}

	double tree_kin_W_qqbar_mid(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pb,
	std::vector<Particle> const & partons,
	CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
	double lambda, ParticleProperties const & Wprop
	){
	CLHEP::HepLorentzVector pq;
	CLHEP::HepLorentzVector pqbar;
	const auto backmidquark = std::find_if(
	begin(partons)+1, end(partons)-1,
	[](Particle const & s){ return s.type != pid::gluon; }
	);

	assert(backmidquark!=end(partons)-1);

	if (is_quark(backmidquark->type)){
	pq = to_HepLorentzVector(*backmidquark);
	pqbar = to_HepLorentzVector(*(backmidquark+1));
	}
	else {
	pqbar = to_HepLorentzVector(*backmidquark);
	pq = to_HepLorentzVector(*(backmidquark+1));
	}

	auto p1 = to_HepLorentzVector(partons.front());
	auto pn = to_HepLorentzVector(partons.back());

	auto q0 = pa - p1;
	// t-channel momentum after qqbar
	auto qqbart = q0;

	bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqbar emit W
	bool wc = !wqq && (aptype==partons.front().type); //leg b emits w
	assert(!wqq \|\| !wc);
	if(wqq){ // emission from qqbar
	qqbart -= pl + plbar;
	} else if(!wc) { // emission from leg a
	q0 -= pl + plbar;
	qqbart -= pl + plbar;
	}

	const auto begin_ladder = cbegin(partons) + 1;
	const auto end_ladder_1 = (backmidquark);
	const auto begin_ladder_2 = (backmidquark+2);
	const auto end_ladder = cend(partons) - 1;
	for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
	qqbart -= to_HepLorentzVector(*parton_it);
	}

	const int nabove = std::distance(begin_ladder, backmidquark);
	const int nbelow = std::distance(begin_ladder_2, end_ladder);

	std::vector<CLHEP::HepLorentzVector> partonsHLV;
	partonsHLV.reserve(partons.size());
	for (std::size_t i = 0; i != partons.size(); ++i) {
	partonsHLV.push_back(to_HepLorentzVector(partons[i]));
	}

	const double current_factor = ME_W_qqbar_mid_current(
	aptype, bptype, nabove, nbelow, pa, pb,
	pq, pqbar, partonsHLV, plbar, pl, wqq, wc, Wprop
	);

	const double ladder_factor = FKL_ladder_weight(
	begin_ladder, end_ladder_1,
	q0, pa, pb, p1, pn,
	lambda
	)*FKL_ladder_weight(
	begin_ladder_2, end_ladder,
	qqbart, pa, pb, p1, pn,
	lambda
	);
	return current_factor*ladder_factor;
	}
	} // namespace

	double MatrixElement::tree_kin_W(Event const & ev) const {
	using namespace event_type;
	auto const & incoming(ev.incoming());

	#ifndef NDEBUG
	// assert that there is exactly one decay corresponding to the W
	assert(ev.decays().size() == 1);
	auto const & w_boson{
	std::find_if(ev.outgoing().cbegin(), ev.outgoing().cend(),
	[] (Particle const & p) -> bool {
	return std::abs(p.type) == ParticleID::Wp;
	}) };
	assert(w_boson != ev.outgoing().cend());
	assert( static_cast<long int>(ev.decays().cbegin()->first)
	== std::distance(ev.outgoing().cbegin(), w_boson) );
	#endif

	// find decay products of W
	auto const & decay{ ev.decays().cbegin()->second };
	assert(decay.size() == 2);
	assert( ( is_anylepton(decay.at(0)) && is_anyneutrino(decay.at(1)) )
	\|\| ( is_anylepton(decay.at(1)) && is_anyneutrino(decay.at(0)) ) );

	// get lepton & neutrino
	CLHEP::HepLorentzVector plbar;
	CLHEP::HepLorentzVector pl;
	if (decay.at(0).type < 0){
	plbar = to_HepLorentzVector(decay.at(0));
	pl = to_HepLorentzVector(decay.at(1));
	}
	else{
	pl = to_HepLorentzVector(decay.at(0));
	plbar = to_HepLorentzVector(decay.at(1));
	}

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	const auto partons = tag_extremal_jet_partons(ev);

	if(ev.type() == FKL){
	return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
	pa, pb, partons, plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Wprop());
	}
	if(ev.type() == unordered_backward){
	return tree_kin_W_uno(cbegin(incoming), cbegin(partons),
	cend(partons), plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Wprop());
	}
	if(ev.type() == unordered_forward){
	return tree_kin_W_uno(crbegin(incoming), crbegin(partons),
	crend(partons), plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Wprop());
	}
	if(ev.type() == extremal_qqbar_backward){
	return tree_kin_W_qqbar(cbegin(incoming), cbegin(partons),
	cend(partons), plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Wprop());
	}
	if(ev.type() == extremal_qqbar_forward){
	return tree_kin_W_qqbar(crbegin(incoming), crbegin(partons),
	crend(partons), plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Wprop());
	}
	assert(ev.type() == central_qqbar);
	return tree_kin_W_qqbar_mid(incoming[0].type, incoming[1].type,
	pa, pb, partons, plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Wprop());
	}

	namespace /* WW */ {
	std::vector <double> tree_kin_WW_FKL(
	ParticleID aptype, ParticleID bptype,
	CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
	std::vector<Particle> const & partons,
	CLHEP::HepLorentzVector const & pl1bar, CLHEP::HepLorentzVector const & pl1,
	CLHEP::HepLorentzVector const & pl2bar, CLHEP::HepLorentzVector const & pl2,
	double lambda, ParticleProperties const & Wprop
	){
	auto p1 = to_HepLorentzVector(partons[0]);
	auto pn = to_HepLorentzVector(partons[partons.size() - 1]);

	const std::vector <double> current_factor = ME_WW_current(
	aptype, bptype, pn, pb, p1, pa,
	pl1bar, pl1, pl2bar, pl2,
	Wprop
	);

	auto const begin_ladder = cbegin(partons) + 1;
	auto const end_ladder = cend(partons) - 1;

	// pa -> W1 p1, pb -> W2 + pn
	const auto q0 = pa - p1 - (pl1 + pl1bar);
	// pa -> W2 p1, pb -> W1 + pn
	const auto q1 = pa - p1 - (pl2 + pl2bar);

	const std::vector <double> ladder_factor = FKL_ladder_weight_mix(
	begin_ladder, end_ladder,
	q0, q1, pa, pb, p1, pn,
	lambda
	);

	std::vector <double> result;
	for(size_t i=0; i<current_factor.size(); ++i){
	result.push_back(current_factor.at(i)*ladder_factor.at(i));
	}
	return result;
	}
	} // namespace

	std::vector <double> MatrixElement::tree_kin_WW(Event const & ev) const {
	using namespace event_type;
	auto const & incoming(ev.incoming());

	auto const pa = to_HepLorentzVector(incoming[0]);
	auto const pb = to_HepLorentzVector(incoming[1]);

	auto const partons = tag_extremal_jet_partons(ev);

	// W1 & W2
	assert(ev.decays().size() == 2);

	std::vector<CLHEP::HepLorentzVector> plbar;
	std::vector<CLHEP::HepLorentzVector> pl;

	for(auto const & decay_pair : ev.decays()) {
	auto const decay = decay_pair.second;
	// TODO: how to label W1, W2
	if(decay.at(0).type < 0) {
	plbar.emplace_back(to_HepLorentzVector(decay.at(0)));
	pl .emplace_back(to_HepLorentzVector(decay.at(1)));
	}
	else {
	pl .emplace_back(to_HepLorentzVector(decay.at(0)));
	plbar.emplace_back(to_HepLorentzVector(decay.at(1)));
	}
	}

	if(ev.type() == FKL) {

	return tree_kin_WW_FKL(
	incoming[0].type, incoming[1].type,
	pa, pb, partons,
	plbar[0], pl[0], plbar[1], pl[1],
	param_.regulator_lambda,
	param_.ew_parameters.Wprop()
	);

	}
	throw std::logic_error("Can only reweight FKL events in WW");
	}

	namespace{
	std::vector <double> tree_kin_Z_FKL(
	const ParticleID aptype, const ParticleID bptype,
	CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
	std::vector<Particle> const & partons,
	CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
	const double lambda, ParticleProperties const & Zprop,
	const double stw2, const double ctw
	){
	const auto p1 = to_HepLorentzVector(partons[0]);
	const auto pn = to_HepLorentzVector(partons[partons.size() - 1]);

	const auto begin_ladder = cbegin(partons) + 1;
	const auto end_ladder = cend(partons) - 1;

	const std::vector <double> current_factor = ME_Z_current(
	aptype, bptype, pn, pb, p1, pa,
	plbar, pl, Zprop, stw2, ctw
	);

	std::vector <double> ladder_factor;
	if(is_gluon(bptype)){
	// This is a qg event
	const auto q0 = pa-p1-plbar-pl;
	ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
	q0, pa, pb, p1, pn, lambda));
	} else if(is_gluon(aptype)){
	// This is a gq event
	const auto q0 = pa-p1;
	ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
	q0, pa, pb, p1, pn, lambda));
	} else {
	// This is a qq event
	const auto q0 = pa-p1-plbar-pl;
	const auto q1 = pa-p1;
	ladder_factor=FKL_ladder_weight_mix(begin_ladder, end_ladder,
	q0, q1, pa, pb, p1, pn, lambda);
	}

	std::vector <double> result;
	for(size_t i=0; i<current_factor.size(); ++i){
	result.push_back(current_factor.at(i)*ladder_factor.at(i));
	}
	return result;
	}

	template<class InIter, class partIter>
	std::vector <double> tree_kin_Z_uno(InIter BeginIn, partIter BeginPart, partIter EndPart,
	const CLHEP::HepLorentzVector & plbar,
	const CLHEP::HepLorentzVector & pl,
	const double lambda, ParticleProperties const & Zprop,
	const double stw2, const double ctw){
	const auto pa = to_HepLorentzVector(*BeginIn);
	const auto pb = to_HepLorentzVector(*(BeginIn+1));

	const auto pg = to_HepLorentzVector(*BeginPart);
	const auto p1 = to_HepLorentzVector(*(BeginPart+1));
	const auto pn = to_HepLorentzVector(*(EndPart-1));

	const ParticleID aptype = (BeginIn)->type;
	const ParticleID bptype = (BeginIn+1)->type;

	const std::vector <double> current_factor = ME_Z_uno_current(
	aptype, bptype, pn, pb, p1, pa, pg,
	plbar, pl, Zprop, stw2, ctw
	);

	std::vector <double> ladder_factor;
	if(is_gluon(bptype)){
	// This is a qg event
	const auto q0 = pa-pg-p1-plbar-pl;
	ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
	q0, pa, pb, p1, pn, lambda));
	}else if(is_gluon(aptype)){
	// This is a gq event
	const auto q0 = pa-pg-p1;
	ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
	q0, pa, pb, p1, pn, lambda));
	}else{
	// This is a qq event
	const auto q0 = pa-pg-p1-plbar-pl;
	const auto q1 = pa-pg-p1;
	ladder_factor=FKL_ladder_weight_mix(BeginPart+2, EndPart-1,
	q0, q1, pa, pb, p1, pn, lambda);
	}

	std::vector <double> result;
	for(size_t i=0; i<current_factor.size(); ++i){
	result.push_back(current_factor.at(i)*ladder_factor.at(i));
	}
	return result;
	}

	} // namespace

	std::vector<double> MatrixElement::tree_kin_Z(Event const & ev) const {
	using namespace event_type;
	auto const & incoming(ev.incoming());

	// find decay products of Z
	auto const & decay{ ev.decays().cbegin()->second };
	assert(decay.size() == 2);
	assert(is_anylepton(decay.at(0)) && !is_anyneutrino(decay.at(0))
	&& decay.at(0).type==-decay.at(1).type);

	// get leptons
	CLHEP::HepLorentzVector plbar;
	CLHEP::HepLorentzVector pl;
	if (decay.at(0).type < 0){
	plbar = to_HepLorentzVector(decay.at(0));
	pl = to_HepLorentzVector(decay.at(1));
	}
	else{
	pl = to_HepLorentzVector(decay.at(0));
	plbar = to_HepLorentzVector(decay.at(1));
	}

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	const auto partons = tag_extremal_jet_partons(ev);

	const double stw2 = param_.ew_parameters.sin2_tw();
	const double ctw = param_.ew_parameters.cos_tw();

	if(ev.type() == FKL){
	return tree_kin_Z_FKL(incoming[0].type, incoming[1].type,
	pa, pb, partons, plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Zprop(),
	stw2, ctw);
	}
	if(ev.type() == unordered_backward){
	return tree_kin_Z_uno(cbegin(incoming), cbegin(partons),
	cend(partons), plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Zprop(),
	stw2, ctw);
	}
	if(ev.type() == unordered_forward){
	return tree_kin_Z_uno(crbegin(incoming), crbegin(partons),
	crend(partons), plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Zprop(),
	stw2, ctw);
	}
	throw std::logic_error("Can only reweight FKL or uno processes in Z+Jets");
	}

	double MatrixElement::tree_kin_Higgs(Event const & ev) const {
	if(ev.outgoing().front().type == pid::Higgs){
	return tree_kin_Higgs_first(ev);
	}
	if(ev.outgoing().back().type == pid::Higgs){
	return tree_kin_Higgs_last(ev);
	}
	return tree_kin_Higgs_between(ev);
	}

	// kinetic matrix element square for backward Higgs emission
	// cf. eq:ME_h_jets_peripheral in developer manual,
	// but without factors \alpha_s and the FKL ladder
	double MatrixElement::MH2_backwardH(
	const ParticleID type_forward,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & pH,
	CLHEP::HepLorentzVector const & pn
	) const {
	using namespace currents;
	const double vev = param_.ew_parameters.vev();
	return K(type_forward, pn, pb)*ME_jgH_j(
	pH, pa, pn, pb,
	param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
	param_.Higgs_coupling.mb, vev
	)/(4.(N_CN_C - 1)(pa-pH).m2()(pb-pn).m2());
	}

	// kinetic matrix element square for backward unordered emission
	// and forward g -> Higgs
	double MatrixElement::MH2_unob_forwardH(
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & pg,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pH
	) const {
	using namespace currents;
	const double vev = param_.ew_parameters.vev();

	constexpr double K_f1 = K_q;

	constexpr double nhel = 4.;

	return K_f1*ME_juno_jgH(
	pg, p1, pa, pH, pb,
	param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
	param_.Higgs_coupling.mb, vev
	)/(nhel(N_CN_C - 1)(pa - p1 - pg).m2()(pb - pH).m2());
	}

	double MatrixElement::tree_kin_Higgs_first(Event const & ev) const {
	auto const & incoming = ev.incoming();
	auto const & outgoing = ev.outgoing();
	assert(outgoing.front().type == pid::Higgs);

	if(is_anyquark(incoming.front())) {
	assert(incoming.front().type == outgoing[1].type);
	return tree_kin_Higgs_between(ev);
	}

	const auto partons = tag_extremal_jet_partons(ev);

	const auto pa = to_HepLorentzVector(incoming.front());
	const auto pb = to_HepLorentzVector(incoming.back());

	const auto pH = to_HepLorentzVector(outgoing.front());

	const auto end_ladder = end(partons) - ((ev.type() == event_type::unof)?2:1);
	const auto pn = to_HepLorentzVector(*end_ladder);

	const double ladder = FKL_ladder_weight(
	begin(partons), end_ladder,
	pa - pH, pa, pb, pa, pn,
	param_.regulator_lambda
	);

	if(ev.type() == event_type::unof) {
	const auto pg = to_HepLorentzVector(outgoing.back());
	return MH2_unob_forwardH(
	pb, pa, pg, pn, pH
	)*ladder;
	}

	return MH2_backwardH(
	incoming.back().type,
	pa, pb, pH, pn
	)*ladder;
	}

	double MatrixElement::tree_kin_Higgs_last(Event const & ev) const {
	auto const & incoming = ev.incoming();
	auto const & outgoing = ev.outgoing();
	assert(outgoing.back().type == pid::Higgs);

	if(is_anyquark(incoming.back())) {
	assert(incoming.back().type == outgoing[outgoing.size()-2].type);
	return tree_kin_Higgs_between(ev);
	}

	const auto partons = tag_extremal_jet_partons(ev);

	const auto pa = to_HepLorentzVector(incoming.front());
	const auto pb = to_HepLorentzVector(incoming.back());

	const auto pH = to_HepLorentzVector(outgoing.back());

	auto begin_ladder = begin(partons) + 1;
	auto q0 = pa - to_HepLorentzVector(partons.front());
	if(ev.type() == event_type::unob) {
	q0 -= to_HepLorentzVector(*begin_ladder);
	++begin_ladder;
	}
	const auto p1 = to_HepLorentzVector(*(begin_ladder - 1));

	const double ladder = FKL_ladder_weight(
	begin_ladder, end(partons),
	q0, pa, pb, p1, pb,
	param_.regulator_lambda
	);

	if(ev.type() == event_type::unob) {
	const auto pg = to_HepLorentzVector(outgoing.front());
	return MH2_unob_forwardH(
	pa, pb, pg, p1, pH
	)*ladder;
	}

	return MH2_backwardH(
	incoming.front().type,
	pb, pa, pH, p1
	)*ladder;
	}

	namespace {
	template<class InIter, class partIter>
	double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
	partIter EndPart,
	CLHEP::HepLorentzVector const & qH,
	CLHEP::HepLorentzVector const & qHp1,
	double mt, bool inc_bot, double mb, double vev
	){

	const auto pa = to_HepLorentzVector(*BeginIn);
	const auto pb = to_HepLorentzVector(*(EndIn-1));

	const auto pg = to_HepLorentzVector(*BeginPart);
	const auto p1 = to_HepLorentzVector(*(BeginPart+1));
	const auto pn = to_HepLorentzVector(*(EndPart-1));

	return ME_Higgs_current_uno(
	(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa,
	qH, qHp1, mt, inc_bot, mb, vev
	);
	}
	} // namespace


	double MatrixElement::tree_kin_Higgs_between(Event const & ev) const {
	using namespace event_type;
	auto const & incoming = ev.incoming();
	auto const & outgoing = ev.outgoing();

	const auto the_Higgs = std::find_if(
	begin(outgoing), end(outgoing),
	[](Particle const & s){ return s.type == pid::Higgs; }
	);
	assert(the_Higgs != end(outgoing));
	const auto pH = to_HepLorentzVector(*the_Higgs);
	const auto partons = tag_extremal_jet_partons(ev);

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	auto p1 = to_HepLorentzVector(
	partons[(ev.type() == unob)?1:0]
	);
	auto pn = to_HepLorentzVector(
	partons[partons.size() - ((ev.type() == unof)?2:1)]
	);

	auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
	assert(
	(first_after_Higgs == end(partons) && (
	(ev.type() == unob)
	\|\| partons.back().type != pid::gluon
	))
	\|\| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
	);
	assert(
	(first_after_Higgs == begin(partons) && (
	(ev.type() == unof)
	\|\| partons.front().type != pid::gluon
	))
	\|\| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
	);
	// always treat the Higgs as if it were in between the extremal FKL partons
	if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
	else if(first_after_Higgs == end(partons)) --first_after_Higgs;

	// t-channel momentum before Higgs
	auto qH = pa;
	for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
	qH -= to_HepLorentzVector(*parton_it);
	}

	auto q0 = pa - p1;
	auto begin_ladder = begin(partons) + 1;
	auto end_ladder = end(partons) - 1;

	double current_factor = NAN;
	if(ev.type() == FKL){
	current_factor = ME_Higgs_current(
	incoming[0].type, incoming[1].type,
	pn, pb, p1, pa, qH, qH - pH,
	param_.Higgs_coupling.mt,
	param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
	param_.ew_parameters.vev()
	);
	}
	else if(ev.type() == unob){
	current_factor = 4C_AC_AC_AC_A/2*tree_kin_Higgs_uno(
	begin(incoming), end(incoming), begin(partons),
	end(partons), qH, qH-pH, param_.Higgs_coupling.mt,
	param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
	param_.ew_parameters.vev()
	);
	const auto p_unob = to_HepLorentzVector(partons.front());
	q0 -= p_unob;
	p1 += p_unob;
	++begin_ladder;
	}
	else if(ev.type() == unof){
	current_factor = 4C_AC_AC_AC_A/2*tree_kin_Higgs_uno(
	rbegin(incoming), rend(incoming), rbegin(partons),
	rend(partons), qH-pH, qH, param_.Higgs_coupling.mt,
	param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
	param_.ew_parameters.vev()
	);
	pn += to_HepLorentzVector(partons.back());
	--end_ladder;
	}
	else{
	throw std::logic_error("Can only reweight FKL or uno processes in H+Jets");
	}

	const double ladder_factor = FKL_ladder_weight(
	begin_ladder, first_after_Higgs,
	q0, pa, pb, p1, pn,
	param_.regulator_lambda
	)*FKL_ladder_weight(
	first_after_Higgs, end_ladder,
	qH - pH, pa, pb, p1, pn,
	param_.regulator_lambda
	);
	return current_factor/(4.(N_CN_C-1.))*ladder_factor;
	}

	namespace {
	double get_AWZH_coupling(Event const & ev, double alpha_s, double alpha_w) {
	std::vector<Particle> bosons = filter_AWZH_bosons(ev.outgoing());

	if(bosons.empty()) {
	return 1.;
	}

	if(bosons.size() == 1) {
	switch(bosons[0].type){
	case pid::Higgs:
	return alpha_s*alpha_s;
	case pid::Wp:
	case pid::Wm:
	return alpha_w*alpha_w;
	case pid::Z_photon_mix:
	return alpha_w*alpha_w;
	// TODO
	case pid::photon:
	case pid::Z:
	default:
	throw not_implemented("Emission of boson of unsupported type");
	}
	}

	if(bosons.size() == 2) {
	if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp) {
	return alpha_walpha_walpha_w*alpha_w;
	}
	else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm) {
	return alpha_walpha_walpha_w*alpha_w;
	}

	throw not_implemented("Emission of bosons of unsupported type");
	}


	throw not_implemented("Emission of >2 bosons is unsupported");
	}
	} // namespace

	double MatrixElement::tree_param(Event const & ev, double mur) const {
	assert(is_resummable(ev.type()));

	const auto begin_partons = ev.begin_partons();
	const auto end_partons = ev.end_partons();
	const auto num_partons = std::distance(begin_partons, end_partons);
	const double alpha_s = alpha_s_(mur);
	const double gs2 = 4.M_PIalpha_s;
	double res = std::pow(gs2, num_partons);
	if(param_.log_correction){
	// use alpha_s(q_perp), evolved to mur
	assert(num_partons >= 2);
	const auto first_emission = std::next(begin_partons);
	const auto last_emission = std::prev(end_partons);
	for(auto parton = first_emission; parton != last_emission; ++parton){
	res = 1. + alpha_s/(2.M_PI)BETA0std::log(mur/parton->perp());
	}
	}
	return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res;
	}

	} // namespace HEJ

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:34 PM (1 d, 10 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805846
Default Alt Text: (105 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions