No OneTemporary
Actions

Size

59 KB

Subscribers

None

View Options

	diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
	index d3fbad8..57df05e 100644
	--- a/src/MatrixElement.cc
	+++ b/src/MatrixElement.cc
	@@ -1,1644 +1,1707 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/MatrixElement.hh"

	#include <algorithm>
	#include <assert.h>
	#include <limits>
	#include <math.h>
	#include <stddef.h>
	#include <unordered_map>
	#include <utility>

	#include "CLHEP/Vector/LorentzVector.h"

	#include "HEJ/Constants.hh"
	#include "HEJ/Wjets.hh"
	#include "HEJ/Hjets.hh"
	#include "HEJ/jets.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/utility.hh"

	namespace HEJ{
	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_PIlog(q_j.perp2()/(lambda*lambda));
	if(! param_.log_correction) return result;
	return (
	1. + alpha_s/(4.M_PI)beta0log(murmur/(q_j.perp()*lambda))
	)*result;
	}

	Weights MatrixElement::operator()(Event const & event) const {
	return tree(event)*virtual_corrections(event);
	}

	Weights MatrixElement::tree(Event const & event) const {
	return tree_param(event)*tree_kin(event);
	}

	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;
	}

	Weights MatrixElement::virtual_corrections(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 = virtual_corrections(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 = virtual_corrections(event, var.mur);
	result.variations.emplace_back(wt);
	known.emplace(var.mur, wt);
	}
	else {
	result.variations.emplace_back(ME_it->second);
	}
	}
	return result;
	}

	double MatrixElement::virtual_corrections_W(
	Event const & event,
	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;
	size_t first_idx = 0;
	size_t last_idx = partons.size() - 1;

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

	// With extremal qqx 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::qqxexb) {
	q -= partons[1].p;
	++first_idx;
	if (abs(partons[0].type) != abs(partons[1].type)){
	q -= WBoson.p;
	#ifndef NDEBUG
	wc=false;
	#endif
	}
	}
	else {
	if(event.type() == event_type::unof
	\|\| event.type() == event_type::qqxexf){
	--last_idx;
	}
	if (in[0].type != partons[0].type ){
	q -= WBoson.p;
	#ifndef NDEBUG
	wc=false;
	#endif
	}
	}


	size_t first_idx_qqx = last_idx;
	size_t last_idx_qqx = last_idx;

	//if qqxMid event, virtual correction do not occur between
	//qqx pair.
	if(event.type() == event_type::qqxmid){
	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(abs(backquark->type) != abs((backquark+1)->type)) {
	wqq=true;
	#ifndef NDEBUG
	wc=false;
	#endif
	}
	last_idx = std::distance(begin(partons), backquark);
	first_idx_qqx = last_idx+1;
	}
	double exponent = 0;
	const double alpha_s = alpha_s_(mur);
	for(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_qqx) q -= partons[last_idx+1].p;
	if (wqq) q -= WBoson.p;

	for(size_t j = first_idx_qqx; j < last_idx_qqx; ++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::qqxexf
	);
	#endif

	return exp(exponent);
	}

	double MatrixElement::virtual_corrections(
	Event const & event,
	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

	const auto AWZH_boson = std::find_if(
	begin(out), end(out),
	[](Particle const & p){ return is_AWZH_boson(p); }
	);

	if(AWZH_boson != end(out) && abs(AWZH_boson->type) == pid::Wp){
	return virtual_corrections_W(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;
	size_t first_idx = 0;
	size_t last_idx = out.size() - 1;

	// if there is a Higgs boson, extremal qqx 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)
	\|\| event.type() == event_type::unob
	\|\| event.type() == event_type::qqxexb){
	q -= out[1].p;
	++first_idx;
	}
	if((out.back().type == pid::Higgs)
	\|\| event.type() == event_type::unof
	\|\| event.type() == event_type::qqxexf){
	--last_idx;
	}

	size_t first_idx_qqx = last_idx;
	size_t last_idx_qqx = last_idx;

	//if qqxMid event, virtual correction do not occur between
	//qqx pair.
	if(event.type() == event_type::qqxmid){
	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_qqx = last_idx+1;
	}
	double exponent = 0;
	const double alpha_s = alpha_s_(mur);
	for(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_qqx) q -= out[last_idx+1].p;

	for(size_t j = first_idx_qqx; j < last_idx_qqx; ++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
	\|\| event.type() == event_type::unof
	\|\| event.type() == event_type::qqxexf
	);
	return exp(exponent);
	}

	namespace {
	//! Lipatov vertex for partons emitted into extremal jets
	double C2Lipatov(
	CLHEP::HepLorentzVector const & qav,
	CLHEP::HepLorentzVector const & qbv,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & p2
	){
	CLHEP::HepLorentzVector temptrans=-(qav+qbv);
	CLHEP::HepLorentzVector p5=qav-qbv;
	CLHEP::HepLorentzVector CL=temptrans
	+ 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.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,
	double lambda
	) {
	double kperp=(qav-qbv).perp();
	if (kperp>lambda)
	return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());

	double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
	return Cls-4./(kperp*kperp);
	}

	//! 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
	){
	CLHEP::HepLorentzVector temptrans=-(qav+qbv);
	CLHEP::HepLorentzVector p5=qav-qbv;
	CLHEP::HepLorentzVector CL=temptrans
	+ 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.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,
	double lambda
	) {
	double kperp=(qav-qbv).perp();
	if (kperp>lambda)
	return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());

	double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
	double temp=Cls-4./(kperp*kperp);
	return temp;
	}

	/** 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(
	int aptype, int bptype,
	CLHEP::HepLorentzVector const & pg,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa
	){
	assert(aptype!=21); // aptype cannot be gluon
	if (bptype==21) {
	if (aptype > 0)
	return ME_unob_qg(pg,p1,pa,pn,pb);
	else
	return ME_unob_qbarg(pg,p1,pa,pn,pb);
	}
	else if (bptype<0) { // ----- \|\| -----
	if (aptype > 0)
	return ME_unob_qQbar(pg,p1,pa,pn,pb);
	else
	return ME_unob_qbarQbar(pg,p1,pa,pn,pb);
	}
	else { //bptype == quark
	if (aptype > 0)
	return ME_unob_qQ(pg,p1,pa,pn,pb);
	else
	return ME_unob_qbarQ(pg,p1,pa,pn,pb);
	}
	}

	/** Matrix element squared for tree-level current-current scattering
	* @param bptype Particle b PDG ID
	* @param pgin Incoming gluon momentum
	* @param pq Quark from splitting Momentum
	* @param pqbar Anti-quark from splitting Momentum
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param swap_q_qx Boolean. Ordering of qqbar pair. False: pqbar extremal.
	* @returns ME Squared for Tree-Level Current-Current Scattering
	*
	* @note The qqxf contribution can be calculated by reversing the argument ordering.
	*/
	double ME_qqx_current(
	int bptype,
	CLHEP::HepLorentzVector const & pgin,
	CLHEP::HepLorentzVector const & pq,
	CLHEP::HepLorentzVector const & pqbar,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	bool const swap_q_qx
	){
	if (bptype==21) {
	if (swap_q_qx) // pq extremal
	return ME_Exqqx_qqbarg(pgin,pq,pqbar,pn,pb);
	else // pqbar extremal
	return ME_Exqqx_qbarqg(pgin,pq,pqbar,pn,pb);
	}
	else { // b leg quark line
	if (swap_q_qx) //extremal pq
	return ME_Exqqx_qqbarQ(pgin,pq,pqbar,pn,pb);
	else
	return ME_Exqqx_qbarqQ(pgin,pq,pqbar,pn,pb);
	}
	throw std::logic_error("Unreachable in ME_Exqqx_current()");
	}

	/* \brief Matrix element squared for central qqx 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 qqxpair
	* @param nbelow Number of gluons emitted after central qqxpair
	* @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 qqxmid Tree-Level Current-Current Scattering
	*/
	double ME_qqxmid_current(
	int aptype, int bptype,
	int nabove, int nbelow,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & pq,
	CLHEP::HepLorentzVector const & pqbar,
	std::vector<HLV> const & partons){
	// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
	const bool swap_q_qx=pqbar.rapidity() < pq.rapidity();
	double wt=1.;

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

	return wt*ME_Cenqqx_qq(pa, pb, partons,(bptype<0),(aptype<0), swap_q_qx, nabove, nbelow);
	}


	/** 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(
	int aptype, int bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa
	){
	if (aptype==21&&bptype==21) {
	return ME_gg(pn,pb,p1,pa);
	} else if (aptype==21&&bptype!=21) {
	if (bptype > 0)
	return ME_qg(pn,pb,p1,pa);
	else
	return ME_qbarg(pn,pb,p1,pa);
	}
	else if (bptype==21&&aptype!=21) { // ----- \|\| -----
	if (aptype > 0)
	return ME_qg(p1,pa,pn,pb);
	else
	return ME_qbarg(p1,pa,pn,pb);
	}
	else { // they are both quark
	if (bptype>0) {
	if (aptype>0)
	return ME_qQ(pn,pb,p1,pa);
	else
	return ME_qQbar(pn,pb,p1,pa);
	}
	else {
	if (aptype>0)
	return ME_qQbar(p1,pa,pn,pb);
	else
	return ME_qbarQbar(pn,pb,p1,pa);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	/** 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(
	int aptype, int 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
	){
	// We know it cannot be gg incoming.
	assert(!(aptype==21 && bptype==21));
	if (aptype==21&&bptype!=21) {
	if (bptype > 0)
	return ME_W_qg(pn,plbar,pl,pb,p1,pa,Wprop);
	else
	return ME_W_qbarg(pn,plbar,pl,pb,p1,pa,Wprop);
	}
	else if (bptype==21&&aptype!=21) { // ----- \|\| -----
	if (aptype > 0)
	return ME_W_qg(p1,plbar,pl,pa,pn,pb,Wprop);
	else
	return ME_W_qbarg(p1,plbar,pl,pa,pn,pb,Wprop);
	}
	else { // they are both quark
	if (wc==true){ // emission off b, (first argument pbout)
	if (bptype>0) {
	if (aptype>0)
	return ME_W_qQ(pn,plbar,pl,pb,p1,pa,Wprop);
	else
	return ME_W_qQbar(pn,plbar,pl,pb,p1,pa,Wprop);
	}
	else {
	if (aptype>0)
	return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa,Wprop);
	else
	return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa,Wprop);
	}
	}
	else{ // emission off a, (first argument paout)
	if (aptype > 0) {
	if (bptype > 0)
	return ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop);
	else
	return ME_W_qQbar(p1,plbar,pl,pa,pn,pb,Wprop);
	}
	else { // a is anti-quark
	if (bptype > 0)
	return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb,Wprop);
	else
	return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb,Wprop);
	}

	}
	}
	throw std::logic_error("unknown particle types");
	}

	/** 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(
	int aptype, int 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
	){
	// we know they are not both gluons
	if (bptype == 21 && aptype != 21) { // b gluon => W emission off a
	if (aptype > 0)
	return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	else
	return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	}

	else { // they are both quark
	if (wc) {// emission off b, i.e. b is first current
	if (bptype>0){
	if (aptype>0)
	return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	else
	return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	}
	else{
	if (aptype>0)
	return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	else
	return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	}
	}
	else {// wc == false, emission off a, i.e. a is first current
	if (aptype > 0) {
	if (bptype > 0) //qq
	return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	else //qqbar
	return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	}
	else { // a is anti-quark
	if (bptype > 0) //qbarq
	return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	else //qbarqbar
	return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
	}
	}
	}
	throw std::logic_error("unknown particle types");
	}

	/** \brief Matrix element squared for backward qqx 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_q_qx Boolean. Ordering of qqbar pair. False: pqbar extremal.
	* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
	* @returns ME Squared for qqxb Tree-Level Current-Current Scattering
	*
	* @note calculate forwards qqx contribution by reversing argument ordering.
	*/
	double ME_W_qqx_current(
	int aptype, int 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_q_qx, bool const wc,
	ParticleProperties const & Wprop
	){
	// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
	const double CFbackward = K_g( (swap_q_qx)?pq:pqbar ,pa)/HEJ::C_F;

	// With qqbar we could have 2 incoming gluons and W Emission
	if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
	// This will be a wqqx emission as there is no other possible W Emission Site.
	if (swap_q_qx)
	return ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb, Wprop)*CFbackward;
	else
	return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)*CFbackward;
	}
	else if (aptype==21&&bptype!=21 ) {//a gluon => W emission off b leg or qqx
	if (!wc){ // W Emitted from backwards qqx
	if (swap_q_qx)
	return ME_WExqqx_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop)*CFbackward;
	else
	return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)*CFbackward;
	}
	else { // W Must be emitted from forwards leg.
	if (swap_q_qx)
	return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, bptype<0, Wprop)*CFbackward;
	else
	return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, bptype<0, Wprop)*CFbackward;
	}
	}
	throw std::logic_error("Incompatible incoming particle types with qqxb");
	}

	/* \brief Matrix element squared for central qqx 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 qqxpair
	* @param nbelow Number of gluons emitted after central qqxpair
	* @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 qqx
	* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
	* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
	*/
	double ME_W_qqxmid_current(
	int aptype, int bptype,
	int nabove, int nbelow,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & pq,
	CLHEP::HepLorentzVector const & pqbar,
	std::vector<HLV> const & partons,
	CLHEP::HepLorentzVector const & plbar,
	CLHEP::HepLorentzVector const & pl,
	bool const wqq, bool const wc,
	ParticleProperties const & Wprop
	){
	// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
	const bool swap_q_qx=pqbar.rapidity() < pq.rapidity();
	double wt=1.;

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

	if(wqq)
	return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons,(bptype<0),(aptype<0),
	swap_q_qx, nabove, Wprop);
	return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, (bptype<0), (aptype<0),
	swap_q_qx, nabove, nbelow, wc, Wprop);
	}

	/** \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(
	int aptype, int 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
	){
	if (aptype==21&&bptype==21) // gg initial state
	return ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev);
	else if (aptype==21&&bptype!=21) {
	if (bptype > 0)
	return ME_H_qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
	else
	return ME_H_qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
	}
	else if (bptype==21&&aptype!=21) {
	if (aptype > 0)
	return ME_H_qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
	else
	return ME_H_qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
	}
	else { // they are both quark
	if (bptype>0) {
	if (aptype>0)
	return ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)4.4./(9.*9.);
	else
	return ME_H_qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)4.4./(9.*9.);
	}
	else {
	if (aptype>0)
	return ME_H_qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)4.4./(9.*9.);
	else
	return ME_H_qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)4.4./(9.*9.);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	/** \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(
	int aptype, int 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
	){
	if (bptype==21&&aptype!=21) {
	if (aptype > 0)
	return ME_H_unob_gQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	else
	return ME_H_unob_gQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	}
	else { // they are both quark
	if (aptype>0) {
	if (bptype>0)
	return ME_H_unob_qQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	else
	return ME_H_unob_qbarQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	}
	else {
	if (bptype>0)
	return ME_H_unob_qQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	else
	return ME_H_unob_qbarQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){
	return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
	}

	void validate(HEJ::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 anonymous

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

	double MatrixElement::tree_kin(
	Event const & ev
	) const {
	if(! is_resummable(ev.type())) return 0.;

	auto AWZH_boson = std::find_if(
	begin(ev.outgoing()), end(ev.outgoing()),
	[](Particle const & p){return is_AWZH_boson(p);}
	);

	if(AWZH_boson == end(ev.outgoing()))
	return tree_kin_jets(ev);

	switch(AWZH_boson->type){
	case pid::Higgs:
	return tree_kin_Higgs(ev);
	case pid::Wp:
	case pid::Wm:
	return tree_kin_W(ev);
	// TODO
	case pid::photon:
	case pid::Z:
	default:
	throw not_implemented("Emission of boson of unsupported type");
	}
	}

	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;
	}

	} // namespace anonymous

	std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
	Event const & ev
	) const{
	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;
	}
	const auto & jets = ev.jets();
	assert(jets.size() >= 2);
	auto most_backward = begin(jets);
	auto most_forward = end(jets) - 1;
	// skip jets caused by unordered emission or qqx
	if(ev.type() == event_type::unob \|\| ev.type() == event_type::qqxexb){
	assert(jets.size() >= 3);
	++most_backward;
	}
	else if(ev.type() == event_type::unof \|\| ev.type() == event_type::qqxexf){
	assert(jets.size() >= 3);
	--most_forward;
	}
	const auto extremal_jet_indices = ev.particle_jet_indices(
	{most_backward, most_forward}
	);
	assert(extremal_jet_indices.size() == out_partons.size());
	for(size_t i = 0; i < out_partons.size(); ++i){
	assert(HEJ::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;
	}

	namespace {
	+ double tree_kin_jets_qqxmid(
	+ int aptype, int bptype, HLV pa, HLV pb,
	+ std::vector<Particle> const & partons,
	+ double lambda
	+ ){
	+ HLV pq,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 qqx
	+ auto qqxt = 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){
	+ qqxt -= to_HepLorentzVector(*parton_it);
	+ }
	+
	+ int nabove = std::distance(begin_ladder, backmidquark);
	+ int nbelow = std::distance(begin_ladder_2, end_ladder);
	+
	+ std::vector<HLV> partonsHLV;
	+ partonsHLV.reserve(partons.size());
	+ for (size_t i = 0; i != partons.size(); ++i) {
	+ partonsHLV.push_back(to_HepLorentzVector(partons[i]));
	+ }
	+
	+ const double current_factor = ME_qqxmid_current(
	+ aptype, bptype, nabove, nbelow, 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,
	+ qqxt, pa, pb, p1, pn,
	+ lambda
	+ );
	+ return current_factorC_AC_A/(N_CN_C-1.)ladder_factor;
	+ }
	+
	+
	template<class InIter, class partIter>
	double tree_kin_jets_qqx(InIter BeginIn, InIter EndIn, partIter BeginPart,
	partIter EndPart, double lambda){
	-
	const bool swap_q_qx = is_quark(*BeginPart);
	const auto pgin = to_HepLorentzVector(*BeginIn);
	const auto pb = to_HepLorentzVector(*(EndIn-1));
	const auto pq = to_HepLorentzVector(*(BeginPart+(swap_q_qx?0:1)));
	const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_q_qx?1:0)));
	const auto p1 = to_HepLorentzVector(*(BeginPart));
	const auto pn = to_HepLorentzVector(*(EndPart-1));

	assert((BeginIn)->type==21); // Incoming a must be gluon.
	const double current_factor = ME_qqx_current(
	(EndIn-1)->type, pgin, pq, pqbar, pn, pb, swap_q_qx
	)/(4.(N_CN_C - 1.));
	const double ladder_factor = FKL_ladder_weight(
	(BeginPart+2), (EndPart-1),
	pgin-pq-pqbar, 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;
	}
	}

	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()==HEJ::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
	);
	}
	else if (ev.type()==HEJ::event_type::unordered_backward){
	return tree_kin_jets_uno(incoming.begin(), incoming.end(),
	partons.begin(), partons.end(),
	param_.regulator_lambda);
	}
	else if (ev.type()==HEJ::event_type::unordered_forward){
	return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
	partons.rbegin(), partons.rend(),
	param_.regulator_lambda);
	}
	else if (ev.type()==HEJ::event_type::extremal_qqxb){
	return tree_kin_jets_qqx(incoming.begin(), incoming.end(),
	partons.begin(), partons.end(),
	param_.regulator_lambda);
	}
	else if (ev.type()==HEJ::event_type::extremal_qqxf){
	return tree_kin_jets_qqx(incoming.rbegin(), incoming.rend(),
	partons.rbegin(), partons.rend(),
	param_.regulator_lambda);
	}
	else {
	throw std::logic_error("Can only reweight FKL, uno or ExQQx processes in Pure Jets");
	}
	}

	namespace{
	double tree_kin_W_FKL(
	int aptype, int bptype, HLV pa, HLV pb,
	std::vector<Particle> const & partons,
	HLV plbar, HLV 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 HLV & plbar, const HLV & 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_factorC_AC_A/(N_CN_C-1.)ladder_factor;
	}

	template<class InIter, class partIter>
	double tree_kin_W_qqx(InIter BeginIn, partIter BeginPart,
	partIter EndPart, const HLV & plbar, const HLV & pl,
	double lambda, ParticleProperties const & Wprop){
	const bool swap_q_qx=is_quark(*BeginPart);
	const auto pa = to_HepLorentzVector(*BeginIn);
	const auto pb = to_HepLorentzVector(*(BeginIn+1));
	const auto pq = to_HepLorentzVector(*(BeginPart+(swap_q_qx?0:1)));
	const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_q_qx?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_qqx_current(
	(BeginIn)->type, (BeginIn+1)->type, pa, pb,
	pq, pqbar, pn, plbar, pl, swap_q_qx, wc, Wprop
	);

	const double ladder_factor = FKL_ladder_weight(
	BeginPart+2, EndPart-1,
	q0, pa, pb, p1, pn,
	lambda
	);
	return current_factorC_AC_A/(N_CN_C-1.)ladder_factor;
	}

	double tree_kin_W_qqxmid(
	int aptype, int bptype, HLV pa, HLV pb,
	std::vector<Particle> const & partons,
	HLV plbar, HLV pl,
	double lambda, ParticleProperties const & Wprop
	){
	HLV pq,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 qqx
	auto qqxt = q0;

	bool wc, wqq;
	if (backmidquark->type == -(backmidquark+1)->type){ // Central qqx does not emit
	wqq=false;
	if (aptype==partons[0].type) {
	wc = true;
	}
	else{
	wc = false;
	q0-=pl+plbar;
	}
	}
	else{
	wqq = true;
	wc = false;
	qqxt-=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){
	qqxt -= to_HepLorentzVector(*parton_it);
	}

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

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

	const double current_factor = ME_W_qqxmid_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,
	qqxt, pa, pb, p1, pn,
	lambda
	);
	return current_factorC_AC_A/(N_CN_C-1.)ladder_factor;
	}
	} // namespace anonymous

	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( (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
	HLV plbar, 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_qqxb){
	return tree_kin_W_qqx(cbegin(incoming), cbegin(partons),
	cend(partons), plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Wprop());
	}
	if(ev.type() == extremal_qqxf){
	return tree_kin_W_qqx(crbegin(incoming), crbegin(partons),
	crend(partons), plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Wprop());
	}
	assert(ev.type() == central_qqx);
	return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
	pa, pb, partons, plbar, pl,
	param_.regulator_lambda,
	param_.ew_parameters.Wprop());
	}

	double MatrixElement::tree_kin_Higgs(Event const & ev) const {
	if(is_uno(ev.type())){
	return tree_kin_Higgs_between(ev);
	}
	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);
	}

	namespace {
	// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
	#ifdef HEJ_BUILD_WITH_QCDLOOP
	// TODO: code duplication with jets.cc
	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());
	}
	double K(
	ParticleID type,
	CLHEP::HepLorentzVector const & pout,
	CLHEP::HepLorentzVector const & pin
	) {
	if(type == ParticleID::gluon) return K_g(pout, pin);
	return C_F;
	}
	#endif
	// Colour factor in strict MRK limit
	double K_MRK(ParticleID type) {
	return (type == ParticleID::gluon)?C_A:C_F;
	}
	}

	double MatrixElement::MH2_forwardH(
	CLHEP::HepLorentzVector const & p1out,
	CLHEP::HepLorentzVector const & p1in,
	ParticleID type2,
	CLHEP::HepLorentzVector const & p2out,
	CLHEP::HepLorentzVector const & p2in,
	CLHEP::HepLorentzVector const & pH,
	double t1, double t2
	) const{
	ignore(p2out, p2in);
	const double shat = p1in.invariantMass2(p2in);
	const double vev = param_.ew_parameters.vev();
	// gluon case
	#ifdef HEJ_BUILD_WITH_QCDLOOP
	if(!param_.Higgs_coupling.use_impact_factors){
	return K(type2, p2out, p2in)C_A1./(16M_PIM_PI)t1/t2ME_Houtside_gq(
	p1out, p1in, p2out, p2in, pH,
	param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
	param_.Higgs_coupling.mb, vev
	)/(4(N_CN_C - 1));
	}
	#endif
	return K_MRK(type2)/C_A9./2.shatshat(
	C2gHgp(p1in,p1out,pH,vev) + C2gHgm(p1in,p1out,pH,vev)
	)/(t1*t2);
	}

	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(outgoing[1].type != pid::gluon) {
	assert(incoming.front().type == outgoing[1].type);
	return tree_kin_Higgs_between(ev);
	}
	const auto pH = to_HepLorentzVector(outgoing.front());
	const auto partons = tag_extremal_jet_partons(
	ev
	);

	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());

	const auto q0 = pa - p1 - pH;

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

	return MH2_forwardH(
	p1, pa, incoming[1].type, pn, pb, pH,
	t1, t2
	)*FKL_ladder_weight(
	begin(partons) + 1, end(partons) - 1,
	q0, pa, pb, p1, pn,
	param_.regulator_lambda
	);
	}

	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(outgoing[outgoing.size()-2].type != pid::gluon) {
	assert(incoming.back().type == outgoing[outgoing.size()-2].type);
	return tree_kin_Higgs_between(ev);
	}
	const auto pH = to_HepLorentzVector(outgoing.back());
	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.front());
	const auto pn = to_HepLorentzVector(partons.back());

	auto q0 = pa - p1;

	const double t1 = q0.m2();
	const double t2 = (pn + pH - pb).m2();

	return MH2_forwardH(
	pn, pb, incoming[0].type, p1, pa, pH,
	t2, t1
	)*FKL_ladder_weight(
	begin(partons) + 1, end(partons) - 1,
	q0, pa, pb, p1, pn,
	param_.regulator_lambda
	);
	}

	namespace {
	template<class InIter, class partIter>
	double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
	partIter EndPart, const HLV & qH, const HLV & 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
	);
	}
	}


	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;
	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 = HEJ::C_AHEJ::C_A/2tree_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 = HEJ::C_AHEJ::C_A/2tree_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_factorC_AC_A/(N_CN_C-1.)ladder_factor;
	}

	namespace {
	double get_AWZH_coupling(Event const & ev, double alpha_s, double alpha_w) {
	const auto AWZH_boson = std::find_if(
	begin(ev.outgoing()), end(ev.outgoing()),
	[](auto const & p){return is_AWZH_boson(p);}
	);
	if(AWZH_boson == end(ev.outgoing())) return 1.;
	switch(AWZH_boson->type){
	case pid::Higgs:
	return alpha_s*alpha_s;
	case pid::Wp:
	case pid::Wm:
	return alpha_w*alpha_w;
	// TODO
	case pid::photon:
	case pid::Z:
	default:
	throw not_implemented("Emission of boson of unsupported type");
	}
	}
	}

	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)beta0log(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, 2:53 PM (1 d, 12 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3796594
Default Alt Text: (59 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions