No OneTemporary
Actions

Size

190 KB

Subscribers

None

View Options

	diff --git a/include/RHEJ/currents.hh b/include/RHEJ/currents.hh
	index 89ed652..789b2f1 100644
	--- a/include/RHEJ/currents.hh
	+++ b/include/RHEJ/currents.hh
	@@ -1,875 +1,883 @@
	//////////////////////////////////////////////////
	//////////////////////////////////////////////////
	// This source code is Copyright (2012) of //
	// Jeppe R. Andersen and Jennifer M. Smillie //
	// and is distributed under the //
	// Gnu Public License version 2 //
	// http://www.gnu.org/licenses/gpl-2.0.html //
	// You are allowed to distribute and alter the //
	// source under the conditions of the GPLv2 //
	// as long as this copyright notice //
	// is unaltered and distributed with the source //
	// Any use should comply with the //
	// MCNET GUIDELINES //
	// for Event Generator Authors and Users //
	// as distributed with this source code //
	//////////////////////////////////////////////////
	//////////////////////////////////////////////////

	/** \file
	* \brief Functions computing the square of current contractions.
	*
	* This file contains all the necessary functions to compute the current
	* contractions for all valid HEJ processes. PJETS, H+JETS and W+JETS along with
	* some unordered counterparts.
	*/


	#pragma once

	#include <CLHEP/Vector/LorentzVector.h>
	#include <complex>
	#include <vector>
	#include <limits>

	typedef std::complex<double> COM;
	typedef COM current[4];
	typedef CLHEP::HepLorentzVector HLV;

	//! The Higgs field vacuum expectation value in GeV
	static constexpr double v = 246.;

	constexpr double infinity = std::numeric_limits<double>::infinity();
	constexpr double mb_default = 4.7;

	void Setup_Currents(void);

	//! Square of qQ->qenuQ W+Jets Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param pe Momentum of final state electron
	* @param pnu Momentum of final state Neutrino
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state quark
	* @param p2in Momentum of intial state quark
	* @returns Square of the current contractions for qQ->qenuQ Scattering
	*
	* This returns the square of the current contractions in qQ->qenuQ scattering
	* with an emission of a W Boson.
	*/
	double jMWqQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
	CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);


	//! Square of qbarQ->qbarenuQ W+Jets Scattering Current
	/**
	* @param p1out Momentum of final state anti-quark
	* @param pe Momentum of final state electron
	* @param pnu Momentum of final state Neutrino
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state quark
	* @param p2in Momentum of intial state quark
	* @returns Square of the current contractions for qbarQ->qbarenuQ Scattering
	*
	* This returns the square of the current contractions in qbarQ->qbarenuQ scattering
	* with an emission of a W Boson.
	*/
	double jMWqbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
	CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);


	//! Square of qQbar->qenuQbar W+Jets Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param pe Momentum of final state electron
	* @param pnu Momentum of final state Neutrino
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @returns Square of the current contractions for qQbar->qenuQbar Scattering
	*
	* This returns the square of the current contractions in qQbar->qenuQbar scattering
	* with an emission of a W Boson.
	*/
	double jMWqQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
	CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);


	//! Square of qbarQbar->qbarenuQbar W+Jets Scattering Current
	/**
	* @param p1out Momentum of final state anti-quark
	* @param pe Momentum of final state electron
	* @param pnu Momentum of final state Neutrino
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @returns Square of the current contractions for qbarQbar->qbarenuQbar Scattering
	*
	* This returns the square of the current contractions in qbarQbar->qbarenuQbar scattering
	* with an emission of a W Boson.
	*/
	double jMWqbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
	CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);

	//! Square of qg->qenug W+Jets Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param pe Momentum of final state electron
	* @param pnu Momentum of final state Neutrino
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @returns Square of the current contractions for qg->qenug Scattering
	*
	* This returns the square of the current contractions in qg->qenug scattering
	* with an emission of a W Boson.
	*/
	double jMWqg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
	CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);

	//! Square of qbarg->qbarenug W+Jets Scattering Current
	/**
	* @param p1out Momentum of final state anti-quark
	* @param pe Momentum of final state electron
	* @param pnu Momentum of final state Neutrino
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @returns Square of the current contractions for qbarg->qbarenug Scattering
	*
	* This returns the square of the current contractions in qbarg->qbarenug scattering
	* with an emission of a W Boson.
	*/
	double jMWqbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
	CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);


	//! Square of qQ->qQ Pure Jets Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state quark
	* @param p2in Momentum of intial state quark
	* @returns Square of the current contractions for qQ->qQ Scattering
	*
	* This returns the square of the current contractions in qQ->qQ Pure Jet Scattering.
	*/
	double jM2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);


	//! Square of qQbar->qQbar Pure Jets Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @returns Square of the current contractions for qQbar->qQbar Scattering
	*
	* This returns the square of the current contractions in qQbar->qQbar Pure Jet Scattering.
	* Note this can be used for qbarQ->qbarQ Scattering by inputting arguments appropriately.
	*/
	double jM2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);


	//! Square of qbarQbar->qbarQbar Pure Jets Scattering Current
	/**
	* @param p1out Momentum of final state anti-quark
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
	*
	* This returns the square of the current contractions in qbarQbar->qbarQbar Pure Jet Scattering.
	*/
	double jM2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);


	//! Square of qg->qg Pure Jets Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @returns Square of the current contractions for qg->qg Scattering
	*
	* This returns the square of the current contractions in qg->qg Pure Jet Scattering.
	* Note this can be used for gq->gq Scattering by inputting arguments appropriately.
	*/
	double jM2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);


	//! Square of qbarg->qbarg Pure Jets Scattering Current
	/**
	* @param p1out Momentum of final state anti-quark
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @returns Square of the current contractions for qbarg->qbarg Scattering
	*
	* This returns the square of the current contractions in qbarg->qbarg Pure Jet Scattering.
	* Note this can be used for gqbar->gqbar Scattering by inputting arguments appropriately.
	*/
	double jM2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);


	//! Square of gg->gg Pure Jets Scattering Current
	/**
	* @param p1out Momentum of final state gluon
	* @param p1in Momentum of initial state gluon
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @returns Square of the current contractions for gg->gg Scattering
	*
	* This returns the square of the current contractions in gg->gg Pure Jet Scattering.
	*/
	double jM2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);


	//! Square of gg->gg Higgs+Jets Scattering Current
	/**
	* @param p1out Momentum of final state gluon
	* @param p1in Momentum of initial state gluon
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @param q1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for gg->gg Scattering
	*
	* This returns the square of the current contractions in gg->gg Higgs+Jet Scattering.
	*
	* g~p1 g~p2
	* should be called with q1 meant to be contracted with p2 in first part of vertex
	* (i.e. if g is backward, q1 is forward)
	*/
	double MH2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);

	//! Square of gq->gq Higgs+Jets Scattering Current with Higgs before Gluon
	/**
	* @param p1out Momentum of final state gluon
	* @param p1in Momentum of initial state gluon
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @param pH Momentum of Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contraction
	*
	*/
	double MH2gq_outsideH(CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector pH,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//! Square of qg->qg Higgs+Jets Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @param q1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qg->qg Scattering
	*
	* This returns the square of the current contractions in qg->qg Higgs+Jet Scattering.
	*
	* q~p1 g~p2 (i.e. ALWAYS p1 for quark, p2 for gluon)
	* should be called with q1 meant to be contracted with p2 in first part of vertex
	* (i.e. if g is backward, q1 is forward)
	*/
	double MH2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//! Square of qbarg->qbarg Higgs+Jets Scattering Current
	/**
	* @param p1out Momentum of final state anti-quark
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @param q1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qbarg->qbarg Scattering
	*
	* This returns the square of the current contractions in qbarg->qbarg Higgs+Jet Scattering.
	*
	* qbar~p1 g~p2 (i.e. ALWAYS p1 for anti-quark, p2 for gluon)
	* should be called with q1 meant to be contracted with p2 in first part of vertex
	* (i.e. if g is backward, q1 is forward)
	*/
	double MH2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//! Square of qQ->qQ Higgs+Jets Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state quark
	* @param p2in Momentum of intial state quark
	* @param q1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qQ->qQ Scattering
	*
	* This returns the square of the current contractions in qQ->qQ Higgs+Jet Scattering.
	*
	* q~p1 Q~p2 (i.e. ALWAYS p1 for quark, p2 for quark)
	* should be called with q1 meant to be contracted with p2 in first part of vertex
	* (i.e. if Q is backward, q1 is forward)
	*/
	double MH2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);

	//! Square of qQbar->qQbar Higgs+Jets Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @param q1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qQ->qQ Scattering
	*
	* This returns the square of the current contractions in qQbar->qQbar Higgs+Jet Scattering.
	*
	* q~p1 Qbar~p2 (i.e. ALWAYS p1 for quark, p2 for anti-quark)
	* should be called with q1 meant to be contracted with p2 in first part of vertex
	* (i.e. if Qbar is backward, q1 is forward)
	*/
	double MH2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);

	//! Square of qbarQ->qbarQ Higgs+Jets Scattering Current
	/**
	* @param p1out Momentum of final state anti-quark
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state quark
	* @param p2in Momentum of intial state quark
	* @param q1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qbarQ->qbarQ Scattering
	*
	* This returns the square of the current contractions in qbarQ->qbarQ Higgs+Jet Scattering.
	*
	* qbar~p1 Q~p2 (i.e. ALWAYS p1 for anti-quark, p2 for quark)
	* should be called with q1 meant to be contracted with p2 in first part of vertex
	* (i.e. if Q is backward, q1 is forward)
	*/
	double MH2qbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//! Square of qbarQbar->qbarQbar Higgs+Jets Scattering Current
	/**
	* @param p1out Momentum of final state anti-quark
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @param q1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
	*
	* This returns the square of the current contractions in qbarQbar->qbarQbar Higgs+Jet Scattering.
	*
	* qbar~p1 Qbar~p2 (i.e. ALWAYS p1 for anti-quark, p2 for anti-quark)
	* should be called with q1 meant to be contracted with p2 in first part of vertex
	* (i.e. if Qbar is backward, q1 is forward)
	*/
	double MH2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	// Unordered f

	//! Square of qQ->gqQ Higgs+Jets Unordered f Scattering Current
	/**
	* @param pg Momentum of unordered gluon
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state quark
	* @param p2in Momentum of intial state quark
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qQ->gqQ Scattering
	*
	* This returns the square of the current contractions in qQ->gqQ Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: pg > p1out >> p2out
	*/
	double jM2unogqHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
	CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//! Square of qQbar->gqQbar Higgs+Jets Unordered f Scattering Current
	/**
	* @param pg Momentum of unordered gluon
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qQbar->gqQbar Scattering
	*
	* This returns the square of the current contractions in qQbar->gqQbar Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: pg > p1out >> p2out
	*/
	double jM2unogqHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
	CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);

	//! Square of qbarQ->gqbarQ Higgs+Jets Unordered f Scattering Current
	/**
	* @param pg Momentum of unordered gluon
	* @param p1out Momentum of final state anti-quark
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state quark
	* @param p2in Momentum of intial state quark
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qbarQ->gqbarQ Scattering
	*
	* This returns the square of the current contractions in qbarQ->gqbarQ Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: pg > p1out >> p2out
	*/
	double jM2unogqbarHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
	CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);

	//! Square of qbarQbar->gqbarQbar Higgs+Jets Unordered f Scattering Current
	/**
	* @param pg Momentum of unordered gluon
	* @param p1out Momentum of final state anti-quark
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qbarQbar->gqbarQbar Scattering
	*
	* This returns the square of the current contractions in qbarQbar->gqbarQbar Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: pg > p1out >> p2out
	*/
	double jM2unogqbarHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
	CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//! Square of qg->gqg Higgs+Jets Unordered f Scattering Current
	/**
	* @param pg Momentum of unordered gluon
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qg->gqg Scattering
	*
	* This returns the square of the current contractions in qg->gqg Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: pg > p1out >> p2out
	*/
	double jM2unogqHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
	CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//! Square of qbarg->gqbarg Higgs+Jets Unordered f Scattering Current
	/**
	* @param pg Momentum of unordered gluon
	* @param p1out Momentum of final state anti-quark
	* @param p1in Momentum of initial state anti-quark
	* @param p2out Momentum of final state gluon
	* @param p2in Momentum of intial state gluon
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qbarg->gbarg Scattering
	*
	* This returns the square of the current contractions in qbarg->gqbarg Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: pg > p1out >> p2out
	*/
	double jM2unogqbarHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
	CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//Unordered b

	//! Square of qbarQ->qbarQg Higgs+Jets Unordered b Scattering Current
	/**
	* @param p1out Momentum of final state anti-quark
	* @param p1in Momentum of initial state anti-quark
	* @param pg Momentum of unordered b gluon
	* @param p2out Momentum of final state quark
	* @param p2in Momentum of intial state quark
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qbarQ->qbarQg Scattering
	*
	* This returns the square of the current contractions in qbarQ->qbarQg Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: p1out >> p2out > pg
	*/
	double jM2unobqbarHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//! Square of qQ->qQg Higgs+Jets Unordered b Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param pg Momentum of unordered b gluon
	* @param p2out Momentum of final state quark
	* @param p2in Momentum of intial state quark
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qQ->qQg Scattering
	*
	* This returns the square of the current contractions in qQ->qQg Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: p1out >> p2out > pg
	*/
	double jM2unobqHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//! Square of qQbar->qQbarg Higgs+Jets Unordered b Scattering Current
	/**
	* @param p1out Momentum of final state quark
	* @param p1in Momentum of initial state quark
	* @param pg Momentum of unordered b gluon
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qQbar->qQbarg Scattering
	*
	* This returns the square of the current contractions in qQbar->qQbarg Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: p1out >> p2out > pg
	*/
	double jM2unobqHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);

	//! Square of qbarQbar->qbarQbarg Higgs+Jets Unordered b Scattering Current
	/**
	* @param p1out Momentum of final state anti-quark
	* @param p1in Momentum of initial state anti-quark
	* @param pg Momentum of unordered b gluon
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for qbarQbar->qbarQbarg Scattering
	*
	* This returns the square of the current contractions in qbarQbar->qbarQbarg Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: p1out >> p2out > pg
	*/
	double jM2unobqbarHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);


	//! Square of gQbar->gQbarg Higgs+Jets Unordered b Scattering Current
	/**
	* @param p1out Momentum of final state gluon
	* @param p1in Momentum of initial state gluon
	* @param pg Momentum of unordered b gluon
	* @param p2out Momentum of final state anti-quark
	* @param p2in Momentum of intial state anti-quark
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for gQbar->gQbarg Scattering
	*
	* This returns the square of the current contractions in gQbar->gQbarg Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: p1out >> p2out > pg
	*/
	double jM2unobgHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);

	//! Square of gQ->gQg Higgs+Jets Unordered b Scattering Current
	/**
	* @param p1out Momentum of final state gluon
	* @param p1in Momentum of initial state gluon
	* @param pg Momentum of unordered b gluon
	* @param p2out Momentum of final state quark
	* @param p2in Momentum of intial state quark
	* @param qH1 Momentum of t-channel propagator before Higgs
	* @param qH2 Momentum of t-channel propagator after Higgs
	* @param mt Top quark mass
	* @param include_bottom Specifies whether bottom corrections are included
	* @param mb Bottom quark mass
	* @returns Square of the current contractions for gQ->gQg Scattering
	*
	* This returns the square of the current contractions in gQ->gQg Higgs+Jet Scattering.
	*
	* This construction is taking rapidity order: p1out >> p2out > pg
	*/
	double jM2unobgHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
	CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
	CLHEP::HepLorentzVector qH2,
	double mt = infinity,
	bool include_bottom = false, double mb = mb_default);

	// impact factors for Higgs + jet


	//! Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta
	/**
	* @param p2 Momentum of Particle 2
	* @param p1 Momentum of Particle 1
	* @param pH Momentum of Higgs
	* @returns Value of Eq. (4.22) in Hep-ph/0301013 with modifications
	*
	* This gives the impact factor. First it determines first whether this is the case
	* p1p\sim php>>p3p or the opposite
	*/
	double C2gHgm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1,
	CLHEP::HepLorentzVector pH);


	//! Implements Eq. (4.23) in hep-ph/0301013 with modifications to incoming plus momenta
	/**
	* @param p2 Momentum of Particle 2
	* @param p1 Momentum of Particle 1
	* @param pH Momentum of Higgs
	* @returns Value of Eq. (4.23) in Hep-ph/0301013
	*
	* This gives the impact factor. First it determines first whether this is the case
	* p1p\sim php>>p3p or the opposite
	*/
	double C2gHgp(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1,
	CLHEP::HepLorentzVector pH);


	//! Implements Eq. (4.22) in hep-ph/0301013
	/**
	* @param p2 Momentum of Particle 2
	* @param p1 Momentum of Particle 1
	* @param pH Momentum of Higgs
	* @returns Value of Eq. (4.22) in Hep-ph/0301013
	*
	* This gives the impact factor. First it determines first whether this is the case
	* p1p\sim php>>p3p or the opposite
	*/
	double C2qHqm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1,
	CLHEP::HepLorentzVector pH);

	/** \class CCurrent currents.hh "include/RHEJ/currents.hh"
	* \brief This is the a new class structure for currents.
	*/
	class CCurrent
	{
	public:
	CCurrent(COM sc0, COM sc1, COM sc2, COM sc3)
	:c0(sc0),c1(sc1),c2(sc2),c3(sc3)
	{};
	CCurrent(const CLHEP::HepLorentzVector p)
	{
	c0=p.e();
	c1=p.px();
	c2=p.py();
	c3=p.pz();
	};
	CCurrent()
	{};
	CCurrent operator+(const CCurrent& other);
	CCurrent operator-(const CCurrent& other);
	CCurrent operator*(const double x);
	CCurrent operator*(const COM x);
	CCurrent operator/(const double x);
	CCurrent operator/(const COM x);

	friend std::ostream& operator<<(std::ostream& os, const CCurrent& cur);
	COM dot(CLHEP::HepLorentzVector p1);
	COM dot(CCurrent p1);
	COM c0,c1,c2,c3;
	private:
	};

	/* std::ostream& operator <<(std::ostream& os, const CCurrent& cur); */
	CCurrent operator * ( double x, CCurrent& m);
	CCurrent operator * ( COM x, CCurrent& m);
	CCurrent operator / ( double x, CCurrent& m);
	CCurrent operator / ( COM x, CCurrent& m);

	-//! Current ???
	+//! Current <outgoing state \| mu \| incoming state>
	/**
	* These functions are a mess. There are many more defined in the source file than declared in the
	* header - and the arguments are mislabelled in some cases. Need to investigate.
	*/
	+//! @TODO remove
	CCurrent j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin);

	+//! Current <outgoing state \| mu \| incoming state>
	+/**
	+ * These functions are a mess. There are many more defined in the source file than declared in the
	+ * header - and the arguments are mislabelled in some cases. Need to investigate.
	+ */
	+CCurrent joi (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin);
	+
	//! Current <incoming state \| mu \| outgoing state>
	/**
	* These functions are a mess. There are many more defined in the source file than declared in the
	* header - and the arguments are mislabelled in some cases. Need to investigate.
	*/
	CCurrent jio (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin);

	//! Current <outgoing state \| mu \| outgoing state>
	/**
	* These functions are a mess. There are many more defined in the source file than declared in the
	* header - and the arguments are mislabelled in some cases. Need to investigate.
	*/
	CCurrent joo (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin);

	/* // Coupling values */
	/* const double stw2 = 0.2222; */
	/* const double ctw = sqrt(1.0 - stw2); */
	/* const double gs = 1.217716; */
	/* const double gw = 0.653232911; */
	/* const double Zem = (-1.0 / 2.0 + stw2) / ctw; */
	/* const double Zep = stw2 / ctw; */
	/* const double Zum = ( 1.0 / 2.0 - 2.0 * stw2 / 3.0) / ctw; */
	/* const double Zup = - 2.0 * stw2 / 3.0 / ctw; */
	/* const double Zdm = (-1.0 / 2.0 + 1.0 / 3.0 * stw2) / ctw; */
	/* const double Zdp = stw2 / 3.0 / ctw; */
	/* const double RWeak = -pow(gw, 2.0); */
	/* const double Strong = pow(gs, 4.0); */
	/* const double ee = pow(gw, 2.0) * stw2; */


	/* std::vector <double> jMZqQ (HLV, HLV, HLV, HLV, HLV, HLV, std::vector <double>, std::vector < std::vector <double> >, int, int, bool, bool); */
	/* std::vector <double> jMZqg (HLV, HLV, HLV, HLV, HLV, HLV, std::vector <double>, std::vector < std::vector <double> >, int, int, bool, bool); */
	/* void jZ (HLV, HLV, HLV, HLV, bool, bool, current); */
	/* void jZbar (HLV, HLV, HLV, HLV, bool, bool, current); */
	/* COM PZ(double); */
	/* double Zq (int, bool); */
	/* double Gq (int); */
	diff --git a/src/currents.cc b/src/currents.cc
	index 00099f7..e5f9522 100644
	--- a/src/currents.cc
	+++ b/src/currents.cc
	@@ -1,3767 +1,3615 @@
	//////////////////////////////////////////////////
	//////////////////////////////////////////////////
	// This source code is Copyright (2012) of //
	// Jeppe R. Andersen and Jennifer M. Smillie //
	// and is distributed under the //
	// Gnu Public License version 2 //
	// http://www.gnu.org/licenses/gpl-2.0.html //
	// You are allowed to distribute and alter the //
	// source under the conditions of the GPLv2 //
	// as long as this copyright notice //
	// is unaltered and distributed with the source //
	// Any use should comply with the //
	// MCNET GUIDELINES //
	// for Event Generator Authors and Users //
	// as distributed with this source code //
	//////////////////////////////////////////////////
	//////////////////////////////////////////////////
	#include "RHEJ/currents.hh"
	//#include "ZJets/Flags.h"
	#include "RHEJ/Constants.hh"
	#include "RHEJ/utility.hh"
	#include "RHEJ/PDG_codes.hh"

	const COM looprwfactor = (COM(0.,1.)M_PIM_PI)/pow((2.*M_PI),4);
	//const double HVE = 246.21845810181637;

	#ifdef RHEJ_BUILD_WITH_QCDLOOP
	#include "qcdloop/qcdloop.h"
	#endif

	#include <iostream>
	+#include <utility>

	namespace {
	// Loop integrals
	#ifdef RHEJ_BUILD_WITH_QCDLOOP

	COM B0DD(CLHEP::HepLorentzVector 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(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector 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];
	}
	COM D0DD(CLHEP::HepLorentzVector q1,CLHEP::HepLorentzVector q2, CLHEP::HepLorentzVector q3, double mq)
	{
	static std::vector<std::complex<double>> result(3);
	static auto ql_D0 = [](){
	ql::Box<std::complex<double>,double,double> ql_D0;
	ql_D0.setCacheSize(100);
	return ql_D0;
	}();
	static std::vector<double> masses(4);
	static std::vector<double> momenta(6);
	for(auto & m: masses) m = mq*mq;
	momenta[0] = q1.m2();
	momenta[1] = q2.m2();
	momenta[2] = q3.m2();
	momenta[3] = (q1+q2+q3).m2();
	momenta[4] = (q1+q2).m2();
	momenta[5] = (q2+q3).m2();
	ql_D0.integral(result, 1, masses, momenta);
	return result[0];
	}

	COM A1(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt)
	// As given in Eq. (B.2) of VDD
	{
	double q12,q22,Q2;
	CLHEP::HepLorentzVector Q;
	double Delta3,mt2;
	COM ans(COM(0.,0.));

	q12=q1.m2();
	q22=q2.m2();
	Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
	Q2=Q.m2();

	Delta3=q12q12+q22q22+Q2Q2-2q12q22-2q12Q2-2q22*Q2;
	if (mt < 0.)
	std::cerr<<"Problem in A1! mt = "<<mt<<std::endl;
	mt2=mt*mt;

	ans=looprwfactorCOM(0,-1)C0DD(q1,q2,mt)( 4.mt2/Delta3*(Q2-q12-q22)
	-1.-4.q12q22/Delta3-12.q12q22Q2/Delta3/Delta3(q12+q22-Q2) )
	- looprwfactorCOM(0,-1)( B0DD(q2,mt)-B0DD(Q,mt) )
	* ( 2.q22/Delta3+12.q12q22/Delta3/Delta3(q22-q12+Q2) )
	- looprwfactorCOM(0,-1)( B0DD(q1,mt)-B0DD(Q,mt) )
	* ( 2.q12/Delta3+12.q12q22/Delta3/Delta3(q12-q22+Q2) )
	- 2./Delta3/16/M_PI/M_PI*(q12+q22-Q2);

	return ans;

	}


	COM A2(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt)
	// As given in Eq. (B.2) of VDD, but with high energy limit
	// of invariants taken.
	{
	double q12,q22,Q2;
	CLHEP::HepLorentzVector Q;
	double Delta3,mt2;
	COM ans(COM(0.,0.));

	if (mt < 0.)
	std::cerr<<"Problem in A2! mt = "<<mt<<std::endl;
	mt2=mt*mt;

	q12=q1.m2();
	q22=q2.m2();
	Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
	Q2=Q.m2();

	Delta3=q12q12+q22q22+Q2Q2-2q12q22-2q12Q2-2q22*Q2;
	ans=looprwfactorCOM(0,-1)C0DD(q1,q2,mt)( 2.mt2+1./2.*(q12+q22-Q2)
	+2.q12q22*Q2/Delta3 )
	+looprwfactorCOM(0,-1)(B0DD(q2,mt)-B0DD(Q,mt))
	q22(q22-q12-Q2)/Delta3
	+looprwfactorCOM(0,-1)(B0DD(q1,mt)-B0DD(Q,mt))
	q12(q12-q22-Q2)/Delta3+1./16/M_PI/M_PI;

	return ans;
	}

	#else // no QCDloop

	COM A1(CLHEP::HepLorentzVector, CLHEP::HepLorentzVector, double) {
	throw std::logic_error{"A1 called without QCDloop support"};
	}

	COM A2(CLHEP::HepLorentzVector, CLHEP::HepLorentzVector, double) {
	throw std::logic_error{"A2 called without QCDloop support"};
	}

	#endif

	void to_current(const CLHEP::HepLorentzVector & q, current & ret){
	ret[0]=q.e();
	ret[1]=q.x();
	ret[2]=q.y();
	ret[3]=q.z();
	}

	constexpr double C_A = 3.;
	constexpr double C_F = 4./3.;
	using ParticleID = RHEJ::pid::ParticleID;

	// Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
	// TODO: this is not a current and should be moved somewhere else
	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());
	}

	} // namespace anonymous

	/// @TODO move all of these functions to header? clean up in general
	COM cdot(const current & j1, const current & j2)
	{
	return j1[0]j2[0]-j1[1]j2[1]-j1[2]j2[2]-j1[3]j2[3];
	}

	COM cdot(const HLV & p, const current & j1) {
	return j1[0]p.e()-j1[1]p.x()-j1[2]p.y()-j1[3]p.z();
	}

	void cmult(const COM & factor, const current & j1, current &cur)
	{
	cur[0]=factor*j1[0];
	cur[1]=factor*j1[1];
	cur[2]=factor*j1[2];
	cur[3]=factor*j1[3];
	}

	// WHY!?!
	void cadd(const current & j1, const current & j2, const current & j3,
	const current & j4, const current & j5, current &sum)
	{
	sum[0]=j1[0]+j2[0]+j3[0]+j4[0]+j5[0];
	sum[1]=j1[1]+j2[1]+j3[1]+j4[1]+j5[1];
	sum[2]=j1[2]+j2[2]+j3[2]+j4[2]+j5[2];
	sum[3]=j1[3]+j2[3]+j3[3]+j4[3]+j5[3];
	}

	void cadd(const current & j1, const current & j2, const current & j3,
	const current & j4, current &sum) {
	sum[0] = j1[0] + j2[0] + j3[0] + j4[0];
	sum[1] = j1[1] + j2[1] + j3[1] + j4[1];
	sum[2] = j1[2] + j2[2] + j3[2] + j4[2];
	sum[3] = j1[3] + j2[3] + j3[3] + j4[3];
	}

	void cadd(const current & j1, const current & j2, const current & j3,
	current &sum)
	{
	sum[0]=j1[0]+j2[0]+j3[0];
	sum[1]=j1[1]+j2[1]+j3[1];
	sum[2]=j1[2]+j2[2]+j3[2];
	sum[3]=j1[3]+j2[3]+j3[3];
	}

	void cadd(const current & j1, const current & j2, current &sum)
	{
	sum[0]=j1[0]+j2[0];
	sum[1]=j1[1]+j2[1];
	sum[2]=j1[2]+j2[2];
	sum[3]=j1[3]+j2[3];
	}

	double abs2(const COM & a)
	{
	return (a*conj(a)).real();
	}

	double vabs2(const CCurrent & cur)
	{
	return abs2(cur.c0)-abs2(cur.c1)-abs2(cur.c2)-abs2(cur.c3);
	}

	double vre(const CCurrent & a, const CCurrent & b)
	{
	return real(a.c0conj(b.c0)-a.c1conj(b.c1)-a.c2conj(b.c2)-a.c3conj(b.c3));
	}

	CCurrent CCurrent::operator+(const CCurrent& other)
	{
	COM result_c0=c0 + other.c0;
	COM result_c1=c1 + other.c1;
	COM result_c2=c2 + other.c2;
	COM result_c3=c3 + other.c3;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	CCurrent CCurrent::operator-(const CCurrent& other)
	{
	COM result_c0=c0 - other.c0;
	COM result_c1=c1 - other.c1;
	COM result_c2=c2 - other.c2;
	COM result_c3=c3 - other.c3;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	CCurrent CCurrent::operator*(const double x)
	{
	COM result_c0=x*CCurrent::c0;
	COM result_c1=x*CCurrent::c1;
	COM result_c2=x*CCurrent::c2;
	COM result_c3=x*CCurrent::c3;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	CCurrent CCurrent::operator/(const double x)
	{
	COM result_c0=CCurrent::c0/x;
	COM result_c1=CCurrent::c1/x;
	COM result_c2=CCurrent::c2/x;
	COM result_c3=CCurrent::c3/x;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	CCurrent CCurrent::operator*(const COM x)
	{
	COM result_c0=x*CCurrent::c0;
	COM result_c1=x*CCurrent::c1;
	COM result_c2=x*CCurrent::c2;
	COM result_c3=x*CCurrent::c3;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	CCurrent CCurrent::operator/(const COM x)
	{
	COM result_c0=(CCurrent::c0)/x;
	COM result_c1=(CCurrent::c1)/x;
	COM result_c2=(CCurrent::c2)/x;
	COM result_c3=(CCurrent::c3)/x;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	std::ostream& operator <<(std::ostream& os, const CCurrent& cur)
	{
	os << "("<<cur.c0<< " ; "<<cur.c1<<" , "<<cur.c2<<" , "<<cur.c3<<")";
	return os;
	}

	CCurrent operator * ( double x, CCurrent& m)
	{
	return m*x;
	}

	CCurrent operator * ( COM x, CCurrent& m)
	{
	return m*x;
	}

	CCurrent operator / ( double x, CCurrent& m)
	{
	return m/x;
	}

	CCurrent operator / ( COM x, CCurrent& m)
	{
	return m/x;
	}

	COM CCurrent::dot(CLHEP::HepLorentzVector p1)
	{
	// Current goes (E,px,py,pz)
	// std::cout<<"current = ("<<c0<<","<<c1<<","<<c2<<","<<c3<<")\n";
	// Vector goes (px,py,pz,E)
	// std::cout<<"vector = ("<<p1[0]<<","<<p1[1]<<","<<p1[2]<<","<<p1[3]<<")\n";
	return p1[3]c0-p1[0]c1-p1[1]c2-p1[2]c3;
	}

	COM CCurrent::dot(CCurrent p1)
	{
	return p1.c0c0-p1.c1c1-p1.c2c2-p1.c3c3;
	}

	+//Current Functions

	-
	-void j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin,current &cur) {
	-
	+// Current for <outgoing state \| mu \| incoming state>
	+/// @TODO always use this instead of "j"
	+/// @TODO isn't this jio with flipt helicities?
	+void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) {
	cur[0]=0.;
	cur[1]=0.;
	cur[2]=0.;
	cur[3]=0.;

	- double sqpop=sqrt(pout.plus());
	- double sqpom=sqrt(pout.minus());
	- COM poperp=pout.x()+COM(0,1)*pout.y();
	-
	- if (helout!=helin) {
	- std::cerr<< "void j : Non-matching helicities at line " << __LINE__ << std::endl;
	- } else if (helout==false) { // negative helicity
	- if (pin.plus()>pin.minus()) { // if forward
	- double sqpip=sqrt(pin.plus());
	- cur[0]=sqpop*sqpip;
	- cur[1]=sqpomsqpippoperp/abs(poperp);
	- cur[2]=-COM(0,1)*cur[1];
	- cur[3]=cur[0];
	+ const double sqpop = sqrt(pout.plus());
	+ const double sqpom = sqrt(pout.minus());
	+ const COM poperp = pout.x() + COM(0, 1) * pout.y();
	+
	+ if (helout != helin) {
	+ throw std::invalid_argument{"Non-matching helicities"};
	+ } else if (helout == false) { // negative helicity
	+ if (pin.plus() > pin.minus()) { // if forward
	+ const double sqpip = sqrt(pin.plus());
	+ cur[0] = sqpop * sqpip;
	+ cur[1] = sqpom * sqpip * poperp / abs(poperp);
	+ cur[2] = -COM(0,1) * cur[1];
	+ cur[3] = cur[0];
	} else { // if backward
	- double sqpim=sqrt(pin.minus());
	- cur[0]=-sqpomsqpimpoperp/abs(poperp);
	- cur[1]=-sqpim*sqpop;
	- cur[2]=COM(0,1)*cur[1];
	- cur[3]=-cur[0];
	+ const double sqpim = sqrt(pin.minus());
	+ cur[0] = -sqpom * sqpim * poperp / abs(poperp);
	+ cur[1] = -sqpim * sqpop;
	+ cur[2] = COM(0,1) * cur[1];
	+ cur[3] = -cur[0];
	}
	} else { // positive helicity
	- if (pin.plus()>pin.minus()) { // if forward
	- double sqpip=sqrt(pin.plus());
	- cur[0]=sqpop*sqpip;
	- cur[1]=sqpomsqpipconj(poperp)/abs(poperp);
	- cur[2]=COM(0,1)*cur[1];
	- cur[3]=cur[0];
	+ if (pin.plus() > pin.minus()) { // if forward
	+ const double sqpip = sqrt(pin.plus());
	+ cur[0] = sqpop * sqpip;
	+ cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp);
	+ cur[2] = COM(0,1) * cur[1];
	+ cur[3] = cur[0];
	} else { // if backward
	- double sqpim=sqrt(pin.minus());
	- cur[0]=-sqpomsqpimconj(poperp)/abs(poperp);
	- cur[1]=-sqpim*sqpop;
	- cur[2]=-COM(0,1)*cur[1];
	- cur[3]=-cur[0];
	+ const double sqpim = sqrt(pin.minus());
	+ cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp);
	+ cur[1] = -sqpim * sqpop;
	+ cur[2] = -COM(0,1) * cur[1];
	+ cur[3] = -cur[0];
	}
	}
	}

	-CCurrent j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin)
	+CCurrent joi (HLV pout, bool helout, HLV pin, bool helin)
	{
	- COM cur[4];
	-
	- cur[0]=0.;
	- cur[1]=0.;
	- cur[2]=0.;
	- cur[3]=0.;
	-
	- double sqpop=sqrt(pout.plus());
	- double sqpom=sqrt(pout.minus());
	- COM poperp=pout.x()+COM(0,1)*pout.y();
	-
	- if (helout!=helin) {
	- std::cerr<< "void j : Non-matching helicities\n";
	- } else if (helout==false) { // negative helicity
	- if (pin.plus()>pin.minus()) { // if forward
	- double sqpip=sqrt(pin.plus());
	- cur[0]=sqpop*sqpip;
	- cur[1]=sqpomsqpippoperp/abs(poperp);
	- cur[2]=-COM(0,1)*cur[1];
	- cur[3]=cur[0];
	- } else { // if backward
	- double sqpim=sqrt(pin.minus());
	- cur[0]=-sqpomsqpimpoperp/abs(poperp);
	- cur[1]=-sqpim*sqpop;
	- cur[2]=COM(0,1)*cur[1];
	- cur[3]=-cur[0];
	- }
	- } else { // positive helicity
	- if (pin.plus()>pin.minus()) { // if forward
	- double sqpip=sqrt(pin.plus());
	- cur[0]=sqpop*sqpip;
	- cur[1]=sqpomsqpipconj(poperp)/abs(poperp);
	- cur[2]=COM(0,1)*cur[1];
	- cur[3]=cur[0];
	- } else { // if backward
	- double sqpim=sqrt(pin.minus());
	- cur[0]=-sqpomsqpimconj(poperp)/abs(poperp);
	- cur[1]=-sqpim*sqpop;
	- cur[2]=-COM(0,1)*cur[1];
	- cur[3]=-cur[0];
	- }
	- }
	- CCurrent temp(cur[0],cur[1],cur[2],cur[3]);
	- return temp;
	+ current cur;
	+ joi(pout, helout, pin, helin, cur);
	+ return CCurrent(cur[0],cur[1],cur[2],cur[3]);
	+}
	+
	+/// @TODO remove this
	+void j (HLV pout, bool helout, HLV pin, bool helin,current &cur) {
	+ joi(pout, helout, pin, helin, cur);
	}

	-CCurrent jio (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout, bool helout)
	+/// @TODO remove this
	+CCurrent j (HLV pout, bool helout, HLV pin, bool helin)
	{
	- COM cur[4];
	-
	- cur[0]=0.;
	- cur[1]=0.;
	- cur[2]=0.;
	- cur[3]=0.;
	-
	- double sqpop=sqrt(pout.plus());
	- double sqpom=sqrt(pout.minus());
	- COM poperp=pout.x()+COM(0,1)*pout.y();
	-
	- if (helout!=helin) {
	- std::cerr<< "void j : Non-matching helicities\n";
	- } else if (helout==false) { // negative helicity
	- if (pin.plus()>pin.minus()) { // if forward
	- double sqpip=sqrt(pin.plus());
	- cur[0]=sqpop*sqpip;
	- cur[1]=sqpomsqpipconj(poperp)/abs(poperp);
	- cur[2]=COM(0,1)*cur[1];
	- cur[3]=cur[0];
	- } else { // if backward
	- double sqpim=sqrt(pin.minus());
	- cur[0]=-sqpomsqpimconj(poperp)/abs(poperp);
	- cur[1]=-sqpim*sqpop;
	- cur[2]=-COM(0,1)*cur[1];
	- cur[3]=-cur[0];
	- }
	- } else { // positive helicity
	- if (pin.plus()>pin.minus()) { // if forward
	- double sqpip=sqrt(pin.plus());
	- cur[0]=sqpop*sqpip;
	- cur[1]=sqpomsqpippoperp/abs(poperp);
	- cur[2]=-COM(0,1)*cur[1];
	- cur[3]=cur[0];
	- } else { // if backward
	- double sqpim=sqrt(pin.minus());
	- cur[0]=-sqpomsqpimpoperp/abs(poperp);
	- cur[1]=-sqpim*sqpop;
	- cur[2]=COM(0,1)*cur[1];
	- cur[3]=-cur[0];
	- }
	- }
	- CCurrent temp(cur[0],cur[1],cur[2],cur[3]);
	- return temp;
	+ return joi(pout, helout, pin, helin);
	}


	// Current for <incoming state \| mu \| outgoing state>
	void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur) {

	cur[0] = 0.0;
	cur[1] = 0.0;
	cur[2] = 0.0;
	cur[3] = 0.0;
	- if(helin!=helout){
	- std::cout<<__LINE__<<" "<<__FILE__<<std::endl;
	- }
	- double sqpop = sqrt(pout.plus());
	- double sqpom = sqrt(pout.minus());
	- COM poperp = pout.x() + COM(0, 1) * pout.y();
	+ const double sqpop = sqrt(pout.plus());
	+ const double sqpom = sqrt(pout.minus());
	+ const COM poperp = pout.x() + COM(0, 1) * pout.y();

	- if (helout == false) {

	+ if (helout != helin) {
	+ throw std::invalid_argument{"Non-matching helicities"};
	+ } else if (helout == false) { // negative helicity
	if (pin.plus() > pin.minus()) { // if forward
	- double sqpip=sqrt(pin.plus());
	+ const double sqpip = sqrt(pin.plus());
	cur[0] = sqpop * sqpip;
	cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp);
	cur[2] = COM(0,1) * cur[1];
	cur[3] = cur[0];
	- }
	- else {
	- double sqpim = sqrt(pin.minus());
	+ } else { // if backward
	+ const double sqpim = sqrt(pin.minus());
	cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp);
	cur[1] = -sqpim * sqpop;
	cur[2] = -COM(0,1) * cur[1];
	cur[3] = -cur[0];
	}
	- }
	-
	- else {
	+ } else { // positive helicity
	if (pin.plus() > pin.minus()) { // if forward
	- double sqpip = sqrt(pin.plus());
	+ const double sqpip = sqrt(pin.plus());
	cur[0] = sqpop * sqpip;
	- cur[1] = sqpom * sqpip*poperp/abs(poperp);
	- cur[2] = -COM(0,1)*cur[1];
	+ cur[1] = sqpom * sqpip * poperp / abs(poperp);
	+ cur[2] = -COM(0,1) * cur[1];
	cur[3] = cur[0];
	- }
	- else {
	- double sqpim = sqrt(pin.minus());
	- cur[0] = -sqpom * sqpim * poperp/abs(poperp);
	+ } else { // if backward
	+ const double sqpim = sqrt(pin.minus());
	+ cur[0] = -sqpom * sqpim * poperp / abs(poperp);
	cur[1] = -sqpim * sqpop;
	- cur[2] = COM(0,1)*cur[1];
	+ cur[2] = COM(0,1) * cur[1];
	cur[3] = -cur[0];
	}
	}
	}

	+CCurrent jio (HLV pin, bool helin, HLV pout, bool helout)
	+{
	+ current cur;
	+ jio(pin, helin, pout, helout, cur);
	+ return CCurrent(cur[0],cur[1],cur[2],cur[3]);
	+}
	+

	// Current for <outgoing state \| mu \| outgoing state>
	void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) {

	// Zero our current
	cur[0] = 0.0;
	cur[1] = 0.0;
	cur[2] = 0.0;
	cur[3] = 0.0;
	- if(helj){
	- std::cout<<__LINE__<<" "<<__FILE__<<std::endl;
	- }
	- // If positive helicity swap momenta
	- if (heli == true) {
	- HLV dummy;
	- dummy = pi;
	- pi = pj;
	- pj = dummy;
	+ if (heli!=helj) {
	+ throw std::invalid_argument{"Non-matching helicities"};
	+ } else if ( heli == true ) { // If positive helicity swap momenta
	+ std::swap(pi,pj);
	}

	- double sqpjp = sqrt(pj.plus());
	- double sqpjm = sqrt(pj.minus());
	- double sqpip = sqrt(pi.plus());
	- double sqpim = sqrt(pi.minus());
	+ const double sqpjp = sqrt(pj.plus());
	+ const double sqpjm = sqrt(pj.minus());
	+ const double sqpip = sqrt(pi.plus());
	+ const double sqpim = sqrt(pi.minus());

	- COM piperp = pi.x() + COM(0,1) * pi.y();
	- COM pjperp = pj.x() + COM(0,1) * pj.y();
	- COM phasei = piperp / abs(piperp);
	- COM phasej = pjperp / abs(pjperp);
	+ const COM piperp = pi.x() + COM(0,1) * pi.y();
	+ const COM pjperp = pj.x() + COM(0,1) * pj.y();
	+ const COM phasei = piperp / abs(piperp);
	+ const COM phasej = pjperp / abs(pjperp);

	cur[0] = sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
	cur[1] = sqpim * sqpjp * phasei + sqpip * sqpjm * conj(phasej);
	cur[2] = -COM(0, 1) * (sqpim * sqpjp * phasei - sqpip * sqpjm * conj(phasej));
	cur[3] = -sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
	}

	-CCurrent joo (CLHEP::HepLorentzVector pi, bool heli, CLHEP::HepLorentzVector pj, bool helj)
	+CCurrent joo (HLV pi, bool heli, HLV pj, bool helj)
	{
	- COM cur[4];
	-
	- if (heli!=helj) {
	- std::cerr<< "void j : Non-matching helicities\n";
	- } else if (heli==true) { // negative helicity
	- CLHEP::HepLorentzVector dummy;
	- dummy=pi;
	- pi=pj;
	- pj=dummy;
	- }
	- double sqpjp=sqrt(pj.plus());
	- double sqpjm=sqrt(pj.minus());
	- double sqpip=sqrt(pi.plus());
	- double sqpim=sqrt(pi.minus());
	- COM piperp=pi.x()+COM(0,1)*pi.y();
	- COM pjperp=pj.x()+COM(0,1)*pj.y();
	- COM phasei=piperp/abs(piperp);
	- COM phasej=pjperp/abs(pjperp);
	-
	- cur[0]=sqpimsqpjmphaseiconj(phasej)+sqpipsqpjp;
	- cur[1]=sqpimsqpjpphasei+sqpipsqpjmconj(phasej);
	- cur[2]=-COM(0,1)(sqpimsqpjpphasei-sqpipsqpjm*conj(phasej));
	- cur[3]=-sqpimsqpjmphaseiconj(phasej)+sqpipsqpjp;
	-
	- CCurrent temp(cur[0],cur[1],cur[2],cur[3]);
	- return temp;
	-}
	-
	-
	-// Current Functions
	-// Current for <outgoing state \| mu \| incoming state>
	-void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) {
	-
	- cur[0] = 0.0;
	- cur[1] = 0.0;
	- cur[2] = 0.0;
	- cur[3] = 0.0;
	- if(helin){
	- std::cout<<__LINE__<<" "<<__FILE__<<std::endl;
	- }
	- double sqpop = sqrt(pout.plus());
	- double sqpom = sqrt(pout.minus());
	- COM poperp = pout.x() + COM(0, 1) * pout.y();
	-
	- if (helout == false) {
	-
	- if (pin.plus() > pin.minus()) { // if forward
	- double sqpip=sqrt(pin.plus());
	- cur[0] = sqpop * sqpip;
	- cur[1] = sqpom * sqpip * poperp/abs(poperp);
	- cur[2] = -COM(0,1)*cur[1];
	- cur[3] = cur[0];
	- }
	- else {
	- double sqpim = sqrt(pin.minus());
	- cur[0] = -sqpomsqpimpoperp/abs(poperp);
	- cur[1] = -sqpim*sqpop;
	- cur[2] = COM(0,1)*cur[1];
	- cur[3] = -cur[0];
	- }
	- }
	- else {
	- if (pin.plus() > pin.minus()) { // if forward
	- double sqpip = sqrt(pin.plus());
	- cur[0] = sqpop * sqpip;
	- cur[1] = sqpom * sqpip*conj(poperp)/abs(poperp);
	- cur[2] = COM(0,1)*cur[1];
	- cur[3] = cur[0];
	- }
	- else {
	- double sqpim = sqrt(pin.minus());
	- cur[0] = -sqpom * sqpim * conj(poperp)/abs(poperp);
	- cur[1] = -sqpim * sqpop;
	- cur[2] = -COM(0,1)*cur[1];
	- cur[3] = -cur[0];
	- }
	- }
	+ current cur;
	+ joo(pi, heli, pj, helj, cur);
	+ return CCurrent(cur[0],cur[1],cur[2],cur[3]);
	}

	namespace {
	/// @TODO unused function
	// double jM2 (CLHEP::HepLorentzVector p1out, bool hel1out, CLHEP::HepLorentzVector p1in, bool hel1in, CLHEP::HepLorentzVector p2out, bool hel2out, CLHEP::HepLorentzVector p2in, bool hel2in)
	// {
	// CLHEP::HepLorentzVector q1=p1in-p1out;
	// CLHEP::HepLorentzVector q2=-(p2in-p2out);
	// current C1,C2;
	// j (p1out,hel1out,p1in,hel1in, C1);
	// j (p2out,hel2out,p2in,hel2in, C2);

	// std::cout << "# From Currents, C1 : ("<<C1[0]<<","<<C1[1]<<","<<C1[2]<<","<<C1[3]<<"\n";
	// std::cout << "# From Currents, C2 : ("<<C2[0]<<","<<C2[1]<<","<<C2[2]<<","<<C2[3]<<"\n";

	// COM M=cdot(C1,C2);

	// return (Mconj(M)).real()/(q1.m2()q2.m2());
	// }

	void jW (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin, current cur)
	{
	// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5)
	// Need to swap e and nu for events with W- --> e- nubar!
	if (helin==helout && hele==helnu) {
	CLHEP::HepLorentzVector qa=pout+pe+pnu;
	CLHEP::HepLorentzVector qb=pin-pe-pnu;
	double ta(qa.m2()),tb(qb.m2());

	current t65,vout,vin,temp2,temp3,temp5;
	joo(pnu,helnu,pe,hele,t65);
	vout[0]=pout.e();
	vout[1]=pout.x();
	vout[2]=pout.y();
	vout[3]=pout.z();
	vin[0]=pin.e();
	vin[1]=pin.x();
	vin[2]=pin.y();
	vin[3]=pin.z();

	COM brac615=cdot(t65,vout);
	COM brac645=cdot(t65,vin);

	// prod1565 and prod6465 are zero for Ws (not Zs)!!
	// noalias(temp)=prod(trans(CurrentOutOut(pout,helout,pnu,helout)),metric);
	joo(pout,helout,pnu,helout,temp2);
	// noalias(temp2)=prod(temp,ctemp);
	COM prod1665=cdot(temp2,t65);
	// noalias(temp)=prod(trans(Current(pe,helin,pin,helin)),metric);
	// noalias(temp2)=prod(temp,ctemp);
	j(pe,helin,pin,helin,temp3);
	COM prod5465=cdot(temp3,t65);
	// noalias(temp)=prod(trans(Current(pnu,helin,pin,helin)),metric);
	// noalias(temp2)=prod(temp,ctemp);

	joo(pout,helout,pe,helout,temp2);
	j(pnu,helnu,pin,helin,temp3);
	j(pout,helout,pin,helin,temp5);

	current term1,term2,term3,sum;
	cmult(2.brac615/ta+2.brac645/tb,temp5,term1);
	cmult(prod1665/ta,temp3,term2);
	cmult(-prod5465/tb,temp2,term3);

	// cur=((2.brac615Current(pout,helout,pin,helin)+prod1565Current(pe,helin,pin,helin)+prod1665Current(pnu,helin,pin,helin))/ta + (2.brac645Current(pout,helout,pin,helin)-prod5465CurrentOutOut(pout,helout,pe,helout)-prod6465CurrentOutOut(pout,helout,pnu,helout))/tb);
	// cur=((2.brac615temp5+prod1565temp3+prod1665temp4)/ta + (2.brac645temp5-prod5465temp1-prod6465temp2)/tb);
	cadd(term1,term2,term3,sum);
	// std::cout<<"sum: ("<<sum[0]<<","<<sum[1]<<","<<sum[2]<<","<<sum[3]<<")\n";
	cur[0]=sum[0];
	cur[1]=sum[1];
	cur[2]=sum[2];
	cur[3]=sum[3];
	}
	}

	void jWbar (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin, current cur)
	{
	// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5)
	// Need to swap e and nu for events with W- --> e- nubar!
	if (helin==helout && hele==helnu) {
	CLHEP::HepLorentzVector qa=pout+pe+pnu;
	CLHEP::HepLorentzVector qb=pin-pe-pnu;
	double ta(qa.m2()),tb(qb.m2());

	current t65,vout,vin,temp2,temp3,temp5;
	joo(pnu,helnu,pe,hele,t65);
	vout[0]=pout.e();
	vout[1]=pout.x();
	vout[2]=pout.y();
	vout[3]=pout.z();
	vin[0]=pin.e();
	vin[1]=pin.x();
	vin[2]=pin.y();
	vin[3]=pin.z();

	COM brac615=cdot(t65,vout);
	COM brac645=cdot(t65,vin);

	// prod1565 and prod6465 are zero for Ws (not Zs)!!
	joo(pe,helout,pout,helout,temp2); // temp2 is <5\|alpha\|1>
	COM prod5165=cdot(temp2,t65);
	jio(pin,helin,pnu,helin,temp3); // temp3 is <4\|alpha\|6>
	COM prod4665=cdot(temp3,t65);

	joo(pnu,helout,pout,helout,temp2); // temp2 is now <6\|mu\|1>
	jio(pin,helin,pe,helin,temp3); // temp3 is now <4\|mu\|5>
	jio(pin,helin,pout,helout,temp5); // temp5 is <4\|mu\|1>

	current term1,term2,term3,sum;
	cmult(-2.brac615/ta-2.brac645/tb,temp5,term1);
	cmult(-prod5165/ta,temp3,term2);
	cmult(prod4665/tb,temp2,term3);

	// cur=((2.brac615Current(pout,helout,pin,helin)+prod1565Current(pe,helin,pin,helin)+prod1665Current(pnu,helin,pin,helin))/ta + (2.brac645Current(pout,helout,pin,helin)-prod5465CurrentOutOut(pout,helout,pe,helout)-prod6465CurrentOutOut(pout,helout,pnu,helout))/tb);
	// cur=((2.brac615temp5+prod1565temp3+prod1665temp4)/ta + (2.brac645temp5-prod5465temp1-prod6465temp2)/tb);
	cadd(term1,term2,term3,sum);
	// std::cout<<"term1: ("<<temp5[0]<<" "<<temp5[1]<<" "<<temp5[2]<<" "<<temp5[3]<<")"<<std::endl;
	// std::cout<<"sum: ("<<sum[0]<<","<<sum[1]<<","<<sum[2]<<","<<sum[3]<<")\n";
	cur[0]=sum[0];
	cur[1]=sum[1];
	cur[2]=sum[2];
	cur[3]=sum[3];
	}
	}
	} // namespace anonymous

	double jMWqQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	// Calculates the square of the current contractions for qQ->qenuQ scattering
	// p1: quark (with W emittance)
	// p2: Quark
	{
	current mj1m,mj2p,mj2m;
	CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);

	jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m);
	j(p2out,true,p2in,true,mj2p);
	j(p2out,false,p2in,false,mj2m);

	// std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n";
	// std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n";
	// std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n";

	// mj1m.mj2p

	COM Mmp=cdot(mj1m,mj2p);

	// mj1m.mj2m
	COM Mmm=cdot(mj1m,mj2m);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);

	// // Leave division by colour and Helicity avg until Tree files
	// Leave multi. of couplings to later
	// Multiply by Cf^2
	return (4./3.)(4./3.)(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());

	}

	double jMWqQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	// Calculates the square of the current contractions for qQ->qenuQ scattering
	// p1: quark (with W emittance)
	// p2: Quark
	{
	current mj1m,mj2p,mj2m;
	CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);

	jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m);
	jio(p2in,true,p2out,true,mj2p);
	jio(p2in,false,p2out,false,mj2m);

	// std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n";
	// std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n";
	// std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n";

	// mj1m.mj2p

	COM Mmp=cdot(mj1m,mj2p);

	// mj1m.mj2m
	COM Mmm=cdot(mj1m,mj2m);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);

	// // Leave division by colour and Helicity avg until Tree files
	// Leave multi. of couplings to later
	// Multiply by Cf^2
	return (4./3.)(4./3.)(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());

	}

	double jMWqbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	// Calculates the square of the current contractions for qQ->qenuQ scattering
	// p1: quark (with W emittance)
	// p2: Quark
	{
	current mj1m,mj2p,mj2m;
	CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);

	jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m);
	j(p2out,true,p2in,true,mj2p);
	j(p2out,false,p2in,false,mj2m);

	// std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n";
	// std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n";
	// std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n";

	// mj1m.mj2p

	COM Mmp=cdot(mj1m,mj2p);

	// mj1m.mj2m
	COM Mmm=cdot(mj1m,mj2m);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);

	// // Leave division by colour and Helicity avg until Tree files
	// Leave multi. of couplings to later
	// Multiply by Cf^2
	return (4./3.)(4./3.)(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());

	}

	double jMWqbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	// Calculates the square of the current contractions for qQ->qenuQ scattering
	// p1: quark (with W emittance)
	// p2: Quark
	{
	current mj1m,mj2p,mj2m;
	CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);

	jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m);
	jio(p2in,true,p2out,true,mj2p);
	jio(p2in,false,p2out,false,mj2m);

	// std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n";
	// std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n";
	// std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n";

	// mj1m.mj2p

	COM Mmp=cdot(mj1m,mj2p);

	// mj1m.mj2m
	COM Mmm=cdot(mj1m,mj2m);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);

	// // Leave division by colour and Helicity avg until Tree files
	// Leave multi. of couplings to later
	// Multiply by Cf^2
	return (4./3.)(4./3.)(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());

	}

	double jMWqg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	// Calculates the square of the current contractions for qg->qenug scattering
	// p1: quark
	// p2: gluon
	{
	CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);
	current mj1m,mj2p,mj2m;

	jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m);

	j(p2out,true,p2in,true,mj2p);
	j(p2out,false,p2in,false,mj2m);

	// mj1m.mj2p

	COM Mmp=cdot(mj1m,mj2p);

	// mj1m.mj2m
	COM Mmm=cdot(mj1m,mj2m);

	const double K = K_g(p2out, p2in);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);
	double sst = K/C_A*(a2Mmp+a2Mmm);
	// double sstsave=sst;

	// // Leave division by colour and Helicity avg until Tree files
	// Leave multi. of couplings to later
	// Multiply by Cf*Ca=4
	return 4.sst/(q1.m2()q2.m2());

	}

	double jMWqbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	// Calculates the square of the current contractions for qg->qenug scattering
	// p1: quark
	// p2: gluon
	{
	CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);
	current mj1m,mj2p,mj2m;

	jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m);

	j(p2out,true,p2in,true,mj2p);
	j(p2out,false,p2in,false,mj2m);

	// mj1m.mj2p

	COM Mmp=cdot(mj1m,mj2p);

	// mj1m.mj2m
	COM Mmm=cdot(mj1m,mj2m);

	const double K = K_g(p2out, p2in);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);
	double sst = K/C_A*(a2Mmp+a2Mmm);
	// double sstsave=sst;

	// // Leave division by colour and Helicity avg until Tree files
	// Leave multi. of couplings to later
	// Multiply by Cf*Ca=4
	return 4.sst/(q1.m2()q2.m2());

	}

	double jM2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	{
	// std::cerr<<"Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<std::endl;

	CLHEP::HepLorentzVector q1=p1in-p1out;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);
	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	current mj1m,mj1p,mj2m,mj2p;
	j(p1out,true,p1in,true,mj1p);
	j(p1out,false,p1in,false,mj1m);
	j(p2out,true,p2in,true,mj2p);
	j(p2out,false,p2in,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	double sst=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);

	// Multiply by Cf^2
	return RHEJ::C_FRHEJ::C_F(sst)/(q1.m2()*q2.m2());
	}

	double jM2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	{
	CLHEP::HepLorentzVector q1=p1in-p1out;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);
	current mj1m,mj1p,mj2m,mj2p;
	j(p1out,true,p1in,true,mj1p);
	j(p1out,false,p1in,false,mj1m);
	jio(p2in,true,p2out,true,mj2p);
	jio(p2in,false,p2out,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	double sumsq=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);

	// Multiply by Cf^2
	return (4./3.)(4./3.)(sumsq)/(q1.m2()*q2.m2());
	}

	double jM2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	{
	CLHEP::HepLorentzVector q1=p1in-p1out;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);
	current mj1m,mj1p,mj2m,mj2p;
	jio(p1in,true,p1out,true,mj1p);
	jio(p1in,false,p1out,false,mj1m);
	jio(p2in,true,p2out,true,mj2p);
	jio(p2in,false,p2out,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	double sumsq=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);

	// Multiply by Cf^2
	return (4./3.)(4./3.)(sumsq)/(q1.m2()*q2.m2());
	}

	double jM2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	// Calculates the square of the current contractions for qg scattering
	// p1: quark
	// p2: gluon
	{
	CLHEP::HepLorentzVector q1=p1in-p1out;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);

	current mj1m,mj1p,mj2m,mj2p;
	j(p1out,true,p1in,true,mj1p);
	j(p1out,false,p1in,false,mj1m);
	j(p2out,true,p2in,true,mj2p);
	j(p2out,false,p2in,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	const double K = K_g(p2out, p2in);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);
	double a2Mpp=abs2(Mpp);
	double a2Mpm=abs2(Mpm);
	double sst = K/C_A*(a2Mpp+a2Mpm+a2Mmp+a2Mmm);
	// double sstsave=sst;

	// std::cout <<"ratio: "<<sst/sstsave<<std::endl;
	// Cf*Ca=4
	return 4.sst/(q1.m2()q2.m2());

	}

	double jM2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	// Calculates the square of the current contractions for qg scattering
	// p1: quark
	// p2: gluon
	{
	CLHEP::HepLorentzVector q1=p1in-p1out;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);

	current mj1m,mj1p,mj2m,mj2p;
	jio(p1in,true,p1out,true,mj1p);
	jio(p1in,false,p1out,false,mj1m);
	j(p2out,true,p2in,true,mj2p);
	j(p2out,false,p2in,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	const double K = K_g(p2out, p2in);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);
	double a2Mpp=abs2(Mpp);
	double a2Mpm=abs2(Mpm);
	double sst = K/C_A*(a2Mpp+a2Mpm+a2Mmp+a2Mmm);
	// double sstsave=sst;

	// std::cout <<"ratio: "<<sst/sstsave<<std::endl;
	// Cf*Ca=4
	return 4.sst/(q1.m2()q2.m2());

	}

	double jM2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
	// Calculates the square of the current contractions for gg scattering
	// p1: gluon
	// p2: gluon
	{
	CLHEP::HepLorentzVector q1=p1in-p1out;
	CLHEP::HepLorentzVector q2=-(p2in-p2out);

	current mj1m,mj1p,mj2m,mj2p;
	j(p1out,true,p1in,true,mj1p);
	j(p1out,false,p1in,false,mj1m);
	j(p2out,true,p2in,true,mj2p);
	j(p2out,false,p2in,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	const double K_g1 = K_g(p1out, p1in);
	const double K_g2 = K_g(p2out, p2in);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);
	double a2Mpp=abs2(Mpp);
	double a2Mpm=abs2(Mpm);
	double sst = K_g1/C_AK_g2/C_A(a2Mpp+a2Mpm+a2Mmp+a2Mmm);
	// double sstsave=sst;
	// std::cout <<"ratio: "<<sst/sstsave<<std::endl;
	// Ca*Ca=9
	return 9.sst/(q1.m2()q2.m2());

	}

	namespace {
	/// @TODO what was this intended to do?
	// double MH2helper(current C1, current C2, current q1, current q2)
	// {
	// COM M;
	// COM temp1,temp2;
	// // First the C1.q2 * C2.q1 - part
	// temp1=cdot(C1,q2);
	// temp2=cdot(C2,q1);
	// M=temp1*temp2;

	// // Then the C1.C2 * q1.q2
	// temp1=cdot(C1,C2);
	// temp2=cdot(q1,q2);
	// M-=temp1*temp2;

	// return (M*conj(M)).real();
	// }

	/**
	* @brief Higgs vertex contracted with current @param C1 and @param C2
	*/
	COM cHdot(const current & C1, const current & C2, const current & q1,
	const current & q2, double mt, bool incBot, double mb)
	{
	if (mt == infinity) {
	return (cdot(C1,C2)cdot(q1,q2)-cdot(C1,q2)cdot(C2,q1))/(6M_PIv);
	}
	else {
	CLHEP::HepLorentzVector vq1,vq2;
	vq1.set(q1[1].real(),q1[2].real(),q1[3].real(),q1[0].real());
	vq2.set(q2[1].real(),q2[2].real(),q2[3].real(),q2[0].real());
	// first minus sign obtained because of q1-difference to VDD
	// std::cout<<"A1 : " << A1(-vq1,vq2)<<std::endl;
	// std::cout<<"A2 : " << A2(-vq1,vq2)<<std::endl;
	if(!(incBot))
	// Factor is because 4 mt^2 g^2/v A1 -> 16 pi mt^2/v alphas,
	// and we divide by a factor 4 at the amp sqaured level later
	// which I absorb here (i.e. I divide by 2)
	/// @TODO move factor 1/2 from S to \|ME\|^2 => consistent with general notation
	return 8.M_PImtmt/v(-cdot(C1,q2)cdot(C2,q1)A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt));
	else
	return 8.M_PImtmt/v(-cdot(C1,q2)cdot(C2,q1)A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt))
	+ 8.M_PImbmb/v(-cdot(C1,q2)cdot(C2,q1)A1(-vq1,vq2,mb)-cdot(C1,C2)*A2(-vq1,vq2,mb));
	}
	}
	} // namespace anonymous

	double MH2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2,
	double mt, bool incBot, double mb)
	{
	// CLHEP::HepLorentzVector q1=p1in-p1out;
	// CLHEP::HepLorentzVector q2=-(p2in-p2out);
	current j1p,j1m,j2p,j2m, q1v, q2v;

	j (p1out,true,p1in,true,j1p);
	j (p1out,false,p1in,false,j1m);

	j (p2out,true,p2in,true,j2p);
	j (p2out,false,p2in,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);
	COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);

	double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm);
	// return (4./3.)(4./3.)sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	return sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
	{
	// CLHEP::HepLorentzVector q1=p1in-p1out;
	// CLHEP::HepLorentzVector q2=-(p2in-p2out);
	current j1p,j1m,j2p,j2m,q1v,q2v;

	j (p1out,true,p1in,true,j1p);
	j (p1out,false,p1in,false,j1m);

	jio (p2in,true,p2out,true,j2p);
	jio (p2in,false,p2out,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);
	COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);

	double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm);
	// return (4./3.)(4./3.)sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	return sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2qbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
	{
	// CLHEP::HepLorentzVector q1=p1in-p1out;
	// CLHEP::HepLorentzVector q2=-(p2in-p2out);
	current j1p,j1m,j2p,j2m,q1v,q2v;

	jio (p1in,true,p1out,true,j1p);
	jio (p1in,false,p1out,false,j1m);

	j (p2out,true,p2in,true,j2p);
	j (p2out,false,p2in,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);
	COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);

	double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm);
	// return (4./3.)(4./3.)sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	return sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
	{
	// CLHEP::HepLorentzVector q1=p1in-p1out;
	// CLHEP::HepLorentzVector q2=-(p2in-p2out);
	current j1p,j1m,j2p,j2m,q1v,q2v;

	jio (p1in,true,p1out,true,j1p);
	jio (p1in,false,p1out,false,j1m);

	jio (p2in,true,p2out,true,j2p);
	jio (p2in,false,p2out,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);
	COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);

	double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm);
	// return (4./3.)(4./3.)sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	return sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
	// q~p1 g~p2 (i.e. ALWAYS p1 for quark, p2 for gluon)
	// should be called with q1 meant to be contracted with p2 in first part of vertex
	// (i.e. if g is backward, q1 is forward)
	{
	current j1p,j1m,j2p,j2m,q1v,q2v;

	j (p1out,true,p1in,true,j1p);
	j (p1out,false,p1in,false,j1m);

	j (p2out,true,p2in,true,j2p);
	j (p2out,false,p2in,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	// First, calculate the non-flipping amplitudes:

	COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);
	COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);

	//cout << "Bits in MH2qg: " << Mpp << " " << Mpm << " " << Mmp << " " << Mmm << endl;

	const double K = K_g(p2out, p2in);

	double sst=K/C_A*(abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm));

	// Cf*Ca=4
	// return 4.sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()q2.m2());
	return sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
	// qbar~p1 g~p2 (i.e. ALWAYS p1 for anti-quark, p2 for gluon)
	// should be called with q1 meant to be contracted with p2 in first part of vertex
	// (i.e. if g is backward, q1 is forward)
	{
	current j1p,j1m,j2p,j2m,q1v,q2v;

	jio (p1in,true,p1out,true,j1p);
	jio (p1in,false,p1out,false,j1m);

	j (p2out,true,p2in,true,j2p);
	j (p2out,false,p2in,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	// First, calculate the non-flipping amplitudes:

	COM amp,amm,apm,app;
	app=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	apm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);
	amp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	amm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);

	double MH2sum = abs2(app)+abs2(amm)+abs2(apm)+abs2(amp);

	const double K = K_g(p2out, p2in);
	MH2sum*=K/C_A;

	// Cf*Ca=4
	// return 4.MH2sum/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()q2.m2());
	return MH2sum/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
	// g~p1 g~p2
	// should be called with q1 meant to be contracted with p2 in first part of vertex
	// (i.e. if g is backward, q1 is forward)
	{
	current j1p,j1m,j2p,j2m,q1v,q2v;

	j (p1out,true,p1in,true,j1p);
	j (p1out,false,p1in,false,j1m);

	j (p2out,true,p2in,true,j2p);
	j (p2out,false,p2in,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	// First, calculate the non-flipping amplitudes:

	COM amp,amm,apm,app;
	app=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	apm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);
	amp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	amm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);

	double MH2sum = abs2(app)+abs2(amm)+abs2(apm)+abs2(amp);

	const double K_g1 = K_g(p1out, p1in);
	const double K_g2 = K_g(p2out, p2in);

	MH2sum=K_g1/C_AK_g2/C_A;

	// Ca*Ca=9
	// return 9.MH2sum/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()q2.m2());
	return MH2sum/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	// // Z's stuff
	// void jZ(HLV pin, HLV pout, HLV pem, HLV pep, bool HelPartons, bool HelLeptons, current cur) {

	// // Init current to zero
	// cur[0] = 0.0;
	// cur[1] = 0.0;
	// cur[2] = 0.0;
	// cur[3] = 0.0;

	// // Temporary variables
	// COM temp;
	// current Term_1, Term_2, Term_3, Term_4, J_temp, TempCur1, TempCur2;

	// // Momentum of virtual gluons aroun weak boson emission site
	// HLV qa = pout + pep + pem;
	// HLV qb = pin - pep - pem;

	// double ta = qa.m2();
	// double tb = qb.m2();

	// // Out-Out currents:
	// current Em_Ep, Out_Em, Out_Ep;

	// // Other currents:
	// current Out_In, Em_In, Ep_In;

	// joi(pout, HelPartons, pin, HelPartons, Out_In);
	// joi(pem, HelLeptons, pin, HelPartons, Em_In);
	// joi(pep, HelLeptons, pin, HelPartons, Ep_In);

	// joo(pem, HelLeptons, pep, HelLeptons, Em_Ep);
	// joo(pout, HelPartons, pem, HelLeptons, Out_Em);
	// joo(pout, HelPartons, pep, HelLeptons, Out_Ep);

	// if (HelLeptons == HelPartons) {

	// temp = 2.0 * cdot(pout, Em_Ep);
	// cmult(temp / ta, Out_In, Term_1);

	// temp = cdot(Out_Em, Em_Ep);
	// cmult(temp / ta , Em_In, Term_2);

	// temp = 2.0 * cdot(pin, Em_Ep);
	// cmult(temp / tb, Out_In, Term_3);

	// temp = -cdot(Ep_In, Em_Ep);
	// cmult(temp / tb, Out_Ep, Term_4);

	// cadd(Term_1, Term_2, Term_3, Term_4, J_temp);

	// cur[0] = J_temp[0];
	// cur[1] = J_temp[1];
	// cur[2] = J_temp[2];
	// cur[3] = J_temp[3];
	// }

	// else {
	// if (HelPartons == true) {
	// temp = 2.0 * cdot(pout, Em_Ep);
	// cmult(temp / ta, Out_In, Term_1);

	// joo(pout, true, pep, true, TempCur1);
	// joi(pep, true, pin, true, TempCur2);

	// temp = cdot(TempCur1, Em_Ep);
	// cmult(temp / ta , TempCur2, Term_2);

	// temp = 2.0 * cdot(pin, Em_Ep);
	// cmult(temp / tb, Out_In, Term_3);

	// joo(pout, true, pem, true, TempCur1);
	// joi(pem, true, pin, true, TempCur2);

	// temp = -cdot(TempCur2, Em_Ep);
	// cmult(temp / tb, TempCur1, Term_4);

	// cadd(Term_1, Term_2, Term_3, Term_4, J_temp);

	// cur[0] = J_temp[0];
	// cur[1] = J_temp[1];
	// cur[2] = J_temp[2];
	// cur[3] = J_temp[3];
	// }

	// else {
	// temp = 2.0 * cdot(pout, Em_Ep);
	// cmult(temp / ta, Out_In, Term_1);

	// joo(pout, false, pep, false, TempCur1);
	// joi(pep, false, pin, false, TempCur2);

	// temp = cdot(TempCur1, Em_Ep);
	// cmult(temp / ta, TempCur2, Term_2);

	// temp = 2.0 * cdot(pin, Em_Ep);
	// cmult(temp / tb, Out_In, Term_3);

	// joo(pout, false, pem, false, TempCur1);
	// joi(pem, false, pin, false, TempCur2);

	// temp = -cdot(TempCur2, Em_Ep);
	// cmult(temp / tb, TempCur1, Term_4);

	// cadd(Term_1, Term_2, Term_3, Term_4, J_temp);

	// cur[0] = J_temp[0];
	// cur[1] = J_temp[1];
	// cur[2] = J_temp[2];
	// cur[3] = J_temp[3];
	// }

	// }
	// }

	// void jZbar(HLV pin, HLV pout, HLV pem, HLV pep, bool HelPartons, bool HelLeptons, current cur) {

	// // Init current to zero
	// cur[0] = 0.0;
	// cur[1] = 0.0;
	// cur[2] = 0.0;
	// cur[3] = 0.0;

	// // Temporary variables
	// COM temp;
	// current Term_1, Term_2, Term_3, Term_4, J_temp, TempCur1, TempCur2;

	// // Transfered 4-momenta
	// HLV qa = pout + pep + pem;
	// HLV qb = pin - pep - pem;

	// // The square of the transfered 4-momenta
	// double ta = qa.m2();
	// double tb = qb.m2();

	// // Out-Out currents:
	// current Em_Ep, Em_Out, Ep_Out;

	// // In-Out currents:
	// current In_Out, In_Em, In_Ep;

	// // Safe to use the currents since helicity structure is ok
	// if (HelPartons == HelLeptons) {
	// jio(pin, HelPartons, pout, HelPartons, In_Out);
	// joo(pem, HelLeptons, pep, HelLeptons, Em_Ep);
	// jio(pin, HelPartons, pem, HelLeptons, In_Em);
	// jio(pin, HelPartons, pep, HelLeptons, In_Ep);
	// joo(pem, HelLeptons, pout, HelPartons, Em_Out);
	// joo(pep, HelLeptons, pout, HelPartons, Ep_Out);
	// }

	// else {
	// jio(pin, HelPartons, pout, HelPartons, In_Out);
	// joo(pem, HelLeptons, pep, HelLeptons, Em_Ep);

	// In_Em[0] = 0.0;
	// In_Em[1] = 0.0;
	// In_Em[2] = 0.0;
	// In_Em[3] = 0.0;

	// In_Ep[0] = 0.0;
	// In_Ep[1] = 0.0;
	// In_Ep[2] = 0.0;
	// In_Ep[3] = 0.0;

	// Em_Out[0] = 0.0;
	// Em_Out[1] = 0.0;
	// Em_Out[2] = 0.0;
	// Em_Out[3] = 0.0;

	// Ep_Out[0] = 0.0;
	// Ep_Out[1] = 0.0;
	// Ep_Out[2] = 0.0;
	// Ep_Out[3] = 0.0;
	// }

	// if (HelLeptons == HelPartons) {

	// temp = 2.0 * cdot(pout, Em_Ep);
	// cmult(temp / ta, In_Out, Term_1);

	// temp = cdot(Ep_Out, Em_Ep);
	// cmult(temp / ta, In_Ep, Term_2);

	// temp = 2.0 * cdot(pin, Em_Ep);
	// cmult(temp / tb, In_Out, Term_3);

	// temp = - cdot(In_Em, Em_Ep);
	// cmult(temp / tb, Em_Out, Term_4);

	// cadd(Term_1, Term_2, Term_3, Term_4, J_temp);

	// cur[0] = J_temp[0];
	// cur[1] = J_temp[1];
	// cur[2] = J_temp[2];
	// cur[3] = J_temp[3];
	// }

	// else {
	// if (HelPartons == true) {

	// temp = 2.0 * cdot(pout, Em_Ep);
	// cmult(temp / ta, In_Out, Term_1);

	// joo(pem, true, pout, true, TempCur1);
	// jio(pin, true, pem, true, TempCur2);

	// temp = cdot(TempCur1, Em_Ep);
	// cmult(temp / ta , TempCur2, Term_2);

	// temp = 2.0 * cdot(pin, Em_Ep);
	// cmult(temp / tb, In_Out, Term_3);

	// joo(pep, true, pout, true, TempCur1);
	// jio(pin, true, pep, true, TempCur2);

	// temp = - cdot(TempCur2, Em_Ep);
	// cmult(temp / tb, TempCur1, Term_4);

	// cadd(Term_1, Term_2, Term_3, Term_4, J_temp);

	// cur[0] = J_temp[0];
	// cur[1] = J_temp[1];
	// cur[2] = J_temp[2];
	// cur[3] = J_temp[3];
	// }

	// else {

	// temp = 2.0 * cdot(pout, Em_Ep);
	// cmult(temp / ta, In_Out, Term_1);

	// joo(pem, false, pout, false, TempCur1);
	// jio(pin, false, pem, false, TempCur2);

	// temp = cdot(TempCur1, Em_Ep);
	// cmult(temp / ta , TempCur2, Term_2);

	// temp = 2.0 * cdot(pin, Em_Ep);
	// cmult(temp / tb, In_Out, Term_3);

	// joo(pep, false, pout, false, TempCur1);
	// jio(pin, false, pep, false, TempCur2);

	// temp = - cdot(TempCur2, Em_Ep);
	// cmult(temp / tb, TempCur1, Term_4);

	// cadd(Term_1, Term_2, Term_3, Term_4, J_temp);

	// cur[0] = J_temp[0];
	// cur[1] = J_temp[1];
	// cur[2] = J_temp[2];
	// cur[3] = J_temp[3];
	// }
	// }
	// }

	// // Progagators
	// COM PZ(double s) {

	// double MZ, GammaZ;

	// MZ = 9.118800e+01; // Mass of the mediating gauge boson
	// GammaZ = 2.441404e+00; // Z peak width

	// // Return Z Prop value
	// return 1.0 / (s - MZ * MZ + COM(0.0, 1.0) * GammaZ * MZ);
	// }

	// COM PG(double s) {
	// return 1.0 / s;
	// }

	// // Non-gluonic with pa emitting
	// std::vector <double> jMZqQ (HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem, std::vector <double> VProducts, std::vector < std::vector <double> > Virtuals, int aptype, int bptype, bool UseVirtuals, bool BottomLineEmit) {

	// std::vector <double> ScaledWeights;

	// double Sum;

	// // Propagator factors
	// COM PZs = PZ((pep + pem).m2());
	// COM PGs = PG((pep + pem).m2());


	// // Emitting current initialisation
	// current j1pptop, j1pmtop; // Emission from top line
	// current j1ppbot, j1pmbot; // Emission from bottom line

	// // Non-emitting current initialisation
	// current j2ptop, j2mtop; // Emission from top line
	// current j2pbot, j2mbot; // Emission from bottom line

	// // Currents for top emission
	// // Upper current calculations
	// // if a is a quark
	// if (aptype > 0) {
	// jZ(pa, p1, pem, pep, true, true, j1pptop);
	// jZ(pa, p1, pem, pep, true, false, j1pmtop);
	// }
	// // if a is an antiquark
	// else {
	// jZbar(pa, p1, pem, pep, true, true, j1pptop);
	// jZbar(pa, p1, pem, pep, true, false, j1pmtop);
	// }

	// // Lower current calculations
	// // if b is a quark
	// if (bptype > 0) {
	// joi(p2, true, pb, true, j2ptop);
	// joi(p2, false, pb, false, j2mtop);
	// }
	// // if b is an antiquark
	// else {
	// jio(pb, true, p2, true, j2ptop);
	// jio(pb, false, p2, false, j2mtop);
	// }

	// // Currents for bottom emission
	// // Lower current calculations
	// if (bptype > 0) {
	// jZ(pb, p2, pem, pep, true, true, j1ppbot);
	// jZ(pb, p2, pem, pep, true, false, j1pmbot);
	// }
	// else {
	// jZbar(pb, p2, pem, pep, true, true, j1ppbot);
	// jZbar(pb, p2, pem, pep, true, false, j1pmbot);
	// }

	// // Upper current calculations
	// if (aptype > 0) {
	// joi(p1, true, pa, true, j2pbot);
	// joi(p1, false, pa, false, j2mbot);
	// }
	// else {
	// jio(pa, true, p1, true, j2pbot);
	// jio(pa, false, p1, false, j2mbot);
	// }

	// COM Coeff[2][8];

	// if (!Interference) {

	// double ZCharge_a_P = Zq(aptype, true);
	// double ZCharge_a_M = Zq(aptype, false);
	// double ZCharge_b_P = Zq(bptype, true);
	// double ZCharge_b_M = Zq(bptype, false);

	// if (BottomLineEmit) {
	// // Emission from top-line quark (pa/p1 line)
	// Coeff[0][0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2ptop);
	// Coeff[0][1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2mtop);
	// Coeff[0][2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2ptop);
	// Coeff[0][3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2mtop);
	// Coeff[0][4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop));
	// Coeff[0][5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop));
	// Coeff[0][6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop));
	// Coeff[0][7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop));
	// }

	// else {
	// // Emission from bottom-line quark (pb/p2 line)
	// Coeff[1][0] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2pbot);
	// Coeff[1][7] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2mbot);
	// Coeff[1][2] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2pbot);
	// Coeff[1][5] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2mbot);
	// Coeff[1][4] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2pbot));
	// Coeff[1][3] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2mbot));
	// Coeff[1][6] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2pbot));
	// Coeff[1][1] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2mbot));
	// }
	// }

	// // Else calculate all the possiblities
	// else {

	// double ZCharge_a_P = Zq(aptype, true);
	// double ZCharge_a_M = Zq(aptype, false);
	// double ZCharge_b_P = Zq(bptype, true);
	// double ZCharge_b_M = Zq(bptype, false);

	// // Emission from top-line quark (pa/p1 line)
	// Coeff[0][0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2ptop);
	// Coeff[0][1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2mtop);
	// Coeff[0][2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2ptop);
	// Coeff[0][3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2mtop);
	// Coeff[0][4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop));
	// Coeff[0][5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop));
	// Coeff[0][6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop));
	// Coeff[0][7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop));

	// // Emission from bottom-line quark (pb/p2 line)
	// Coeff[1][0] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2pbot);
	// Coeff[1][7] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2mbot);
	// Coeff[1][2] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2pbot);
	// Coeff[1][5] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2mbot);
	// Coeff[1][4] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2pbot));
	// Coeff[1][3] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2mbot));
	// Coeff[1][6] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2pbot));
	// Coeff[1][1] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2mbot));
	// }

	// // Find the numbers of scales
	// int ScaleCount;
	// #if calcscaleunc
	// ScaleCount = 20;
	// #else
	// ScaleCount = 1;
	// #endif

	// // For each scale...
	// for (int j = 0; j < ScaleCount; j++) {

	// Sum = 0.0;

	// // If we want to compare back to the W's code only emit from one quark and only couple to left handed particles
	// // virtuals arent here since they are calculated and included in weight() call.
	// if (!Interference) {

	// if (BottomLineEmit) for (int i = 0; i < 8; i++) Sum += abs2(Coeff[1][i]) * VProducts.at(1);
	// else for (int i = 0; i < 8; i++) Sum += abs2(Coeff[0][i]) * VProducts.at(0);
	// }

	// // Else work out the full interference
	// else {

	// // For the full calculation...
	// if (UseVirtuals) {
	// for (int i = 0; i < 8; i++) {
	// Sum += abs2(Coeff[0][i]) * VProducts.at(0) * Virtuals.at(j).at(0)
	// + abs2(Coeff[1][i]) * VProducts.at(1) * Virtuals.at(j).at(1)
	// + 2.0 * real(Coeff[0][i] * conj(Coeff[1][i])) * VProducts.at(2) * Virtuals.at(j).at(2);
	// }
	// }

	// // For the tree level calculation...
	// else {
	// for (int i = 0; i < 8; i++) {
	// Sum += abs2(Coeff[0][i]) * VProducts.at(0)
	// + abs2(Coeff[1][i]) * VProducts.at(1)
	// + 2.0 * real(Coeff[0][i] * conj(Coeff[1][i])) * VProducts.at(2);
	// }
	// }
	// }

	// // Add this to the vector to be returned with the other factors of C_A and the helicity sum/average factors.
	// ScaledWeights.push_back(Sum / 18.0);
	// }

	// // Return all the scale values
	// return ScaledWeights;
	// }

	// // Semi-gluonic with pa emitting
	// std::vector <double> jMZqg (HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem, std::vector <double> VProducts, std::vector < std::vector <double> > Virtuals, int aptype, int bptype, bool UseVirtuals, bool BottomLineEmit) {

	// COM Coeff[8];

	// double Sum;

	// std::vector <double> ScaledWeights;

	// COM PZs = PZ((pep + pem).m2());
	// COM PGs = PG((pep + pem).m2());

	// // Emitting current initialisation - Emission from top line
	// current j1pptop, j1pmtop;

	// // Non-emitting current initialisation - Emission from top line
	// current j2ptop, j2mtop;

	// // Currents for top emission
	// // Upper current calculations
	// if (aptype > 0) {
	// jZ (pa, p1, pem, pep, true, true, j1pptop);
	// jZ (pa, p1, pem, pep, true, false, j1pmtop);
	// }
	// else {
	// jZbar(pa, p1, pem, pep, true, true, j1pptop);
	// jZbar(pa, p1, pem, pep, true, false, j1pmtop);
	// }

	// // Lower current calculations
	// joi(p2, true, pb, true, j2ptop);
	// joi(p2, false, pb, false, j2mtop);

	// // Calculate all the possiblities
	// double ZCharge_a_P = Zq(aptype, true);
	// double ZCharge_a_M = Zq(aptype, false);

	// // Emission from top-line quark (pa/p1 line)
	// Coeff[0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2ptop);
	// Coeff[1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2mtop);
	// Coeff[2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2ptop);
	// Coeff[3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2mtop);
	// Coeff[4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop));
	// Coeff[5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop));
	// Coeff[6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop));
	// Coeff[7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop));

	// // Calculate gluon colour accelerated factor
	// double CAMFactor, z;

	// // If b is a forward moving gluon define z (C.F. multiple jets papers)
	// if (pb.pz() > 0) z = p2.plus() / pb.plus();
	// else z = p2.minus() / pb.minus();

	// CAMFactor = (1.0 - 1.0 / 9.0) / 2.0 * (z + 1.0 / z) + 1.0 / 9.0;

	// // Find the numbers of scales
	// int ScaleCount;
	// #if calcscaleunc
	// ScaleCount = 20;
	// #else
	// ScaleCount = 1;
	// #endif

	// // For each scale...
	// for (int j = 0; j < ScaleCount; j++) {

	// Sum = 0.0;

	// // If we dont want the interference
	// if (!Interference) for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0);

	// // Else work out the full interference
	// else {
	// if (UseVirtuals) {
	// for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0) * Virtuals.at(j).at(0);
	// }
	// else {
	// for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0);
	// }
	// }

	// // Add this to the vector to be returned with the other factors of C_A, the colour accelerated factor and the helicity sum/average factors.: (4/3)*3/32
	// ScaledWeights.push_back(CAMFactor * Sum / 8.0);
	// }

	// return ScaledWeights;
	// }

	// // Electroweak Charge Functions
	// double Zq (int PID, bool Helcitiy) {

	// double temp;

	// // Positive Spin
	// if (Helcitiy == true) {
	// if (PID == 1 \|\| PID == 3 \|\| PID == 5) temp = (+ 1.0 * stw2 / 3.0) / ctw;
	// if (PID == 2 \|\| PID == 4) temp = (- 2.0 * stw2 / 3.0) / ctw;
	// if (PID == -1 \|\| PID == -3 \|\| PID == -5) temp = (- 1.0 * stw2 / 3.0) / ctw;
	// if (PID == -2 \|\| PID == -4) temp = (+ 2.0 * stw2 / 3.0) / ctw;

	// // If electron or positron
	// if (PID == 7 \|\| PID == -7) temp = Zep;
	// }

	// // Negative Spin
	// else {
	// if (PID == 1 \|\| PID == 3 \|\| PID == 5) temp = (-0.5 + 1.0 * stw2 / 3.0) / ctw;
	// if (PID == 2 \|\| PID == 4) temp = ( 0.5 - 2.0 * stw2 / 3.0) / ctw;
	// if (PID == -1 \|\| PID == -3 \|\| PID == -5) temp = ( 0.5 - 1.0 * stw2 / 3.0) / ctw;
	// if (PID == -2 \|\| PID == -4) temp = (-0.5 + 2.0 * stw2 / 3.0) / ctw;

	// // If electron or positron
	// if (PID == 7 \|\| PID == -7) temp = Zem;

	// }

	// return temp;
	// }

	// double Gq (int PID) {

	// if (!VirtualPhoton) return 0.0;

	// if (PID == -1) return 1.0 * ee / 3.0;
	// if (PID == -2) return -2.0 * ee / 3.0;
	// if (PID == -3) return 1.0 * ee / 3.0;
	// if (PID == -4) return -2.0 * ee / 3.0;
	// if (PID == -5) return 1.0 * ee / 3.0;
	// if (PID == 1) return -1.0 * ee / 3.0;
	// if (PID == 2) return 2.0 * ee / 3.0;
	// if (PID == 3) return -1.0 * ee / 3.0;
	// if (PID == 4) return 2.0 * ee / 3.0;
	// if (PID == 5) return -1.0 * ee / 3.0;

	// std::cout << "ERROR! No Electroweak Charge Found at line " << __LINE__ << "..." << std::endl;
	// return 0.0;
	// }


	namespace {

	//@{
	/// @brief Higgs vertex contracted with one current

	CCurrent jH (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin,
	bool helin, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2,
	double mt, bool incBot, double mb)
	{

	CCurrent j2 = j(pout,helout,pin,helin);
	CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz());

	if(mt == infinity)
	return ((q1.dot(q2))j2 - j2.dot(q1)jq2)/(3M_PIv);
	else
	{
	if(incBot)
	return (-16.M_PImbmb/vj2.dot(q1)jq2A1(-q1,q2,mb)-16.M_PImbmb/vj2*A2(-q1,q2,mb))
	+ (-16.M_PImtmt/vj2.dot(q1)jq2A1(-q1,q2,mt)-16.M_PImtmt/vj2*A2(-q1,q2,mt));
	else
	return (-16.M_PImtmt/vj2.dot(q1)jq2A1(-q1,q2,mt)-16.M_PImtmt/vj2*A2(-q1,q2,mt));
	}
	}

	CCurrent jioH (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout,
	bool helout, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2,
	double mt, bool incBot, double mb)
	{

	CCurrent j2 = jio(pin,helin,pout,helout);
	CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz());

	if(mt == infinity)
	return ((q1.dot(q2))j2 - j2.dot(q1)jq2)/(3M_PIv);
	else
	{
	if(incBot)
	return (-16.M_PImbmb/vj2.dot(q1)jq2A1(-q1,q2,mb)-16.M_PImbmb/vj2*A2(-q1,q2,mb))
	+ (-16.M_PImtmt/vj2.dot(q1)jq2A1(-q1,q2,mt)-16.M_PImtmt/vj2*A2(-q1,q2,mt));
	else
	return (-16.M_PImtmt/vj2.dot(q1)jq2A1(-q1,q2,mt)-16.M_PImtmt/vj2*A2(-q1,q2,mt));
	}
	}

	CCurrent jHtop (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin,
	bool helin, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2,
	double mt, bool incBot, double mb)
	{

	CCurrent j1 = j(pout,helout,pin,helin);
	CCurrent jq1(q1.e(),q1.px(),q1.py(),q1.pz());

	if(mt == infinity)
	return ((q1.dot(q2))j1 - j1.dot(q2)jq1)/(3M_PIv);
	else
	{
	if(incBot)
	return (-16.M_PImbmb/vj1.dot(q2)jq1A1(-q1,q2,mb)-16.M_PImbmb/vj1*A2(-q1,q2,mb))
	+ (-16.M_PImtmt/vj1.dot(q2)jq1A1(-q1,q2,mt)-16.M_PImtmt/vj1*A2(-q1,q2,mt));
	else
	return (-16.M_PImtmt/vj1.dot(q2)jq1A1(-q1,q2,mt)-16.M_PImtmt/vj1*A2(-q1,q2,mt));
	}
	}


	CCurrent jioHtop (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout,
	bool helout, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2,
	double mt, bool incBot, double mb)
	{

	CCurrent j1 = jio(pin,helin,pout,helout);
	CCurrent jq1(q1.e(),q1.px(),q1.py(),q1.pz());

	if(mt == infinity)
	return ((q1.dot(q2))j1 - j1.dot(q2)jq1)/(3M_PIv);
	else
	{
	if(incBot)
	return (-16.M_PImbmb/vj1.dot(q2)jq1A1(-q1,q2,mb)-16.M_PImbmb/vj1*A2(-q1,q2,mb))
	+ (-16.M_PImtmt/vj1.dot(q2)jq1A1(-q1,q2,mt)-16.M_PImtmt/vj1*A2(-q1,q2,mt));
	else
	return (-16.M_PImtmt/vj1.dot(q2)jq1A1(-q1,q2,mt)-16.M_PImtmt/vj1*A2(-q1,q2,mt));
	}
	}
	//@}
	} // namespace anonymous

	double jM2unogqHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{
	// This construction is taking rapidity order: pg > p1out >> p2out
	// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
	CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon

	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=j(p1out,true,p1in,true);
	mj1m=j(p1out,false,p1in,false);
	mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// std::cout<< p1out.rapidity() << " " << p2out.rapidity()<< " " << qH1 << " " << qH2 << "\n" <<MHmm << " " << MHmp << " " << MHpm << " " << MHpp << std::endl;

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(p1out,true,pg,true);
	j2gm=joo(p1out,false,pg,false);
	jgap=j(pg,true,p1in,true);
	jgam=j(pg,false,p1in,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	Lmm=(qsum(MHmm) + (-2.mjH2m.dot(pg))mj1m+2.mj1m.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmm/2.))/q1.m2();
	Lmp=(qsum(MHmp) + (-2.mjH2p.dot(pg))mj1m+2.mj1m.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmp/2.))/q1.m2();
	Lpm=(qsum(MHpm) + (-2.mjH2m.dot(pg))mj1p+2.mj1p.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpm/2.))/q1.m2();
	Lpp=(qsum(MHpp) + (-2.mjH2p.dot(pg))mj1p+2.mj1p.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpp/2.))/q1.m2();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=RHEJ::C_FRHEJ::C_F/RHEJ::C_A/RHEJ::C_A; // Factor of (Cf/Ca) for each quark to match MH2qQ.
	//Higgs coupling is included in Hjets.C

	return ampsq;
	}

	double jM2unogqbarHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{
	// This construction is taking rapidity order: pg > p1out >> p2out
	// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
	CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon

	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=jio(p1in,true,p1out,true);
	mj1m=jio(p1in,false,p1out,false);
	mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(pg,true,p1out,true);
	j2gm=joo(pg,false,p1out,false);
	jgap=jio(p1in,true,pg,true);
	jgam=jio(p1in,false,pg,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	Lmm=(qsum(MHmm) + (-2.mjH2m.dot(pg))mj1m+2.mj1m.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmm/2.))/q1.m2();
	Lmp=(qsum(MHmp) + (-2.mjH2p.dot(pg))mj1m+2.mj1m.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmp/2.))/q1.m2();
	Lpm=(qsum(MHpm) + (-2.mjH2m.dot(pg))mj1p+2.mj1p.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpm/2.))/q1.m2();
	Lpp=(qsum(MHpp) + (-2.mjH2p.dot(pg))mj1p+2.mj1p.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpp/2.))/q1.m2();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
	//Higgs coupling is included in Hjets.C

	return ampsq;
	}

	double jM2unogqHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{
	// This construction is taking rapidity order: pg > p1out >> p2out
	// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
	CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon

	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=j(p1out,true,p1in,true);
	mj1m=j(p1out,false,p1in,false);
	mjH2p=jioH(p2in,true,p2out,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jioH(p2in,false,p2out,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(p1out,true,pg,true);
	j2gm=joo(p1out,false,pg,false);
	jgap=j(pg,true,p1in,true);
	jgam=j(pg,false,p1in,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	Lmm=(qsum(MHmm) + (-2.mjH2m.dot(pg))mj1m+2.mj1m.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmm/2.))/q1.m2();
	Lmp=(qsum(MHmp) + (-2.mjH2p.dot(pg))mj1m+2.mj1m.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmp/2.))/q1.m2();
	Lpm=(qsum(MHpm) + (-2.mjH2m.dot(pg))mj1p+2.mj1p.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpm/2.))/q1.m2();
	Lpp=(qsum(MHpp) + (-2.mjH2p.dot(pg))mj1p+2.mj1p.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpp/2.))/q1.m2();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
	//Higgs coupling is included in Hjets.C

	return ampsq;
	}

	double jM2unogqbarHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{
	// This construction is taking rapidity order: pg > p1out >> p2out
	// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
	CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon

	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=jio(p1in,true,p1out,true);
	mj1m=jio(p1in,false,p1out,false);
	mjH2p=jioH(p2in,true,p2out,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jioH(p2in,false,p2out,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(pg,true,p1out,true);
	j2gm=joo(pg,false,p1out,false);
	jgap=jio(p1in,true,pg,true);
	jgam=jio(p1in,false,pg,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	Lmm=(qsum(MHmm) + (-2.mjH2m.dot(pg))mj1m+2.mj1m.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmm/2.))/q1.m2();
	Lmp=(qsum(MHmp) + (-2.mjH2p.dot(pg))mj1m+2.mj1m.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmp/2.))/q1.m2();
	Lpm=(qsum(MHpm) + (-2.mjH2m.dot(pg))mj1p+2.mj1p.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpm/2.))/q1.m2();
	Lpp=(qsum(MHpp) + (-2.mjH2p.dot(pg))mj1p+2.mj1p.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpp/2.))/q1.m2();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	//Higgs coupling is included in Hjets.C
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.

	return ampsq;
	}

	double jM2unogqHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{
	// This construction is taking rapidity order: pg > p1out >> p2out
	// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
	CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon

	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=j(p1out,true,p1in,true);
	mj1m=j(p1out,false,p1in,false);
	mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(p1out,true,pg,true);
	j2gm=joo(p1out,false,pg,false);
	jgap=j(pg,true,p1in,true);
	jgam=j(pg,false,p1in,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	Lmm=(qsum(MHmm) + (-2.mjH2m.dot(pg))mj1m+2.mj1m.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmm/2.))/q1.m2();
	Lmp=(qsum(MHmp) + (-2.mjH2p.dot(pg))mj1m+2.mj1m.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmp/2.))/q1.m2();
	Lpm=(qsum(MHpm) + (-2.mjH2m.dot(pg))mj1p+2.mj1p.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpm/2.))/q1.m2();
	Lpp=(qsum(MHpp) + (-2.mjH2p.dot(pg))mj1p+2.mj1p.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpp/2.))/q1.m2();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4./9.4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ.
	// here we need 2 to match with the normalization
	// gq is 9./4. times the qQ
	//Higgs coupling is included in Hjets.C

	const double K = K_g(p2out, p2in);

	return ampsqK/C_A9./4.; //ca/cf = 9/4
	}

	double jM2unogqbarHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{
	// This construction is taking rapidity order: pg > p1out >> p2out
	// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
	CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon

	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=jio(p1in,true,p1out,true);
	mj1m=jio(p1in,false,p1out,false);
	mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(pg,true,p1out,true);
	j2gm=joo(pg,false,p1out,false);
	jgap=jio(p1in,true,pg,true);
	jgam=jio(p1in,false,pg,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	Lmm=(qsum(MHmm) + (-2.mjH2m.dot(pg))mj1m+2.mj1m.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmm/2.))/q1.m2();
	Lmp=(qsum(MHmp) + (-2.mjH2p.dot(pg))mj1m+2.mj1m.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHmp/2.))/q1.m2();
	Lpm=(qsum(MHpm) + (-2.mjH2m.dot(pg))mj1p+2.mj1p.dot(pg)mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpm/2.))/q1.m2();
	Lpp=(qsum(MHpp) + (-2.mjH2p.dot(pg))mj1p+2.mj1p.dot(pg)mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()*MHpp/2.))/q1.m2();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4./9.4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ.
	// here we need 2 to match with the normalization
	// gq is 9./4. times the qQ
	//Higgs coupling is included in Hjets.C

	const double K = K_g(p2out, p2in);

	return ampsq*K/C_F;
	}

	double jM2unobqHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{
	// std::cout << "####################\n";
	// std::cout << "# p1in : "<<p1in<< " "<<p1in.plus()<<" "<<p1in.minus()<<std::endl;
	// std::cout << "# p2in : "<<p2in<< " "<<p2in.plus()<<" "<<p2in.minus()<<std::endl;
	// std::cout << "# p1out : "<<p1out<< " "<<p1out.rapidity()<<std::endl;
	// std::cout << "# (qH1-qH2) : "<<(qH1-qH2)<< " "<<(qH1-qH2).rapidity()<<std::endl;
	// std::cout << "# pg : "<<pg<< " "<<pg.rapidity()<<std::endl;
	// std::cout << "# p2out : "<<p2out<< " "<<p2out.rapidity()<<std::endl;
	// std::cout << "####################\n";

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End

	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2, mt, incBot, mb);
	mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2, mt, incBot, mb);
	mj2p=j(p2out,true,p2in,true);
	mj2m=j(p2out,false,p2in,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);

	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(p2out,true,pg,true);
	j2gm=joo(p2out,false,pg,false);
	jgbp=j(pg,true,p2in,true);
	jgbm=j(pg,false,p2in,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	// COM test=pminus.dot(p1in);

	Lmm=((-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m+2.mj2m.dot(pg)mjH1m
	+ (p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmm/2.))/q3.m2();
	Lmp=((-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p+2.mj2p.dot(pg)mjH1m
	+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmp/2.))/q3.m2();
	Lpm=((-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m+2.mj2m.dot(pg)mjH1p
	+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpm/2.))/q3.m2();
	Lpp=((-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p+2.mj2p.dot(pg)mjH1p
	+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpp/2.))/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	// 1/3. = 1/C_A ?
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	const double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=RHEJ::C_FRHEJ::C_F/(RHEJ::C_A*RHEJ::C_A); // Factor of (Cf/Ca) for each quark to match MH2qQ.

	return ampsq;
	}



	double jM2unobqbarHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End

	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jioHtop(p1in,true,p1out,true,qH1,qH2, mt, incBot, mb);
	mjH1m=jioHtop(p1in,false,p1out,false,qH1,qH2, mt, incBot, mb);
	mj2p=j(p2out,true,p2in,true);
	mj2m=j(p2out,false,p2in,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);


	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(p2out,true,pg,true);
	j2gm=joo(p2out,false,pg,false);
	jgbp=j(pg,true,p2in,true);
	jgbm=j(pg,false,p2in,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	// COM test=pminus.dot(p1in);


	Lmm=((-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m+2.mj2m.dot(pg)mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmm/2.))/q3.m2();
	Lmp=((-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p+2.mj2p.dot(pg)mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmp/2.))/q3.m2();
	Lpm=((-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m+2.mj2m.dot(pg)mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpm/2.))/q3.m2();
	Lpp=((-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p+2.mj2p.dot(pg)mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpp/2.))/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
	//Higgs coupling is included in Hjets.C

	return ampsq;
	}



	double jM2unobqHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End

	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2,mt, incBot, mb);
	mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2,mt, incBot, mb);
	mj2p=jio(p2in,true,p2out,true);
	mj2m=jio(p2in,false,p2out,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);


	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(pg,true,p2out,true);
	j2gm=joo(pg,false,p2out,false);
	jgbp=jio(p2in,true,pg,true);
	jgbm=jio(p2in,false,pg,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	// COM test=pminus.dot(p1in);


	Lmm=((-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m+2.mj2m.dot(pg)mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmm/2.))/q3.m2();
	Lmp=((-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p+2.mj2p.dot(pg)mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmp/2.))/q3.m2();
	Lpm=((-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m+2.mj2m.dot(pg)mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpm/2.))/q3.m2();
	Lpp=((-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p+2.mj2p.dot(pg)mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpp/2.))/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
	//Higgs coupling is included in Hjets.C

	return ampsq;
	}



	double jM2unobqbarHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End

	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jioHtop(p1in,true,p1out,true,qH1,qH2,mt, incBot, mb);
	mjH1m=jioHtop(p1in,false,p1out,false,qH1,qH2,mt, incBot, mb);
	mj2p=jio(p2in,true,p2out,true);
	mj2m=jio(p2in,false,p2out,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);


	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(pg,true,p2out,true);
	j2gm=joo(pg,false,p2out,false);
	jgbp=jio(p2in,true,pg,true);
	jgbm=jio(p2in,false,pg,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	// COM test=pminus.dot(p1in);


	Lmm=((-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m+2.mj2m.dot(pg)mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmm/2.))/q3.m2();
	Lmp=((-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p+2.mj2p.dot(pg)mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmp/2.))/q3.m2();
	Lpm=((-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m+2.mj2m.dot(pg)mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpm/2.))/q3.m2();
	Lpp=((-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p+2.mj2p.dot(pg)mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpp/2.))/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
	//Higgs coupling is included in Hjets.C



	return ampsq;
	}

	double jM2unobgHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{
	// std::cout << "####################\n";
	// std::cout << "# p1in : "<<p1in<< " "<<p1in.plus()<<" "<<p1in.minus()<<std::endl;
	// std::cout << "# p2in : "<<p2in<< " "<<p2in.plus()<<" "<<p2in.minus()<<std::endl;
	// std::cout << "# p1out : "<<p1out<< " "<<p1out.rapidity()<<std::endl;
	// std::cout << "# (qH1-qH2) : "<<(qH1-qH2)<< " "<<(qH1-qH2).rapidity()<<std::endl;
	// std::cout << "# pg : "<<pg<< " "<<pg.rapidity()<<std::endl;
	// std::cout << "# p2out : "<<p2out<< " "<<p2out.rapidity()<<std::endl;
	// std::cout << "####################\n";

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End

	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2,mt, incBot, mb);
	mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2,mt, incBot, mb);
	mj2p=j(p2out,true,p2in,true);
	mj2m=j(p2out,false,p2in,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);

	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(p2out,true,pg,true);
	j2gm=joo(p2out,false,pg,false);
	jgbp=j(pg,true,p2in,true);
	jgbm=j(pg,false,p2in,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	// COM test=pminus.dot(p1in);

	Lmm=((-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m+2.mj2m.dot(pg)mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmm/2.))/q3.m2();
	Lmp=((-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p+2.mj2p.dot(pg)mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmp/2.))/q3.m2();
	Lpm=((-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m+2.mj2m.dot(pg)mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpm/2.))/q3.m2();
	Lpp=((-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p+2.mj2p.dot(pg)mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpp/2.))/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4./9.4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ.
	// need twice to match the normalization
	//Higgs coupling is included in Hjets.C

	const double K = K_g(p1out, p1in);

	return ampsq*K/C_F;
	}



	double jM2unobgHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
	{

	CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
	CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End

	// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2,mt, incBot, mb);
	mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2,mt, incBot, mb);
	mj2p=jio(p2in,true,p2out,true);
	mj2m=jio(p2in,false,p2out,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);


	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(pg,true,p2out,true);
	j2gm=joo(pg,false,p2out,false);
	jgbp=jio(p2in,true,pg,true);
	jgbm=jio(p2in,false,pg,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	// COM test=pminus.dot(p1in);


	Lmm=((-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m+2.mj2m.dot(pg)mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmm/2.))/q3.m2();
	Lmp=((-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p+2.mj2p.dot(pg)mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmp/2.))/q3.m2();
	Lpm=((-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m+2.mj2m.dot(pg)mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpm/2.))/q3.m2();
	Lpp=((-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p+2.mj2p.dot(pg)mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpp/2.))/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=RHEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
	// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
	// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;

	// Now add the t-channels for the Higgs
	double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4./9.4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ.
	//Higgs coupling is included in Hjets.C

	const double K = K_g(p1out, p1in);

	return ampsq*K/C_F; //ca/cf = 9/4

	}

	// Begin finite mass stuff
	#ifdef RHEJ_BUILD_WITH_QCDLOOP
	namespace {


	// All the stuff needed for the box functions in qg->qgH now...

	//COM E1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM E1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//CLHEP::HepLorentzVector q2=k3+k4;
	CLHEP::HepLorentzVector q2=-(k1+k2+kh);
	double Delta, Sigma, S1, S2, s12, s34;
	S1 = 2.*k1.dot(q2);
	S2 = 2.*k2.dot(q2);
	s12 = 2.*k1.dot(k2);
	//s34 = 2.*k3.dot(k4);
	s34 = q2.m2();
	Delta = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor(-s12D0DD(k2, k1, q2, mq)(1 - 8.mqmq/s12 + S2/(2.s12) +
	S2(s12 - 8.mqmq)(s34 + S1)/(2.s12Delta) +
	2.(s34 + S1)(s34 + S1)/Delta +
	S2pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)C0DD(k2,
	k1 + q2, mq) -
	s12C0DD(k1, k2, mq) + (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S1*C0DD(k1, q2,
	mq))(S2(s12 - 4.mqmq)/(2.s12Delta) +
	2.*(s34 + S1)/Delta +
	S2*pow((s34 + S1),2)/Delta/Delta) + (C0DD(k1, q2, mq) -
	C0DD(k1 + k2, q2, mq))(1. - 4.mq*mq/s12) -
	C0DD(k1 + k2, q2, mq)2.s34/
	S1 - (B0DD(k1 + q2, mq) -
	B0DD(k1 + k2 + q2, mq))2.s34*(s34 +
	S1)/(S1*Delta) + (B0DD(q2, mq) -
	B0DD(k1 + k2 + q2, mq) +
	s12*C0DD(k1 + k2, q2,
	mq))(2.s34*(s34 +
	S1)(S1 - S2)/(DeltaSigma) +
	2.s34(s34 + S1)/(S1*Delta)) + (B0DD(k1 + k2, mq) -
	B0DD(k1 + k2 + q2,
	mq) - (s34 + S1 + S2)C0DD(k1 + k2, q2, mq))2.*(s34 +
	S1)(2.s12*s34 -
	S2(S1 + S2))/(DeltaSigma));
	}
	//COM F1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM F1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//CLHEP::HepLorentzVector q2=k3+k4;
	CLHEP::HepLorentzVector q2 = -(k1+k2+kh);
	double Delta, Sigma, S1, S2, s12, s34;
	S1 = 2.*k1.dot(q2);
	S2 = 2.*k2.dot(q2);
	s12 = 2.*k1.dot(k2);
	//s34 = 2.*k3.dot(k4);
	s34 = q2.m2();
	Delta = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor(-S2D0DD(k1, k2, q2,
	mq)(0.5 - (s12 - 8.mqmq)(s34 + S2)/(2.*Delta) -
	s12pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)C0DD(k1,
	k2 + q2, mq) -
	s12C0DD(k1, k2, mq) - (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S2*C0DD(k2, q2,
	mq))(S2(s12 - 4.mqmq)/(2.s12Delta) +
	S2pow((s34 + S2),2)/Delta/Delta) - (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))(1. - 4.mqmq/s12) -
	C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) -
	B0DD(k1 + k2 + q2,
	mq))2.pow((s34 + S2),2)/((s12 + S1)*Delta) - (B0DD(
	q2, mq) - B0DD(k1 + k2 + q2, mq) +
	s12C0DD(k1 + k2, q2, mq))2.s34(s34 +
	S2)(S2 - S1)/(DeltaSigma) + (B0DD(
	k1 + k2, mq) -
	B0DD(k1 + k2 + q2,
	mq) - (s34 + S1 + S2)C0DD(k1 + k2, q2, mq))2.*(s34 +
	S2)(2.s12*s34 -
	S2(S1 + S2))/(DeltaSigma));
	}
	//COM G1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM G1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//CLHEP::HepLorentzVector q2=k3+k4;
	CLHEP::HepLorentzVector q2 = -(k1+k2+kh);
	double Delta, S1, S2, s12, s34;
	S1 = 2.*k1.dot(q2);
	S2 = 2.*k2.dot(q2);
	s12 = 2.*k1.dot(k2);
	//s34 = 2.*k3.dot(k4);
	s34 = q2.m2();
	Delta = s12s34 - S1S2;

	return looprwfactor(S2D0DD(k1, q2, k2,
	mq)(Delta/s12/s12 - 4.mq*mq/s12) -
	S2((s12 + S1)C0DD(k1, k2 + q2, mq) -
	S1C0DD(k1, q2, mq))(1./
	s12/s12 - (s12 - 4.mqmq)/(2.s12Delta)) -
	S2((s12 + S2)C0DD(k1 + q2, k2, mq) -
	S2C0DD(k2, q2, mq))(1./
	s12/s12 + (s12 - 4.mqmq)/(2.s12Delta)) -
	C0DD(k1, q2, mq) - (C0DD(k1, k2 + q2, mq) -
	C0DD(k1, q2, mq))4.mq*mq/
	s12 + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2./
	s12 + (B0DD(k1 + q2, mq) -
	B0DD(q2, mq))2.s34/(s12*S1) + (B0DD(k2 + q2, mq) -
	B0DD(k1 + k2 + q2, mq))2.(s34 + S2)/(s12*(s12 + S1)));
	}
	//COM E4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM E4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//CLHEP::HepLorentzVector q2=k3+k4;
	CLHEP::HepLorentzVector q2 = -(k1+k2+kh);
	double Delta, Sigma, S1, S2, s12, s34;
	S1 = 2.*k1.dot(q2);
	S2 = 2.*k2.dot(q2);
	s12 = 2.*k1.dot(k2);
	//s34 = 2.*k3.dot(k4);
	s34 = q2.m2();
	Delta = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor* (-s12*D0DD(k2, k1, q2,
	mq)(0.5 - (S1 - 8.mqmq)(s34 + S1)/(2.*Delta) -
	s12pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)C0DD(k2,
	k1 + q2, mq) -
	s12C0DD(k1, k2, mq) + (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S1C0DD(k1, q2, mq))((S1 - 4.mqmq)/(2.*Delta) +
	s12*pow((s34 + S1),2)/Delta/Delta) -
	C0DD(k1 + k2, q2, mq) + (B0DD(k1 + q2, mq) -
	B0DD(k1 + k2 + q2, mq))(2.s34/Delta +
	2.s12(s34 + S1)/((s12 + S2)*Delta)) - (B0DD(
	q2, mq) - B0DD(k1 + k2 + q2, mq) +
	s12*C0DD(k1 + k2, q2,
	mq))((2.s34(2.s12s34 - S2(S1 + S2) +
	s12*(S1 -
	S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) -
	B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)C0DD(k1 + k2, q2, mq))((2.s12(2.s12s34 - S1*(S1 + S2) +
	s34(S2 - S1)))/(DeltaSigma)));
	}
	//COM F4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM F4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//CLHEP::HepLorentzVector q2=k3+k4;
	CLHEP::HepLorentzVector q2 = -(k1+k2+kh);
	double Delta, Sigma, S1, S2, s12, s34;
	S1 = 2.*k1.dot(q2);
	S2 = 2.*k2.dot(q2);
	s12 = 2.*k1.dot(k2);
	//s34 = 2.*k3.dot(k4);
	s34 = q2.m2();
	Delta = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor* (-s12*D0DD(k1, k2, q2,
	mq)(0.5 + (S1 - 8.mqmq)(s34 + S2)/(2.*Delta) +
	s12pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)C0DD(k1,
	k2 + q2, mq) -
	s12C0DD(k1, k2, mq) - (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S2C0DD(k2, q2, mq))((S1 - 4.mqmq)/(2.*Delta) +
	s12*pow((s34 + S2),2)/Delta/Delta) -
	C0DD(k1 + k2, q2, mq) - (B0DD(k2 + q2, mq) -
	B0DD(k1 + k2 + q2, mq))2.(s34 +
	S2)/Delta + (B0DD(q2, mq) -
	B0DD(k1 + k2 + q2, mq) +
	s12C0DD(k1 + k2, q2, mq))2.s34(2.s12s34 -
	S1*(S1 + S2) +
	s12(S2 - S1))/(DeltaSigma) - (B0DD(k1 + k2, mq) -
	B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)C0DD(k1 + k2, q2, mq))(2.s12(2.s12s34 - S2*(S1 + S2) +
	s34(S1 - S2))/(DeltaSigma)));
	}
	//COM G4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM G4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//CLHEP::HepLorentzVector q2=k3+k4;
	CLHEP::HepLorentzVector q2 = -(k1+k2+kh);
	double Delta, S1, S2, s12, s34;
	S1 = 2.*k1.dot(q2);
	S2 = 2.*k2.dot(q2);
	s12 = 2.*k1.dot(k2);
	//s34 = 2.*k3.dot(k4);
	s34 = q2.m2();
	Delta = s12s34 - S1S2;

	return looprwfactor* (-D0DD(k1, q2, k2,
	mq)*(Delta/s12 + (s12 + S1)/2. -
	4.mqmq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
	S1C0DD(k1, q2, mq))(1./
	s12 - (S1 - 4.mqmq)/(2.Delta)) + ((s12 + S2)C0DD(
	k1 + q2, k2, mq) -
	S2C0DD(k2, q2, mq))(1./
	s12 + (S1 - 4.mqmq)/(2.*Delta)) + (B0DD(
	k1 + k2 + q2, mq) -
	B0DD(k1 + q2, mq))*2./(s12 + S2));
	}
	//COM E10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM E10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//CLHEP::HepLorentzVector q2=k3+k4;
	CLHEP::HepLorentzVector q2 = -(k1+k2+kh);
	double Delta, Sigma, S1, S2, s12, s34;
	S1 = 2.*k1.dot(q2);
	S2 = 2.*k2.dot(q2);
	s12 = 2.*k1.dot(k2);
	//s34 = 2.*k3.dot(k4);
	s34 = q2.m2();
	Delta = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor(-s12D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta +
	12.mqmqS1(s34 + S1)/Delta/Delta -
	4.s12S1pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)C0DD(k2, k1 + q2, mq) -
	s12C0DD(k1, k2, mq) + (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S1C0DD(k1, q2, mq))(1./Delta +
	4.mqmq*S1/Delta/Delta -
	4.s12S1*pow((s34 + S1),2)/Delta/Delta/Delta) +
	C0DD(k1 + k2, q2, mq)(4.s12s34(S1 - S2)/(Delta*Sigma) -
	4.*(s12 -
	2.mqmq)(2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma)) + (B0DD(k1 + q2, mq) -
	B0DD(k1 + k2 + q2, mq))(4.(s34 + S1)/((s12 + S2)*Delta) +
	8.S1(s34 + S1)/Delta/Delta) + (B0DD(q2, mq) -
	B0DD(k1 + k2 + q2, mq) +
	s12C0DD(k1 + k2, q2, mq))(12.s34(2.*s12 + S1 +
	S2)(2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma*Sigma) -
	4.s34(4.s12 + 3.S1 +
	S2)/(Delta*Sigma) +
	8.s12s34(s34(s12 + S2) -
	S1*(s34 +
	S1))/(DeltaDeltaSigma)) + (B0DD(k1 + k2, mq) -
	B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2,
	mq))(12.s12(2.s34 + S1 +
	S2)(2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma*Sigma) +
	8.s12S1(s34(s12 + S2) -
	S1*(s34 +
	S1))/(DeltaDeltaSigma))) + (COM(0.,1.)/(4.M_PIM_PI))((2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma));
	}

	//COM F10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM F10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//CLHEP::HepLorentzVector q2=k3+k4;
	CLHEP::HepLorentzVector q2 = -(k1+k2+kh);
	double Delta, Sigma, S1, S2, s12, s34;
	S1 = 2.*k1.dot(q2);
	S2 = 2.*k2.dot(q2);
	s12 = 2.*k1.dot(k2);
	//s34 = 2.*k3.dot(k4);
	s34 = q2.m2();
	Delta = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor* (s12*D0DD(k1, k2, q2,
	mq)((s34 + S2)/Delta - 4.mq*mq/Delta +
	12.mqmqs34(s12 + S1)/Delta/Delta -
	4.s12pow((s34 + S2),2)/Delta/Delta -
	4.s12S1pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)C0DD(k1, k2 + q2, mq) -
	s12C0DD(k1, k2, mq) - (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S2C0DD(k2, q2, mq))(1./Delta +
	4.mqmq*S1/Delta/Delta -
	4.s12(s34 + S2)/Delta/Delta -
	4.s12S1*pow((s34 + S2),2)/Delta/Delta/Delta) -
	C0DD(k1 + k2, q2, mq)(4.s12s34/(S2Delta) +
	4.s12s34(S2 - S1)/(DeltaSigma) +
	4.*(s12 -
	2.mqmq)(2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma)) - (B0DD(
	k2 + q2, mq) -
	B0DD(k1 + k2 + q2, mq))(4.s34/(S2*Delta) +
	8.s34(s12 + S1)/Delta/Delta) - (B0DD(q2, mq) -
	B0DD(k1 + k2 + q2, mq) +
	s12*C0DD(k1 + k2, q2,
	mq))(-12s34(2s12 + S1 +
	S2)(2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma*Sigma) -
	4.s12s34s34/(S2Delta*Delta) +
	4.s34S1/(Delta*Sigma) -
	4.s34(s12s34(2.*s12 + S2) -
	S1S1(2.*s12 +
	S1))/(DeltaDeltaSigma)) - (B0DD(k1 + k2, mq) -
	B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)C0DD(k1 + k2, q2, mq))(-12.s12(2.*s34 + S1 +
	S2)(2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma*Sigma) +
	8.s12(2.s34 + S1)/(DeltaSigma) -
	8.s12s34(2.s12s34 - S1(S1 + S2) +
	s12*(S2 -
	S1))/(DeltaDeltaSigma))) + (COM(0.,1.)/(4.M_PIM_PI))((2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma));

	}

	//COM G10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM G10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//CLHEP::HepLorentzVector q2=k3+k4;
	CLHEP::HepLorentzVector q2 = -(k1+k2+kh);
	double Delta, S1, S2, s12, s34;
	S1 = 2.*k1.dot(q2);
	S2 = 2.*k2.dot(q2);
	s12 = 2.*k1.dot(k2);
	//s34 = 2.*k3.dot(k4);
	s34 = q2.m2();
	Delta = s12s34 - S1S2;

	return looprwfactor* (-D0DD(k1, q2, k2, mq)*(1. +
	4.S1mqmq/Delta) + ((s12 + S1)C0DD(k1,
	k2 + q2, mq) -
	S1C0DD(k1, q2, mq))(1./Delta +
	4.S1mqmq/Delta/Delta) - ((s12 + S2)C0DD(k1 + q2,
	k2, mq) - S2C0DD(k2, q2, mq))(1./Delta +
	4.S1mq*mq/Delta/Delta) + (B0DD(k1 + k2 + q2, mq) -
	B0DD(k1 + q2, mq))4.(s34 +
	S1)/(Delta*(s12 + S2)) + (B0DD(q2, mq) -
	B0DD(k2 + q2, mq))4.s34/(Delta*S2));
	}

	//COM H1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM H1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//return E1(k1,k2,k3,k4,mq)+F1(k1,k2,k3,k4,mq)+G1(k1,k2,k3,k4,mq);
	return E1(k1,k2,kh,mq)+F1(k1,k2,kh,mq)+G1(k1,k2,kh,mq);
	}
	//COM H4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM H4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//return E4(k1,k2,k3,k4,mq)+F4(k1,k2,k3,k4,mq)+G4(k1,k2,k3,k4,mq);
	return E4(k1,k2,kh,mq)+F4(k1,k2,kh,mq)+G4(k1,k2,kh,mq);
	}
	//COM H10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM H10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//return E10(k1,k2,k3,k4,mq)+F10(k1,k2,k3,k4,mq)+G10(k1,k2,k3,k4,mq);
	return E10(k1,k2,kh,mq)+F10(k1,k2,kh,mq)+G10(k1,k2,kh,mq);
	}
	//COM H2(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM H2(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//return -1.*H1(k2,k1,k3,k4,mq);
	return -1.*H1(k2,k1,kh,mq);
	}
	//COM H5(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM H5(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//return -1.*H4(k2,k1,k3,k4,mq);
	return -1.*H4(k2,k1,kh,mq);
	}
	//COM H12(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
	COM H12(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
	{
	//return -1.*H10(k2,k1,k3,k4,mq);
	return -1.*H10(k2,k1,kh,mq);
	}

	// FL and FT functions
	COM FL(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mq)
	{
	CLHEP::HepLorentzVector Q = q1 + q2;
	double detQ2 = q1.m2()q2.m2() - q1.dot(q2)q1.dot(q2);
	return -1./(2.detQ2)((2.-
	3.q1.m2()q2.dot(Q)/detQ2)*(B0DD(q1, mq) -
	B0DD(Q, mq)) + (2. -
	3.q2.m2()q1.dot(Q)/detQ2)*(B0DD(q2, mq) -
	B0DD(Q, mq)) - (4.mqmq + q1.m2() + q2.m2() +
	Q.m2() - 3.q1.m2()q2.m2()Q.m2()/detQ2)C0DD(
	q1, q2, mq) - 2.);
	}
	COM FT(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mq)
	{
	CLHEP::HepLorentzVector Q = q1 + q2;
	double detQ2 = q1.m2()q2.m2() - q1.dot(q2)q1.dot(q2);
	return -1./(2.detQ2)(Q.m2()(B0DD(q1, mq) + B0DD(q2, mq) - 2.B0DD(Q, mq) -
	2.q1.dot(q2)C0DD(q1, q2, mq)) + (q1.m2() -
	q2.m2()) *(B0DD(q1, mq) - B0DD(q2, mq))) -
	q1.dot(q2)*FL(q1, q2, mq);
	}

	CLHEP::HepLorentzVector ParityFlip(CLHEP::HepLorentzVector p)
	{
	CLHEP::HepLorentzVector flippedVector;
	flippedVector.setE(p.e());
	flippedVector.setX(-p.x());
	flippedVector.setY(-p.y());
	flippedVector.setZ(-p.z());
	return flippedVector;
	}
	/// @brief HC amp for qg->qgH with finite top (i.e. j^{++}_H)
	void g_gH_HC(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1,
	CLHEP::HepLorentzVector pH, double mq, current &retAns)
	{
	current cura1,pacur,p1cur,pHcur,conjeps1,conjepsH1,epsa,epsHa,epsHapart1,
	epsHapart2,conjepsH1part1,conjepsH1part2;
	COM ang1a,sqa1;

	const double F = 4.mqmq/v;
	// Easier to have the whole thing as current object so I can use cdot functionality.
	// Means I need to write pa,p1 as current objects
	to_current(pa, pacur);
	to_current(p1,p1cur);
	to_current(pH,pHcur);
	bool gluonforward = true;
	if(pa.z() < 0)
	gluonforward = false;
	//HEJ gauge
	jio(pa,false,p1,false,cura1);

	if(gluonforward){
	// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
	ang1a = sqrt(pa.plus()p1.minus())(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
	// sqrt(2pa_-/p1_-)p1_perp/abs(p1_perp)
	sqa1 = sqrt(pa.plus()p1.minus())(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
	} else {
	ang1a = sqrt(pa.minus()*p1.plus());
	sqa1 = sqrt(pa.minus()*p1.plus());
	}


	const double prop = (pa-p1-pH).m2();

	cmult(-1./sqrt(2)/ang1a,cura1,conjeps1);
	cmult(1./sqrt(2)/sqa1,cura1,epsa);

	const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
	const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();

	const COM h4 = H4(p1,-pa,pH,mq);
	const COM h5 = H5(p1,-pa,pH,mq);
	const COM h10 = H10(p1,-pa,pH,mq);
	const COM h12 = H12(p1,-pa,pH,mq);


	cmult(Fta*pa.dot(pH), epsa, epsHapart1);
	cmult(-1.Ftacdot(pHcur,epsa), pacur, epsHapart2);
	cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
	cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
	cadd(epsHapart1, epsHapart2, epsHa);
	cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
	const COM aH1 = cdot(pHcur, cura1);

	current T1,T2,T3,T4,T5,T6,T7,T8,T9,T10;

	if(gluonforward){
	cmult(sqrt(2.)sqrt(p1.plus()/pa.plus())prop/sqa1, conjepsH1, T1);
	cmult(-sqrt(2.)sqrt(pa.plus()/p1.plus())prop/ang1a, epsHa, T2);
	}
	else{
	cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())
	((p1.x()-COM(0.,1.)p1.y())/p1.perp())*prop/sqa1, conjepsH1, T1);
	cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())
	((p1.x()-COM(0.,1.)p1.y())/p1.perp())*prop/ang1a, epsHa, T2);
	}

	cmult(sqrt(2.)/ang1a*aH1, epsHa, T3);
	cmult(sqrt(2.)/sqa1*aH1, conjepsH1, T4);

	cmult(-sqrt(2.)Ftapa.dot(p1)*aH1/sqa1, conjeps1, T5);
	cmult(-sqrt(2.)Ft1pa.dot(p1)*aH1/ang1a, epsa, T6);

	cmult(-aH1/sqrt(2.)/sqa1h48.COM(0.,1.)M_PI*M_PI, conjeps1, T7);
	cmult(aH1/sqrt(2.)/ang1ah58.COM(0.,1.)M_PI*M_PI, epsa, T8);
	cmult(aH1aH1/2./ang1a/sqa1h108.COM(0.,1.)M_PIM_PI, pacur, T9);
	cmult(-aH1aH1/2./ang1a/sqa1h128.COM(0.,1.)M_PIM_PI, p1cur, T10);

	current ans;
	for(int i=0;i<4;i++)
	{
	ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5[i]+T6[i]+T7[i]+T8[i]+T9[i]+T10[i];
	}

	retAns[0] = F/prop*ans[0];
	retAns[1] = F/prop*ans[1];
	retAns[2] = F/prop*ans[2];
	retAns[3] = F/prop*ans[3];
	}

	/// @brief HNC amp for qg->qgH with finite top (i.e. j^{+-}_H)
	void g_gH_HNC(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH, double mq, current &retAns)
	{
	const double F = 4.mqmq/v;
	COM ang1a,sqa1;
	current conjepsH1,epsHa,p1cur,pacur,pHcur,conjeps1,epsa,paplusp1cur,
	p1minuspacur,cur1a,cura1,epsHapart1,epsHapart2,conjepsH1part1,
	conjepsH1part2;
	// Find here if pa, meaning the gluon, is forward or backward
	bool gluonforward = true;
	if(pa.z() < 0)
	gluonforward = false;

	jio(pa,true,p1,true,cura1);
	j(p1,true,pa,true,cur1a);

	to_current(pa,pacur);
	to_current(p1,p1cur);
	to_current(pH,pHcur);
	to_current(pa+p1,paplusp1cur);
	to_current(p1-pa,p1minuspacur);
	const COM aH1 = cdot(pHcur,cura1);
	const COM oneHa = std::conj(aH1); // = cdot(pHcur,cur1a)

	if(gluonforward){
	// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
	ang1a = sqrt(pa.plus()p1.minus())(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
	// sqrt(2pa_-/p1_-)p1_perp/abs(p1_perp)
	sqa1 = sqrt(pa.plus()p1.minus())(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
	}
	else {
	ang1a = sqrt(pa.minus()*p1.plus());
	sqa1 = sqrt(pa.minus()*p1.plus());
	}

	const double prop = (pa-p1-pH).m2();

	cmult(1./sqrt(2)/sqa1, cur1a, epsa);
	cmult(-1./sqrt(2)/sqa1, cura1, conjeps1);
	const COM phase = cdot(conjeps1, epsa);
	const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
	const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
	const COM Falpha = FT(p1-pa,pa-p1-pH,mq);
	const COM Fbeta = FL(p1-pa,pa-p1-pH,mq);

	const COM h1 = H1(p1,-pa, pH, mq);
	const COM h2 = H2(p1,-pa, pH, mq);
	const COM h4 = H4(p1,-pa, pH, mq);
	const COM h5 = H5(p1,-pa, pH, mq);
	const COM h10 = H10(p1,-pa, pH, mq);
	const COM h12 = H12(p1,-pa, pH, mq);

	cmult(Fta*pa.dot(pH), epsa, epsHapart1);
	cmult(-1.Ftacdot(pHcur,epsa), pacur, epsHapart2);
	cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
	cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
	cadd(epsHapart1, epsHapart2, epsHa);
	cadd(conjepsH1part1, conjepsH1part2, conjepsH1);

	current T1,T2,T3,T4,T5a,T5b,T6,T7,T8a,T8b,T9,T10,T11a,
	T11b,T12a,T12b,T13;

	if(gluonforward){
	cmult(sqrt(2.)sqrt(p1.plus()/pa.plus())prop/sqa1, conjepsH1, T1);
	cmult(-sqrt(2.)sqrt(pa.plus()/p1.plus())prop/sqa1, epsHa, T2);
	}
	else{
	cmult(-sqrt(2.)sqrt(p1.minus()/pa.minus())((p1.x()-COM(0.,1.)*p1.y())/p1.perp())
	*prop/sqa1, conjepsH1, T1);
	cmult(sqrt(2.)sqrt(pa.minus()/p1.minus())((p1.x()+COM(0.,1.)*p1.y())/p1.perp())
	*prop/sqa1, epsHa, T2);
	}

	const COM boxdiagFact = 8.COM(0.,1.)M_PI*M_PI;

	cmult(aH1*sqrt(2.)/sqa1, epsHa, T3);
	cmult(oneHa*sqrt(2.)/sqa1, conjepsH1, T4);
	cmult(-2.phaseFta*pa.dot(pH), p1cur, T5a);
	cmult(2.phaseFt1*p1.dot(pH), pacur, T5b);
	cmult(-sqrt(2.)Ftap1.dot(pa)*oneHa/sqa1, conjeps1, T6);
	cmult(-sqrt(2.)Ft1pa.dot(p1)*aH1/sqa1, epsa, T7);

	cmult(-boxdiagFactphaseh2, pacur, T8a);
	cmult(boxdiagFactphaseh1, p1cur, T8b);
	cmult(boxdiagFactaH1/sqrt(2.)/sqa1h5, epsa, T9);
	cmult(-boxdiagFactoneHa/sqrt(2.)/sqa1h4, conjeps1, T10);
	cmult(boxdiagFactaH1oneHa/2./sqa1/sqa1*h10, pacur, T11a);
	cmult(-boxdiagFactaH1oneHa/2./sqa1/sqa1*h12, p1cur, T11b);

	cmult(-phase/(pa-p1).m2()Falpha(p1-pa).dot(pa-p1-pH), paplusp1cur, T12a);
	cmult(phase/(pa-p1).m2()Falpha(pa+p1).dot(pa-p1-pH), p1minuspacur, T12b);
	cmult(-phaseFbeta(pa-p1-pH).m2(), paplusp1cur, T13);

	current ans;
	for(int i=0;i<4;i++)
	{
	ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5a[i]+T5b[i]+T6[i]+T7[i]+T8a[i]+T8b[i]+T9[i]+T10[i]+T11a[i]+T11b[i]+T12a[i]+T12b[i]+T13[i];
	}

	retAns[0] = F/prop*ans[0];
	retAns[1] = F/prop*ans[1];
	retAns[2] = F/prop*ans[2];
	retAns[3] = F/prop*ans[3];
	}

	} // namespace anonymous
	// JDC - new amplitude with Higgs emitted close to gluon with full mt effects. Keep usual HEJ-style function call
	double MH2gq_outsideH(CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pH, double mq, bool includeBottom, double mq2)
	{

	current cur2bplus,cur2bminus, cur2bplusFlip, cur2bminusFlip;
	current retAns,retAnsb;
	j(p2out,true,p2in,true,cur2bplus);
	j(p2out,false,p2in,false,cur2bminus);
	j(ParityFlip(p2out),true,ParityFlip(p2in),true,cur2bplusFlip);
	j(ParityFlip(p2out),false,ParityFlip(p2in),false,cur2bminusFlip);

	COM app1,app2,apm1,apm2;
	COM app3, app4, apm3, apm4;

	if(!includeBottom)
	{
	g_gH_HC(p1in,p1out,pH,mq,retAns);
	app1=cdot(retAns,cur2bplus);
	app2=cdot(retAns,cur2bminus);

	g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
	app3=cdot(retAns,cur2bplusFlip);
	app4=cdot(retAns,cur2bminusFlip);

	// And non-conserving bits
	g_gH_HNC(p1in,p1out,pH,mq,retAns);
	apm1=cdot(retAns,cur2bplus);
	apm2=cdot(retAns,cur2bminus);

	g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
	apm3=cdot(retAns,cur2bplusFlip);
	apm4=cdot(retAns,cur2bminusFlip);
	} else {
	g_gH_HC(p1in,p1out,pH,mq,retAns);
	g_gH_HC(p1in,p1out,pH,mq2,retAnsb);
	app1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
	app2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);

	g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
	g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,retAnsb);
	app3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
	app4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);

	// And non-conserving bits
	g_gH_HNC(p1in,p1out,pH,mq,retAns);
	g_gH_HNC(p1in,p1out,pH,mq2,retAnsb);
	apm1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
	apm2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);

	g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
	g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,retAnsb);
	apm3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
	apm4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
	}

	return abs2(app1) + abs2(app2) + abs2(app3) + abs2(app4) + abs2(apm1)
	+ abs2(apm2) + abs2(apm3) + abs2(apm4);
	}
	#endif // RHEJ_BUILD_WITH_QCDLOOP

	double C2gHgm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH)
	{
	static double A=1./(3.M_PIv);
	// Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta
	double s12,p1p,p2p;
	COM p1perp,p3perp,phperp;
	// Determine first whether this is the case p1p\sim php>>p3p og the opposite
	s12=p1.invariantMass2(-p2);
	if (p2.pz()>0.) { // case considered in hep-ph/0301013
	p1p=p1.plus();
	p2p=p2.plus();
	} else { // opposite case
	p1p=p1.minus();
	p2p=p2.minus();
	}
	p1perp=p1.px()+COM(0,1)*p1.py();
	phperp=pH.px()+COM(0,1)*pH.py();
	p3perp=-(p1perp+phperp);

	COM temp=COM(0,1)A/(2.s12)(p2p/p1pconj(p1perp)p3perp+p1p/p2pp1perp*conj(p3perp));
	temp=temp*conj(temp);
	return temp.real();
	}

	double C2gHgp(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH)
	{
	static double A=1./(3.M_PIv);
	// Implements Eq. (4.23) in hep-ph/0301013
	double s12,php,p1p,phm;
	COM p1perp,p3perp,phperp;
	// Determine first whether this is the case p1p\sim php>>p3p or the opposite
	s12=p1.invariantMass2(-p2);
	if (p2.pz()>0.) { // case considered in hep-ph/0301013
	php=pH.plus();
	phm=pH.minus();
	p1p=p1.plus();
	} else { // opposite case
	php=pH.minus();
	phm=pH.plus();
	p1p=p1.minus();
	}
	p1perp=p1.px()+COM(0,1)*p1.py();
	phperp=pH.px()+COM(0,1)*pH.py();
	p3perp=-(p1perp+phperp);

	COM temp=-COM(0,1)A/(2.s12)( conj(p1perpp3perp)*pow(php/p1p,2)/(1.+php/p1p)
	+s12(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1pphm)
	-pow(conj(p3perp)
	+(1.+php/p1p)conj(p1perp),2)/((1.+php/p1p)(pH.m2()+2.*p1.dot(pH)))) );
	temp=temp*conj(temp);
	return temp.real();
	}

	double C2qHqm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH)
	{
	static double A=1./(3.M_PIv);
	// Implements Eq. (4.22) in hep-ph/0301013
	double s12,p2p,p1p;
	COM p1perp,p3perp,phperp;
	// Determine first whether this is the case p1p\sim php>>p3p or the opposite
	s12=p1.invariantMass2(-p2);
	if (p2.pz()>0.) { // case considered in hep-ph/0301013
	p2p=p2.plus();
	p1p=p1.plus();
	} else { // opposite case
	p2p=p2.minus();
	p1p=p1.minus();
	}
	p1perp=p1.px()+COM(0,1)*p1.py();
	phperp=pH.px()+COM(0,1)*pH.py();
	p3perp=-(p1perp+phperp);

	COM temp=A/(2.s12)( sqrt(p2p/p1p)p3perpconj(p1perp)
	+sqrt(p1p/p2p)p1perpconj(p3perp) );
	temp=temp*conj(temp);
	return temp.real();
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 8:32 PM (1 d, 4 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806049
Default Alt Text: (190 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions