diff --git a/src/Wjets.cc b/src/Wjets.cc
index 57dca4f..ba1da6c 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,1137 +1,1666 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Wjets.hh"
#include <array>
#include <iostream>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/jets.hh"
#include "HEJ/Tensor.hh"
using HEJ::Tensor;
using HEJ::init_sigma_index;
using HEJ::metric;
using HEJ::rank3_current;
using HEJ::rank5_current;
using HEJ::eps;
using HEJ::to_tensor;
using HEJ::Helicity;
using HEJ::angle;
using HEJ::square;
using HEJ::flip;
using HEJ::ParticleProperties;
namespace helicity = HEJ::helicity;
namespace { // Helper Functions
// FKL W Helper Functions
double WProp (const HLV & plbar, const HLV & pl, ParticleProperties const & wprop){
COM propW = COM(0.,-1.)/( (pl+plbar).m2() - wprop.mass*wprop.mass
+ COM(0.,1.)*wprop.mass*wprop.width);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
namespace {
// FKL current including W emission off negative helicities
// See eq. (87) {eq:jW-} in developer manual
// Note that the terms are rearranged
Tensor<1> jW_minus(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl
){
using HEJ::helicity::minus;
const double tWin = (pa-pl-plbar).m2();
const double tWout = (p1+pl+plbar).m2();
// C++ arithmetic operators are evaluated left-to-right,
// so the following first computes complex scalar coefficients,
// which then multiply a current, reducing the number
// of multiplications
return 2.*(
+ angle(p1, pl)*square(p1, plbar)/tWout
+ square(pa, plbar)*angle(pa, pl)/tWin
)*HEJ::current(p1, pa, helicity::minus)
+ 2.*angle(p1, pl)*square(pl, plbar)/tWout
*HEJ::current(pl, pa, helicity::minus)
+ 2.*square(pa, plbar)*angle(pl, plbar)/tWin
*HEJ::current(p1, plbar, helicity::minus);
}
}
// FKL current including W emission
// see eqs. (87), (88) {eq:jW-}, {eq:jW+} in developer manual
Tensor<1> jW(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl,
Helicity h
){
if(h == helicity::minus) {
return jW_minus(pa, p1, plbar, pl);
}
return jW_minus(pa, p1, pl, plbar).complex_conj();
}
/**
+ * @brief Contraction of the \tilde{U}_1 tensor
+ *
+ * This is the contraction of the tensor defined in eq:U1tensor
+ * in the developer manual with the uno gluon polarisation vector,
+ * the leptonic and the partonic current (see eq:wunocurrent) in the
+ * developer manual)
+ */
+ COM U1contr(
+ HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
+ HLV const & pa, Helicity h1,
+ HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
+ ) {
+ const HLV pW = pl+plbar;
+ const HLV q1g = pa-pW-p1-pg;
+ const HLV q1 = pa-p1-pW;
+ const HLV q2 = p2-pb;
+
+ const double taW = (pa-pW).m2();
+ const double s1W = (p1+pW).m2();
+ const double s1gW = (p1+pW+pg).m2();
+ const double s1g = (p1+pg).m2();
+
+ if(h1 == helicity::plus) {
+ if(h2 == helicity::plus) {
+ if(hg == helicity::plus) {
+ // helicity +++
+ return (4.*sqrt(2.)*(-1.*taW*angle(pa,pb)*square(p1,plbar)*
+ (s1W*angle(p1,pg)*angle(pg,pl)*square(p1,p2)*square(p2,pg) -
+ 1.*(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar))*
+ ((s1g + s1W)*angle(p1,pg)*square(p1,p2) +
+ s1g*(angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar)))) +
+ s1gW*s1W*angle(p1,pg)*angle(pa,pl)*pow(square(p1,p2),2.)*
+ (angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))))/
+ (s1g*s1gW*s1W*taW*square(p2,pg));
+ }
+ else {
+ // helicity ++-
+ return (sqrt(2.)*(taW*angle(pa,pb)*(4.*s1g*square(p1,plbar)*
+ (angle(p1,pl)*square(p1,pg)*(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
+ angle(pb,plbar)*square(p2,plbar)) +
+ angle(pl,plbar)*(angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg) +
+ angle(pb,pl)*square(p2,pl) + angle(pb,plbar)*square(p2,plbar))*square(pg,plbar)) +
+ square(p1,pg)*(4.*angle(p1,pb)*square(p1,plbar)*
+ ((s1g + s1W)*angle(p1,pl)*square(p1,p2) - 1.*s1W*angle(pg,pl)*square(p2,pg) +
+ s1W*angle(pl,plbar)*square(p2,plbar)) -
+ 4.*s1W*angle(pb,pg)*(angle(p1,pl)*square(p1,p2) - 1.*angle(pg,pl)*square(p2,pg) +
+ angle(pl,plbar)*square(p2,plbar))*square(pg,plbar))) +
+ 4.*s1gW*s1W*angle(pa,pl)*square(p1,pg)*
+ (angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg))*
+ (angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))))/
+ (s1g*s1gW*s1W*taW*angle(pb,pg));
+ }
+ }
+ else {
+ if(hg == helicity::plus) {
+ // helicity +-+
+ return (sqrt(2.)*(4.*taW*angle(p2,pa)*square(p1,plbar)*
+ (s1W*angle(p1,pg)*angle(pg,pl)*square(p1,pb)*square(pb,pg) -
+ 1.*(angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar))*
+ ((s1g + s1W)*angle(p1,pg)*square(p1,pb) +
+ s1g*(angle(pg,pl)*square(pb,pl) + angle(pg,plbar)*square(pb,plbar)))) -
+ 4.*s1gW*s1W*angle(p1,pg)*angle(pa,pl)*pow(square(p1,pb),2.)*
+ (angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
+ (s1g*s1gW*s1W*taW*square(pb,pg));
+ }
+ else {
+ // helicity +--
+ return (-1.*sqrt(2.)*(taW*angle(p2,pa)*(4.*s1g*square(p1,plbar)*
+ (angle(p1,pl)*square(p1,pg)*(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
+ angle(p2,plbar)*square(pb,plbar)) +
+ angle(pl,plbar)*(angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg) +
+ angle(p2,pl)*square(pb,pl) + angle(p2,plbar)*square(pb,plbar))*square(pg,plbar)) +
+ square(p1,pg)*(4.*angle(p1,p2)*square(p1,plbar)*
+ ((s1g + s1W)*angle(p1,pl)*square(p1,pb) - 1.*s1W*angle(pg,pl)*square(pb,pg) +
+ s1W*angle(pl,plbar)*square(pb,plbar)) -
+ 4.*s1W*angle(p2,pg)*(angle(p1,pl)*square(p1,pb) - 1.*angle(pg,pl)*square(pb,pg) +
+ angle(pl,plbar)*square(pb,plbar))*square(pg,plbar))) +
+ 4.*s1gW*s1W*angle(pa,pl)*square(p1,pg)*
+ (angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg))*
+ (angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
+ (s1g*s1gW*s1W*taW*angle(p2,pg));
+ }
+ }
+ }
+ else {
+ if(h2 == helicity::plus) {
+ if(hg == helicity::plus) {
+ // helicity -++
+ return (
+ sqrt(2.)*(-4.*s1gW*s1W*angle(p1,pg)*(angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg))*
+ (angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*square(pa,plbar) +
+ taW*angle(pg,pl)*square(p2,pa)*(4.*s1g*angle(p1,pl)*
+ (angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
+ angle(pb,plbar)*square(p2,plbar))*square(pl,plbar) +
+ 4.*s1W*angle(p1,pg)*square(p2,pg)*
+ (angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pg)*square(pg,plbar) -
+ 1.*angle(pb,pl)*square(pl,plbar))) -
+ 4.*taW*angle(p1,pg)*angle(p1,pl)*square(p2,pa)*
+ (square(p1,plbar)*((s1g + s1W)*angle(p1,pb)*square(p1,p2) +
+ s1g*(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
+ angle(pb,plbar)*square(p2,plbar))) -
+ 1.*s1W*square(p1,p2)*(angle(pb,pg)*square(pg,plbar) + angle(pb,pl)*square(pl,plbar))))
+ )/(s1g*s1gW*s1W*taW*square(p2,pg));
+ }
+ else {
+ // helicity -+-
+ return (-1.*sqrt(2.)*(4.*s1W*angle(p1,pb)*square(p1,pg)*
+ (s1gW*angle(p1,pb)*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
+ square(pa,plbar) - 1.*taW*angle(p1,pl)*angle(pb,pg)*square(p2,pa)*square(pg,plbar)) +
+ taW*angle(p1,pl)*square(p2,pa)*((s1g + s1W)*angle(p1,pb)*square(p1,pg) +
+ s1g*(angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar)))*
+ (4.*angle(p1,pb)*square(p1,plbar) - 4.*angle(pb,pl)*square(pl,plbar))))/
+ (s1g*s1gW*s1W*taW*angle(pb,pg));
+ }
+ }
+ else {
+ if(hg == helicity::plus) {
+ // helicity --+
+ return (sqrt(2.)*(4.*s1gW*s1W*angle(p1,pg)*square(pa,plbar)*
+ (angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg))*
+ (angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)) +
+ taW*square(pa,pb)*(-1.*angle(pg,pl)*
+ (4.*s1g*angle(p1,pl)*(angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg) +
+ angle(p2,pl)*square(pb,pl) + angle(p2,plbar)*square(pb,plbar))*square(pl,plbar) +
+ 4.*s1W*angle(p1,pg)*square(pb,pg)*
+ (angle(p1,p2)*square(p1,plbar) - 1.*angle(p2,pg)*square(pg,plbar) -
+ 1.*angle(p2,pl)*square(pl,plbar))) +
+ 4.*angle(p1,pg)*angle(p1,pl)*(square(p1,plbar)*
+ ((s1g + s1W)*angle(p1,p2)*square(p1,pb) +
+ s1g*(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
+ angle(p2,plbar)*square(pb,plbar))) -
+ 1.*s1W*square(p1,pb)*(angle(p2,pg)*square(pg,plbar) + angle(p2,pl)*square(pl,plbar)))))
+ )/(s1g*s1gW*s1W*taW*square(pb,pg));
+ }
+ else {
+ // helicity ---
+ return (sqrt(2.)*(s1W*angle(p1,p2)*square(p1,pg)*
+ (4.*s1gW*angle(p1,p2)*square(pa,plbar)*
+ (angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)) -
+ 4.*taW*angle(p1,pl)*angle(p2,pg)*square(pa,pb)*square(pg,plbar)) +
+ taW*angle(p1,pl)*square(pa,pb)*((s1g + s1W)*angle(p1,p2)*square(p1,pg) +
+ s1g*(angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))*
+ (4.*angle(p1,p2)*square(p1,plbar) - 4.*angle(p2,pl)*square(pl,plbar))))/
+ (s1g*s1gW*s1W*taW*angle(p2,pg));
+ }
+ }
+ }
+ }
+
+ /**
+ * @brief Contraction of the \tilde{U}_2 tensor
+ *
+ * This is the contraction of the tensor defined in eq:U2tensor
+ * in the developer manual with the uno gluon polarisation vector,
+ * the leptonic and the partonic current (see eq:wunocurrent) in the
+ * developer manual)
+ */
+ COM U2contr(
+ HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
+ HLV const & pa, Helicity h1,
+ HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
+ ) {
+ const HLV pW = pl+plbar;
+ const HLV q1g=pa-pW-p1-pg;
+ const HLV q1 = pa-p1-pW;
+ const HLV q2 = p2-pb;
+
+ const double taW = (pa-pW).m2();
+ const double s1W = (p1+pW).m2();
+ const double tag = (pa-pg).m2();
+ const double taWg = (pa-pW-pg).m2();
+
+ if(h1 == helicity::plus) {
+ if(h2 == helicity::plus) {
+ if(hg == helicity::plus) {
+ // helicity +++
+ return (-4.*sqrt(2.)*(taW*angle(pa,pg)*(taWg*square(p1,plbar)*
+ (angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg))*
+ (angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) -
+ 1.*s1W*angle(pg,pl)*square(p1,p2)*square(p2,pg)*
+ (angle(pa,pb)*square(pa,plbar) + angle(pb,pg)*square(pg,plbar) +
+ angle(pb,pl)*square(pl,plbar))) +
+ s1W*angle(pa,pl)*square(p1,p2)*(tag*
+ (angle(pa,pg)*(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
+ angle(pb,plbar)*square(p2,plbar))*square(pa,plbar) +
+ angle(pg,pl)*(angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg) +
+ angle(pb,pl)*square(p2,pl) + angle(pb,plbar)*square(p2,plbar))*square(pl,plbar)) +
+ angle(pa,pg)*square(p2,pa)*((tag + taW)*angle(pa,pb)*square(pa,plbar) +
+ taW*(angle(pb,pg)*square(pg,plbar) + angle(pb,pl)*square(pl,plbar))))))/
+ (s1W*tag*taW*taWg*square(p2,pg));
+ }
+ else {
+ // helicity ++-
+ return (-4.*sqrt(2.)*(s1W*tag*angle(pa,pl)*square(p1,p2)*
+ (angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar))*
+ (angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)) -
+ 1.*angle(pa,pb)*square(pa,pg)*(taW*taWg*angle(pa,pb)*square(p1,plbar)*
+ (angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) +
+ s1W*angle(pa,pl)*square(p1,p2)*
+ (taW*angle(pb,pg)*square(pg,plbar) +
+ (tag + taW)*(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))))))/
+ (s1W*tag*taW*taWg*angle(pb,pg));
+ }
+ }
+ else {
+ if(hg == helicity::plus) {
+ // helicity +-+
+ return (-4.*sqrt(2.)*(taW*angle(pa,pg)*(taWg*square(p1,plbar)*
+ (angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg))*
+ (angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar)) -
+ 1.*s1W*square(p1,pb)*(angle(pa,pl)*square(pa,pb) + angle(pg,pl)*square(pb,pg))*
+ (angle(p2,pg)*square(pg,plbar) + angle(p2,pl)*square(pl,plbar))) +
+ s1W*square(p1,pb)*(tag*angle(pa,pl)*
+ (angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
+ angle(p2,plbar)*square(pb,plbar))*
+ (angle(pa,pg)*square(pa,plbar) + angle(pg,pl)*square(pl,plbar)) +
+ angle(p2,pa)*(angle(pa,pg)*square(pa,plbar)*
+ ((tag + taW)*angle(pa,pl)*square(pa,pb) + taW*angle(pg,pl)*square(pb,pg)) +
+ tag*angle(pa,pl)*angle(pg,pl)*square(pa,pb)*square(pl,plbar)))))/
+ (s1W*tag*taW*taWg*square(pb,pg));
+ }
+ else {
+ // helicity +--
+ return (sqrt(2.)*(4.*taW*angle(p2,pa)*square(pa,pg)*
+ (taWg*angle(p2,pa)*square(p1,plbar)*
+ (angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar)) -
+ 1.*s1W*angle(p2,pg)*angle(pa,pl)*square(p1,pb)*square(pg,plbar)) +
+ s1W*angle(pa,pl)*square(p1,pb)*((tag + taW)*angle(p2,pa)*square(pa,pg) +
+ tag*(angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))*
+ (4.*angle(p2,pa)*square(pa,plbar) - 4.*angle(p2,pl)*square(pl,plbar))))/
+ (s1W*tag*taW*taWg*angle(p2,pg));
+ }
+ }
+ }
+ else {
+ if(h2 == helicity::plus) {
+ if(hg == helicity::plus) {
+ // helicity -++
+ return (sqrt(2.)*(-4.*s1W*angle(p1,pb)*(taW*angle(pa,pg)*angle(pg,pl)*square(p2,pa)*square(p2,pg) -
+ 1.*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
+ ((tag + taW)*angle(pa,pg)*square(p2,pa) +
+ tag*(angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))))*square(pa,plbar)
+ + 4.*taW*taWg*angle(p1,pl)*angle(pa,pg)*pow(square(p2,pa),2.)*
+ (angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar))))/
+ (s1W*tag*taW*taWg*square(p2,pg));
+ }
+ else {
+ // helicity -+-
+ return (sqrt(2.)*(4.*s1W*angle(p1,pb)*(tag*angle(pl,plbar)*
+ (angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
+ angle(pb,plbar)*square(p2,plbar))*square(pa,plbar)*square(pg,plbar) +
+ taW*angle(pg,pl)*square(p2,pg)*square(pa,pg)*
+ (angle(pa,pb)*square(pa,plbar) + angle(pb,pg)*square(pg,plbar)) -
+ 1.*square(pa,pg)*((tag*angle(pa,pl)*
+ (angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
+ angle(pb,plbar)*square(p2,plbar)) +
+ angle(pa,pb)*((tag + taW)*angle(pa,pl)*square(p2,pa) +
+ taW*angle(pl,plbar)*square(p2,plbar)))*square(pa,plbar) +
+ taW*angle(pb,pg)*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
+ square(pg,plbar))) - 4.*taW*taWg*angle(p1,pl)*
+ (angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg))*square(pa,pg)*
+ (angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar))))/
+ (s1W*tag*taW*taWg*angle(pb,pg));
+ }
+ }
+ else {
+ if(hg == helicity::plus) {
+ // helicity --+
+ return (4.*sqrt(2.)*(angle(p1,p2)*(taW*taWg*angle(p1,pl)*angle(pa,pg)*pow(square(pa,pb),2.)*
+ square(p1,plbar) + s1W*square(pa,plbar)*
+ (tag*(angle(pg,pl)*square(pb,pl) + angle(pg,plbar)*square(pb,plbar))*
+ (-1.*angle(pa,pl)*square(pa,pb) + angle(pl,plbar)*square(pb,plbar)) +
+ angle(pa,pg)*square(pa,pb)*((tag + taW)*angle(pa,pl)*square(pa,pb) +
+ taW*angle(pg,pl)*square(pb,pg) - 1.*(tag + taW)*angle(pl,plbar)*square(pb,plbar))))
+ - 1.*taW*taWg*angle(p1,pl)*angle(p2,pl)*angle(pa,pg)*pow(square(pa,pb),2.)*square(pl,plbar))
+ )/(s1W*tag*taW*taWg*square(pb,pg));
+ }
+ else {
+ // helicity ---
+ return (sqrt(2.)*(s1W*angle(p1,p2)*(-4.*square(pa,pg)*square(pa,plbar)*
+ (tag*angle(pa,pl)*(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
+ angle(p2,plbar)*square(pb,plbar)) +
+ angle(p2,pa)*((tag + taW)*angle(pa,pl)*square(pa,pb) +
+ taW*angle(pg,pl)*square(pb,pg) - 1.*taW*angle(pl,plbar)*square(pb,plbar))) +
+ 4.*tag*angle(pl,plbar)*square(pa,plbar)*
+ (angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
+ angle(p2,plbar)*square(pb,plbar))*square(pg,plbar) +
+ 4.*taW*angle(p2,pg)*square(pa,pg)*
+ (angle(pa,pl)*square(pa,pb) + angle(pg,pl)*square(pb,pg) -
+ 1.*angle(pl,plbar)*square(pb,plbar))*square(pg,plbar)) -
+ 4.*taW*taWg*angle(p1,pl)*square(pa,pg)*
+ (angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg))*
+ (angle(p1,p2)*square(p1,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
+ (s1W*tag*taW*taWg*angle(p2,pg));
+ }
+ }
+ }
+ }
+
+ /**
+ * @brief Contraction of the \tilde{L} tensor
+ *
+ * This is the contraction of the tensor defined in eq:Ltensor
+ * in the developer manual with the uno gluon polarisation vector,
+ * the leptonic and the partonic current (see eq:wunocurrent) in the
+ * developer manual)
+ */
+ COM Lcontr(
+ HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
+ HLV const & pa, Helicity h1,
+ HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
+ ) {
+ const HLV pW = pl+plbar;
+ const HLV q1g=pa-pW-p1-pg;
+ const HLV q1 = pa-p1-pW;
+ const HLV q2 = p2-pb;
+
+ const double taW = (pa-pW).m2();
+ const double taW1 = (pa-pW-p1).m2();
+ const double s1W = (p1+pW).m2();
+ const double s2g = (p2+pg).m2();
+ const double sbg = (pb+pg).m2();
+
+ if(h1 == helicity::plus) {
+ if(h2 == helicity::plus) {
+ if(hg == helicity::plus) {
+ // helicity +++
+ return (-2.*sqrt(2.)*(angle(pb,pg)*q1g.m2()*square(p2,pb)*
+ (taW*angle(pa,pb)*square(p1,plbar)*
+ (angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) +
+ s1W*angle(pa,pl)*square(p1,p2)*
+ (angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))) +
+ 2.*sbg*(taW*square(p1,plbar)*(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar))*
+ (angle(pa,pg)*angle(pb,pg)*square(p2,pg) +
+ angle(pa,pb)*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
+ angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))) +
+ s1W*angle(pa,pl)*square(p1,p2)*
+ ((angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) + angle(pg,pl)*square(p2,pl) +
+ angle(pg,plbar)*square(p2,plbar))*
+ (angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)) +
+ angle(pb,pg)*square(p2,pg)*(angle(pa,pg)*square(pa,plbar) +
+ angle(pg,pl)*square(pl,plbar))))))/(s1W*sbg*taW*taW1*square(p2,pg));
+ }
+ else {
+ // helicity ++-
+ return (-1.*sqrt(2.)*(s2g*taW*angle(pa,pb)*square(p1,plbar)*
+ (4.*(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar))*
+ (angle(p1,pb)*square(p1,pg) - 1.*angle(pa,pb)*square(pa,pg) +
+ angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar)) +
+ 4.*angle(pb,pg)*square(p2,pg)*(angle(p1,pl)*square(p1,pg) +
+ angle(pl,plbar)*square(pg,plbar))) +
+ s1W*s2g*angle(pa,pl)*(4.*angle(pb,pg)*square(p1,pg)*square(p2,pg) -
+ 4.*angle(pa,pb)*square(p1,p2)*square(pa,pg) +
+ 4.*square(p1,p2)*(angle(p1,pb)*square(p1,pg) + angle(pb,pl)*square(pg,pl) +
+ angle(pb,plbar)*square(pg,plbar)))*
+ (angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)) -
+ 2.*angle(p2,pb)*q1g.m2()*square(p2,pg)*
+ (taW*angle(pa,pb)*square(p1,plbar)*
+ (angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) +
+ s1W*angle(pa,pl)*square(p1,p2)*
+ (angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)))))/
+ (s1W*s2g*taW*taW1*angle(pb,pg));
+ }
+ }
+ else {
+ if(hg == helicity::plus) {
+ // helicity +-+
+ return (-1.*sqrt(2.)*(2.*taW*square(p1,plbar)*(angle(p1,pl)*square(p1,pb) +
+ angle(pl,plbar)*square(pb,plbar))*
+ (2.*s2g*angle(pa,pg)*(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg)) +
+ angle(p2,pa)*(angle(p2,pg)*q1g.m2()*square(p2,pb) -
+ 2.*s2g*(angle(p1,pg)*square(p1,pb) + angle(pg,pl)*square(pb,pl) +
+ angle(pg,plbar)*square(pb,plbar)))) +
+ s1W*angle(pa,pl)*square(p1,pb)*(4.*s2g*square(pa,plbar)*
+ (angle(pa,pg)*(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg)) -
+ 1.*angle(p2,pa)*(angle(p1,pg)*square(p1,pb) + angle(pg,pl)*square(pb,pl) +
+ angle(pg,plbar)*square(pb,plbar))) +
+ s2g*(-4.*angle(p2,pl)*angle(pa,pg)*square(pa,pb) +
+ 4.*angle(p2,pg)*angle(pg,pl)*square(pb,pg) +
+ 4.*angle(p2,pl)*(angle(p1,pg)*square(p1,pb) + angle(pg,pl)*square(pb,pl) +
+ angle(pg,plbar)*square(pb,plbar)))*square(pl,plbar) +
+ 2.*angle(p2,pg)*q1g.m2()*square(p2,pb)*
+ (angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar)))))/
+ (s1W*s2g*taW*taW1*square(pb,pg));
+ }
+ else {
+ // helicity +--
+ return (sqrt(2.)*(2.*taW*angle(p2,pa)*square(p1,plbar)*
+ (2.*sbg*angle(p2,pg)*square(pb,pg)*
+ (angle(p1,pl)*square(p1,pg) + angle(pl,plbar)*square(pg,plbar)) +
+ (angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar))*
+ (angle(p2,pb)*q1g.m2()*square(pb,pg) +
+ 2.*sbg*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
+ angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))) +
+ 2.*s1W*angle(pa,pl)*(2.*sbg*angle(p1,p2)*square(p1,pb)*square(p1,pg) +
+ 2.*sbg*angle(p2,pa)*square(p1,pb)*square(pa,pg) +
+ angle(p2,pb)*q1g.m2()*square(p1,pb)*square(pb,pg) +
+ 2.*sbg*angle(p2,pg)*square(p1,pg)*square(pb,pg) +
+ 2.*sbg*angle(p2,pl)*square(p1,pb)*square(pg,pl) +
+ 2.*sbg*angle(p2,plbar)*square(p1,pb)*square(pg,plbar))*
+ (angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
+ (s1W*sbg*taW*taW1*angle(p2,pg));
+ }
+ }
+ }
+ else {
+ if(h2 == helicity::plus) {
+ if(hg == helicity::plus) {
+ // helicity -++
+ return (sqrt(2.)*(2.*s1W*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
+ (2.*sbg*angle(p1,pg)*angle(pb,pg)*square(p2,pg) +
+ angle(p1,pb)*(angle(pb,pg)*q1g.m2()*square(p2,pb) +
+ 2.*sbg*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
+ angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))))*square(pa,plbar) +
+ taW*angle(p1,pl)*square(p2,pa)*(4.*sbg*square(p1,plbar)*
+ (angle(p1,pg)*angle(pb,pg)*square(p2,pg) +
+ angle(p1,pb)*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
+ angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))) -
+ 4.*sbg*(angle(pb,pg)*angle(pg,pl)*square(p2,pg) +
+ angle(pb,pl)*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
+ angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar)))*square(pl,plbar) +
+ 2.*angle(pb,pg)*q1g.m2()*square(p2,pb)*
+ (angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar)))))/
+ (s1W*sbg*taW*taW1*square(p2,pg));
+ }
+ else {
+ // helicity -+-
+ return (sqrt(2.)*((angle(p1,pb)*square(pa,plbar)*
+ (((angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
+ (4.*angle(p1,pb)*square(p1,pg) - (2.*angle(p2,pb)*q1g.m2()*square(p2,pg))/s2g -
+ 4.*angle(pa,pb)*square(pa,pg) + 4.*angle(pb,pl)*square(pg,pl) +
+ 4.*angle(pb,plbar)*square(pg,plbar)))/angle(pb,pg) +
+ square(p2,pg)*(-4.*angle(pa,pl)*square(pa,pg) + 4.*angle(pl,plbar)*square(pg,plbar))))/
+ taW + (2.*angle(p1,pl)*(2.*s2g*angle(p1,pb)*square(p1,pg)*square(p2,pa) -
+ 1.*angle(p2,pb)*q1g.m2()*square(p2,pa)*square(p2,pg) +
+ 2.*s2g*(-1.*angle(pa,pb)*square(p2,pa)*square(pa,pg) -
+ 1.*angle(pb,pg)*square(p2,pg)*square(pa,pg) +
+ square(p2,pa)*(angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar))))*
+ (angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar)))/(s1W*s2g*angle(pb,pg))
+ ))/taW1;
+ }
+ }
+ else {
+ if(hg == helicity::plus) {
+ // helicity --+
+ return (2.*sqrt(2.)*(-1.*angle(p1,p2)*(2.*s2g*angle(p1,pg)*square(p1,pb) -
+ 1.*angle(p2,pg)*q1g.m2()*square(p2,pb) +
+ 2.*s2g*(-1.*angle(pa,pg)*square(pa,pb) + angle(pg,pl)*square(pb,pl) +
+ angle(pg,plbar)*square(pb,plbar)))*
+ (taW*angle(p1,pl)*square(p1,plbar)*square(pa,pb) +
+ s1W*square(pa,plbar)*(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)))
+ + taW*angle(p1,pl)*square(pa,pb)*
+ (angle(p2,pg)*(-1.*angle(p2,pl)*q1g.m2()*square(p2,pb) +
+ 2.*s2g*angle(pg,pl)*square(pb,pg)) +
+ 2.*s2g*angle(p2,pl)*(-1.*angle(pa,pg)*square(pa,pb) + angle(pg,pl)*square(pb,pl) +
+ angle(pg,plbar)*square(pb,plbar)))*square(pl,plbar) -
+ 2.*s2g*angle(p1,pg)*(s1W*angle(p2,pg)*square(pa,plbar)*square(pb,pg)*
+ (angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)) +
+ taW*angle(p1,pl)*square(pa,pb)*
+ (angle(p2,pg)*square(p1,plbar)*square(pb,pg) -
+ 1.*angle(p2,pl)*square(p1,pb)*square(pl,plbar)))))/(s1W*s2g*taW*taW1*square(pb,pg));
+ }
+ else {
+ // helicity ---
+ return (2.*sqrt(2.)*(-1.*angle(p1,p2)*(angle(p2,pb)*q1g.m2()*square(pb,pg)*
+ (taW*angle(p1,pl)*square(p1,plbar)*square(pa,pb) +
+ s1W*square(pa,plbar)*(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar))
+ ) + 2.*sbg*(taW*angle(p1,pl)*square(p1,plbar) + s1W*angle(pa,pl)*square(pa,plbar))*
+ (angle(p2,pg)*square(pa,pg)*square(pb,pg) +
+ square(pa,pb)*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
+ angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar))) -
+ 2.*s1W*sbg*angle(pl,plbar)*square(pa,plbar)*
+ (angle(p2,pg)*square(pb,pg)*square(pg,plbar) +
+ square(pb,plbar)*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
+ angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))) +
+ taW*angle(p1,pl)*angle(p2,pl)*(2.*sbg*angle(p2,pg)*square(pa,pg)*square(pb,pg) +
+ square(pa,pb)*(angle(p2,pb)*q1g.m2()*square(pb,pg) +
+ 2.*sbg*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
+ angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar))))*square(pl,plbar)))/
+ (s1W*sbg*taW*taW1*angle(p2,pg));
+ }
+ }
+ }
+ }
+
+ /**
* @brief W+Jets Unordered Contribution Helper Functions
* @returns result of equation (4.1.28) in Helen's Thesis (p.100)
*/
double jM2Wuno(HLV pg, HLV p1,HLV plbar, HLV pl, HLV pa, Helicity h1,
HLV p2, HLV pb, Helicity h2, Helicity pol,
ParticleProperties const & wprop
){
//@TODO Simplify the below (less Tensor class?)
init_sigma_index();
HLV pW = pl+plbar;
HLV q1g=pa-pW-p1-pg;
HLV q1 = pa-p1-pW;
HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double taW1 = (pa-pW-p1).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+pg).m2();
const double s2g = (p2+pg).m2();
const double sbg = (pb+pg).m2();
const double tag = (pa-pg).m2();
const double taWg = (pa-pW-pg).m2();
//use p1 as ref vec in pol tensor
Tensor<1> epsg = eps(pg,p2,pol);
Tensor<1> epsW = HEJ::current(pl,plbar,helicity::minus);
Tensor<1> j2b = HEJ::current(p2,pb,h2);
Tensor<1> Tq1q2 = to_tensor(
(pb/sbg + p2/s2g)*(q1 - pg).m2() + 2*q1 - pg
);
Tensor<3> J31a = rank3_current(p1, pa, h1);
Tensor<2> J2_qaW =J31a.contract((pa-pW)/taW/taW1, 2);
Tensor<2> J2_p1W =J31a.contract((p1+pW)/s1W/taW1, 2);
Tensor<3> L1a = outer(Tq1q2, J2_qaW);
Tensor<3> L1b = outer(Tq1q2, J2_p1W);
Tensor<3> L2a = outer(-2*pg, J2_qaW);
Tensor<3> L2b = outer(-2*pg, J2_p1W);
Tensor<3> L3 = outer(metric(), J2_qaW.contract(2*pg-q1,1)+J2_p1W.contract(2*pg-q1,2));
Tensor<3> L(0.);
Tensor<5> J51a = rank5_current(p1, pa, h1);
Tensor<4> J_qaW = J51a.contract((pa-pW),4);
Tensor<4> J_qag = J51a.contract(pa-pg,4);
Tensor<4> J_p1gW = J51a.contract(p1+pg+pW,4);
Tensor<3> U1a = J_qaW.contract(p1+pg,2);
Tensor<3> U1b = J_p1gW.contract(p1+pg,2);
Tensor<3> U1c = J_p1gW.contract(p1+pW,2);
Tensor<3> U1(0.);
Tensor<3> U2a = J_qaW.contract(pa-pg-pW,2);
Tensor<3> U2b = J_qag.contract(pa-pg-pW,2);
Tensor<3> U2c = J_qag.contract(p1+pW,2);
Tensor<3> U2(0.);
for(int nu=0; nu<4;nu++){
for(int mu=0;mu<4;mu++){
for(int rho=0;rho<4;rho++){
L(nu, mu, rho) = L1a(nu,mu,rho) + L1b(nu,rho,mu)
+ L2a(mu,nu,rho) + L2b(mu,rho,nu) + L3(mu,nu,rho);
U1(nu, mu, rho) = U1a(nu, mu, rho) / (s1g*taW)
+ U1b(nu,rho,mu) / (s1g*s1gW) + U1c(rho,nu,mu) / (s1W*s1gW);
U2(nu,mu,rho) = U2a(mu,nu,rho) / (taWg*taW)
+ U2b(mu,rho,nu) / (taWg*tag) + U2c(rho,mu,nu) / (s1W*tag);
}
}
}
COM X = ((((U1-L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
COM Y = ((((U2+L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
double amp = HEJ::C_A*HEJ::C_F*HEJ::C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
double t1 = q1g.m2();
double t2 = q2.m2();
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
//Divide by t-channels
amp/=(t1*t2);
return amp;
}
+ /**
+ * @brief New implementation of unordered W+Jets current
+ *
+ * See eq:wunocurrent in the developer manual for the definition
+ * of this current
+ *
+ * The aim is to replace jM2Wuno and eventually all uses of the Tensor class.
+ * Explicit tensor contractions are rather slow and the initialisation code
+ * in Tensor.cc is not very transparent.
+ *
+ * At the moment, this function _only_ works for positive-energy spinors,
+ * so jM2Wuno has to be kept for qqbar emission.
+ */
+ double jM2Wuno_pos_energy(
+ HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
+ HLV const & pa, Helicity const h1,
+ HLV const & p2, HLV const & pb, Helicity const h2, Helicity const pol,
+ ParticleProperties const & wprop
+ ){
+ using HEJ::C_A;
+ using HEJ::C_F;
+
+ const COM U1 = U1contr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
+ const COM U2 = U2contr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
+ const COM L = Lcontr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
+
+ const COM X = U1 - L;
+ const COM Y = U2 + L;
+
+ double amp = C_A*C_F*C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
+ amp *= WProp(plbar, pl, wprop);
+
+ const HLV q1g = pa-pl-plbar-p1-pg;
+ const HLV q2 = p2-pb;
+
+ const double t1 = q1g.m2();
+ const double t2 = q2.m2();
+
+ amp /= (t1*t2);
+ return amp;
+ }
+
+
// Relevant Wqqx Helper Functions.
//g->qxqlxl (Calculates gluon to qqx Current. See JV_\mu in WSubleading Notes)
Tensor <1> gtqqxW(HLV pq,HLV pqbar,HLV pl,HLV plbar){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
//blank 3 Gamma Current
Tensor<3> JV23 = rank3_current(pq,pqbar,helicity::minus);
// Components of g->qqW before W Contraction
Tensor<2> JV1 = JV23.contract((pq + pl + plbar),2)/(s2AB);
Tensor<2> JV2 = JV23.contract((pqbar + pl + plbar),2)/(s3AB);
// g->qqW Current. Note Minus between terms due to momentum flow.
// Also note: (-I)^2 from W vert. (I) from Quark prop.
Tensor<1> JVCur = (JV1.contract(ABCur,1) - JV2.contract(ABCur,2))*COM(0.,-1.);
return JVCur;
}
// Helper Functions Calculate the Crossed Contribution
Tensor <2> MCrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar, HLV pl,
HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCross?
// Useful propagator factors
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV q1, q3;
q1=pa;
for(int i=0; i<nabove+1;i++){
q1=q1-partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double tcro1=(q3+pq).m2();
double tcro2=(q1-pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar,helicity::minus);
//Blank 5 gamma Current
Tensor<5> J523 = rank5_current(pq,pqbar,helicity::minus);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_q3q = J523.contract((q3 + pq),2);
Tensor<4> J_2AB = J523.contract((pq + pl + plbar),2);
// Components of Crossed Vertex Contribution
Tensor<3> Xcro1 = J_q3q.contract((pqbar + pl + plbar),3);
Tensor<3> Xcro2 = J_q3q.contract((q1 - pqbar),3);
Tensor<3> Xcro3 = J_2AB.contract((q1 - pqbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xcro1Cont = Xcro1.contract(ABCur,3)/(tcro1*s3AB);
Tensor<2> Xcro2Cont = Xcro2.contract(ABCur,2)/(tcro1*tcro2);
Tensor<2> Xcro3Cont = Xcro3.contract(ABCur,1)/(s2AB*tcro2);
//Initialise the Crossed Vertex Object
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro(mu,nu) = -(-Xcro1Cont(nu,mu)+Xcro2Cont(nu,mu)+Xcro3Cont(nu,mu));
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar,
HLV pl, HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncross?
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV q1, q3;
q1=pa;
for(int i=0; i<nabove+1;i++){
q1=q1-partons.at(i);
}
q3 = q1 - pl - plbar - pq - pqbar;
double tunc1 = (q1-pq).m2();
double tunc2 = (q3+pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
//Blank 5 gamma Current
Tensor<5> J523 = rank5_current(pq,pqbar,helicity::minus);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_2AB = J523.contract((pq + pl + plbar),2);
Tensor<4> J_q1q = J523.contract((q1 - pq),2);
// 2 Contractions taken care of.
Tensor<3> Xunc1 = J_2AB.contract((q3 + pqbar),3);
Tensor<3> Xunc2 = J_q1q.contract((q3 + pqbar),3);
Tensor<3> Xunc3 = J_q1q.contract((pqbar + pl + plbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xunc1Cont = Xunc1.contract(ABCur,1)/(s2AB*tunc2);
Tensor<2> Xunc2Cont = Xunc2.contract(ABCur,2)/(tunc1*tunc2);
Tensor<2> Xunc3Cont = Xunc3.contract(ABCur,3)/(tunc1*s3AB);
//Initialise the Uncrossed Vertex Object
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc(mu,nu) = -(- Xunc1Cont(mu,nu)+Xunc2Cont(mu,nu) +Xunc3Cont(mu,nu));
}
}
return Xunc;
}
// Helper Functions Calculate the g->qqxW (Eikonal) Contributions
Tensor <2> MSymW(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
HLV pl,HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MSym?
double sa2=(pa+pq).m2();
double s12=(p1+pq).m2();
double sa3=(pa+pqbar).m2();
double s13=(p1+pqbar).m2();
double saA=(pa+pl).m2();
double s1A=(p1+pl).m2();
double saB=(pa+plbar).m2();
double s1B=(p1+plbar).m2();
double sb2=(pb+pq).m2();
double s42=(p4+pq).m2();
double sb3=(pb+pqbar).m2();
double s43=(p4+pqbar).m2();
double sbA=(pb+pl).m2();
double s4A=(p4+pl).m2();
double sbB=(pb+plbar).m2();
double s4B=(p4+plbar).m2();
double s23AB=(pl+plbar+pq+pqbar).m2();
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3=q1-pq-pqbar-plbar-pl;
double t1 = (q1).m2();
double t3 = (q3).m2();
// g->qqW Current (Factors of sqrt2 dealt with in this function.)
Tensor<1> JV = gtqqxW(pq,pqbar,pl,plbar);
// 1a gluon emisson Contribution
Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13+s1A+s1B))
+ pa*(t1/(sa2+sa3+saA+saB)) );
Tensor<2> X1aCont = X1a.contract(JV,3);
//4b gluon emission Contribution
Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43+s4A+s4B))
+ pb*(t3/(sb2+sb3+sbA+sbB)) );
Tensor<2> X4bCont = X4b.contract(JV,3);
//Set up each term of 3G diagram.
Tensor<3> X3g1 = outer(q1+pq+pqbar+pl+plbar, metric());
Tensor<3> X3g2 = outer(q3-pq-pqbar-pl-plbar, metric());
Tensor<3> X3g3 = outer(q1+q3, metric());
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(JV,3);
Tensor<2> X3g2Cont = X3g2.contract(JV,2);
Tensor<2> X3g3Cont = X3g3.contract(JV,1);
// XSym is an amalgamation of x1a, X4b and X3g.
// Makes sense from a colour factor point of view.
Tensor<2>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym(mu,nu) = (X3g1Cont(nu,mu) + X3g2Cont(mu,nu) - X3g3Cont(nu,mu))
+ (X1aCont(mu,nu) - X4bCont(mu,nu));
}
}
return Xsym/s23AB;
}
Tensor <2> MCross(HLV pa, HLV pq, HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2=(q1-pqbar).m2();
//Blank 3 gamma Current
Tensor<3> J323 = rank3_current(pq,pqbar,hq);
// 2 gamma current (with 1 contraction already).
Tensor<2> XCroCont = J323.contract((q1-pqbar),2)/(t2);
//Initialise the Crossed Vertex
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro(mu,nu) = XCroCont(nu,mu);
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncross(HLV pa, HLV pq,HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2 = (q1-pq).m2();
//Blank 3 gamma Current
Tensor<3> J323 = rank3_current(pq,pqbar,hq);
// 2 gamma currents (with 1 contraction already).
Tensor<2> XUncCont = J323.contract((q1-pq),2)/t2;
//Initialise the Uncrossed Vertex
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc(mu,nu) = -XUncCont(mu,nu);
}
}
return Xunc;
}
// Helper Functions Calculate the Eikonal Contributions
Tensor <2> MSym(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
std::vector<HLV> partons, Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MsymW?
HLV q1, q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1-pq-pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
double s23 = (pq+pqbar).m2();
double sa2 = (pa+pq).m2();
double sa3 = (pa+pqbar).m2();
double s12 = (p1+pq).m2();
double s13 = (p1+pqbar).m2();
double sb2 = (pb+pq).m2();
double sb3 = (pb+pqbar).m2();
double s42 = (p4+pq).m2();
double s43 = (p4+pqbar).m2();
Tensor<1> qqxCur = HEJ::current(pq, pqbar, hq);
// // 1a gluon emisson Contribution
Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13))+ pa*(t1/(sa2+sa3)));
Tensor<2> X1aCont = X1a.contract(qqxCur,3);
// //4b gluon emission Contribution
Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43)) + pb*(t3/(sb2+sb3)));
Tensor<2> X4bCont = X4b.contract(qqxCur,3);
// New Formulation Corresponding to New Analytics
Tensor<3> X3g1 = outer(q1+pq+pqbar, metric());
Tensor<3> X3g2 = outer(q3-pq-pqbar, metric());
Tensor<3> X3g3 = outer(q1+q3, metric());
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(qqxCur,3);
Tensor<2> X3g2Cont = X3g2.contract(qqxCur,2);
Tensor<2> X3g3Cont = X3g3.contract(qqxCur,1);
Tensor<2>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym(mu, nu) = COM(0,1) * ( (X3g1Cont(nu,mu) + X3g2Cont(mu,nu)
- X3g3Cont(nu,mu)) + (X1aCont(mu,nu) - X4bCont(mu,nu)) );
}
}
return Xsym/s23;
}
//! W+Jets FKL Contributions
/**
* @brief W+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_\mu.
* Handles all possible incoming states.
*/
double jW_j( HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
bool aqlineb, bool /* aqlinef */,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out);
const double WPropfact = WProp(plbar, pl, wprop);
const auto j_W = COM{0,-1}*jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
double Msqr = 0.;
for(const auto h: {plus, minus}) {
const auto j = HEJ::current(p2out, p2in, h);
Msqr += abs2(j_W.contract(j, 1));
}
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*Msqr/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
} // Anonymous Namespace
double ME_W_qQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop);
}
double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true, wprop);
}
double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop);
}
double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true, wprop);
}
double ME_W_qg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param pg Unordered Gluon momenta
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_{uno}_\mu. Ie, unordered with W emission opposite side.
* Handles all possible incoming states.
*/
double jW_juno(HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out,
HLV p2in, HLV pg, bool aqlineb, bool aqlinef,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out-pg);
const HLV q3=-(p2in-p2out);
const Helicity fhel = aqlinef?plus:minus;
const auto j_W = jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
const auto mj2p = HEJ::current(p2out, p2in, flip(fhel));
const auto mj2m = HEJ::current(p2out, p2in, fhel);
const auto jgbp = HEJ::current(pg, p2in, flip(fhel));
const auto jgbm = HEJ::current(pg, p2in, fhel);
const auto j2gp = HEJ::current(p2out, pg, flip(fhel));
const auto j2gm = HEJ::current(p2out, pg, fhel);
// Dot products of these which occur again and again
COM MWmp=j_W.dot(mj2p); // And now for the Higgs ones
COM MWmm=j_W.dot(mj2m);
const auto qsum = to_tensor(q2+q3);
const auto p1o = to_tensor(p1out);
const auto p1i = to_tensor(p1in);
const auto p2o = to_tensor(p2out);
const auto p2i = to_tensor(p2in);
const auto Lmm=( (-1.)*qsum*(MWmm) + (-2.*COM{j_W.dot(pg)})*mj2m + 2.*COM{mj2m.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmm/2. ) )/q3.m2();
const auto Lmp=( (-1.)*qsum*(MWmp) + (-2.*COM{j_W.dot(pg)})*mj2p + 2.*COM{mj2p.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmp/2. ) )/q3.m2();
const auto U1mm=(COM{jgbm.dot(j_W)}*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
const auto U1mp=(COM{jgbp.dot(j_W)}*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
const auto U2mm=((-1.)*COM{j2gm.dot(j_W)}*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
const auto U2mp=((-1.)*COM{j2gp.dot(j_W)}*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
double amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
ampsq*=WProp(plbar, pl, wprop);
return ampsq/((16)*(q2.m2()*q1.m2()));
}
}
double ME_W_unob_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false, wprop);
}
double ME_W_unob_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true, wprop);
}
double ME_W_unob_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false, wprop);
}
double ME_W_unob_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true, wprop);
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (Quark - W and Uno emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (Quark - W and Uno emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
*
* Calculates j_W_{uno} ^\mu j_\mu. Ie, unordered with W emission same side.
* Handles all possible incoming states. Note this handles both forward and back-
* -ward Wuno emission. For forward, ensure p1out is the uno and W emission parton.
* @TODO: Include separate wrapper functions for forward and backward to clean up
* ME_W_unof_current in `MatrixElement.cc`.
*/
double jWuno_j(HLV pg, HLV p1out, HLV plbar, HLV pl, HLV p1in,
HLV p2out, HLV p2in, bool aqlineb,
ParticleProperties const & wprop
){
//Calculate different Helicity choices
const Helicity h = aqlineb?helicity::plus:helicity::minus;
- double ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
+ double ME2mpp = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::plus,helicity::plus, wprop);
- double ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
+ double ME2mpm = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::plus,helicity::minus, wprop);
- double ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
+ double ME2mmp = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::minus,helicity::plus, wprop);
- double ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
+ double ME2mmm = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::minus,helicity::minus, wprop);
//Helicity sum and average over initial states
return (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
}
}
double ME_Wuno_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_Wuno_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
/**
* @brief W+Jets Extremal qqx Contributions, function to handle all incoming types.
* @param pgin Incoming gluon which will split into qqx.
* @param pqout Quark of extremal qqx outgoing (W-Emission).
* @param plbar Outgoing anti-lepton momenta
* @param pl Outgoing lepton momenta
* @param pqbarout Anti-quark of extremal qqx pair. (W-Emission)
* @param pout Outgoing Particle 2 (end of FKL chain)
* @param p2in Incoming Particle 2
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W_{qqx} ^\mu j_\mu. Ie, Ex-QQX with W emission same side.
* Handles all possible incoming states. Calculated via crossing symmetry from jWuno_j.
*/
double jWqqx_j(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in, bool aqlinef,
ParticleProperties const & wprop
){
//Calculate Different Helicity Configurations.
const Helicity h = aqlinef?helicity::plus:helicity::minus;
double ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::plus,helicity::plus, wprop);
double ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::plus,helicity::minus, wprop);
double ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::minus,helicity::plus, wprop);
double ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::minus,helicity::minus, wprop);
//Helicity sum and average over initial states.
double ME2 = (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop);
}
double ME_WExqqx_qqbarQ(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop);
}
double ME_WExqqx_qbarqg(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop)
*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_WExqqx_qqbarg(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop)
*K_g(p2out,p2in)/HEJ::C_F;
}
namespace {
//Function to calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm1(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below. (Less Tensor class use?)
double t1 = (p3-pb)*(p3-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p3-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,2)/t1;
return gqqCur*(-1);
}
//Function to calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm2(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double t1 = (p2-pb)*(p2-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb,refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p2-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,1)/t1;
return gqqCur;
}
//Function to calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm3(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s23 = (p2+p3)*(p2+p3);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank1 = outer(p2+p3, metric())/s23;
Tensor<3> qqCurBlank2 = outer(pb, metric())/s23;
Tensor<1> Cur23 = HEJ::current(p2, p3,hel2);
Tensor<2> qqCur1 = qqCurBlank1.contract(Cur23,3);
Tensor<2> qqCur2 = qqCurBlank2.contract(Cur23,3);
Tensor<2> qqCur3 = qqCurBlank2.contract(Cur23,1);
Tensor<1> gqqCur = (qqCur1.contract(epsg,1)
- qqCur2.contract(epsg,2)
+ qqCur3.contract(epsg,1))*2*COM(0,1);
return gqqCur;
}
}
// no wqq emission
double ME_W_Exqqx_QQq(HLV pa, HLV pb, HLV p1, HLV p2,
HLV p3,HLV plbar, HLV pl, bool aqlinepa,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
// 2 independent helicity choices (complex conjugation related).
Tensor<1> TMmmm1 = qggm1(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm2 = qggm2(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm3 = qggm3(pb,p2,p3,minus,minus, pa);
Tensor<1> TMpmm1 = qggm1(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm2 = qggm2(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm3 = qggm3(pb,p2,p3,minus,plus, pa);
// Build the external quark line W Emmision
Tensor<1> cur1a = jW(pa,p1,plbar,pl, aqlinepa?plus:minus);
//Contract with the qqxCurrent.
COM Mmmm1 = TMmmm1.contract(cur1a,1);
COM Mmmm2 = TMmmm2.contract(cur1a,1);
COM Mmmm3 = TMmmm3.contract(cur1a,1);
COM Mpmm1 = TMpmm1.contract(cur1a,1);
COM Mpmm2 = TMpmm2.contract(cur1a,1);
COM Mpmm3 = TMpmm3.contract(cur1a,1);
//Colour factors:
COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
cm1m1=8./3.;
cm2m2=8./3.;
cm3m3=6.;
cm1m2 =-1./3.;
cm1m3 = -3.*COM(0.,1.);
cm2m3 = 3.*COM(0.,1.);
//Sqaure and sum for each helicity config:
double Mmmm = real( cm1m1*pow(abs(Mmmm1),2) + cm2m2*pow(abs(Mmmm2),2)
+ cm3m3*pow(abs(Mmmm3),2) + 2.*real(cm1m2*Mmmm1*conj(Mmmm2))
+ 2.*real(cm1m3*Mmmm1*conj(Mmmm3))
+ 2.*real(cm2m3*Mmmm2*conj(Mmmm3)) );
double Mpmm = real( cm1m1*pow(abs(Mpmm1),2) + cm2m2*pow(abs(Mpmm2),2)
+ cm3m3*pow(abs(Mpmm3),2) + 2.*real(cm1m2*Mpmm1*conj(Mpmm2))
+ 2.*real(cm1m3*Mpmm1*conj(Mpmm3))
+ 2.*real(cm2m3*Mpmm2*conj(Mpmm3)) );
// Divide by WProp
const double WPropfact = WProp(plbar, pl, wprop);
return (2*WPropfact*(Mmmm+Mpmm)/24./4.)/(pa-p1-pl-plbar).m2()/(p2+p3-pb).m2();
}
// W+Jets qqxCentral
double ME_WCenqqx_qq(HLV pa, HLV pb,HLV pl, HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
ParticleProperties const & wprop
){
init_sigma_index();
HLV pq, pqbar, p1, p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am, T4bm, T1ap, T4bp;
if(!(aqlinepa)){
T1ap = HEJ::current(p1, pa, helicity::plus);
T1am = HEJ::current(p1, pa, helicity::minus);}
else if(aqlinepa){
T1ap = HEJ::current(pa, p1, helicity::plus);
T1am = HEJ::current(pa, p1, helicity::minus);}
if(!(aqlinepb)){
T4bp = HEJ::current(p4, pb, helicity::plus);
T4bm = HEJ::current(p4, pb, helicity::minus);}
else if(aqlinepb){
T4bp = HEJ::current(pb, p4, helicity::plus);
T4bm = HEJ::current(pb, p4, helicity::minus);}
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc = MUncrossW(pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xcro = MCrossW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xsym = MSymW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is
// also the choice in qqbar helicity.
// (- - hel choice)
COM M_mmUnc = (((Xunc).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc).contract(T1am,1)).contract(T4bp,1));
COM M_mpCro = (((Xcro).contract(T1am,1)).contract(T4bp,1));
COM M_mpSym = (((Xsym).contract(T1am,1)).contract(T4bp,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc).contract(T1ap,1)).contract(T4bp,1));
COM M_ppCro = (((Xcro).contract(T1ap,1)).contract(T4bp,1));
COM M_ppSym = (((Xsym).contract(T1ap,1)).contract(T4bp,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
return amp;
}
// no wqq emission
double ME_W_Cenqqx_qq(HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
int nbelow, bool forwards, ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
if (!forwards){ //If Emission from Leg a instead, flip process.
std::swap(pa, pb);
std::reverse(partons.begin(),partons.end());
std::swap(aqlinepa, aqlinepb);
qqxmarker = !qqxmarker;
std::swap(nabove,nbelow);
}
HLV pq, pqbar, p1,p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am(0.), T1ap(0.);
if(!(aqlinepa)){
T1ap = HEJ::current(p1, pa, plus);
T1am = HEJ::current(p1, pa, minus);}
else if(aqlinepa){
T1ap = HEJ::current(pa, p1, plus);
T1am = HEJ::current(pa, p1, minus);}
Tensor <1> T4bm = jW(pb, p4, plbar, pl, aqlinepb?plus:minus);
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc_m = MUncross(pa, pq, pqbar,partons, minus, nabove);
Tensor<2> Xcro_m = MCross( pa, pq, pqbar,partons, minus, nabove);
Tensor<2> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, minus, nabove);
Tensor<2> Xunc_p = MUncross(pa, pq, pqbar,partons, plus, nabove);
Tensor<2> Xcro_p = MCross( pa, pq, pqbar,partons, plus, nabove);
Tensor<2> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, plus, nabove);
// (- - hel choice)
COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
return amp;
}
diff --git a/t/ME_data/ME_Wm.dat b/t/ME_data/ME_Wm.dat
index b1ad4a2..dfb391f 100644
--- a/t/ME_data/ME_Wm.dat
+++ b/t/ME_data/ME_Wm.dat
@@ -1,1308 +1,1308 @@
2.928617337e-06
1.647083496e-06
6.792718307e-07
8.309052795e-05
1.188563005e-09
7.418703203e-09
4.64176585e-07
4.469210588e-07
6.108953218e-06
2.937369558e-06
4.834662872e-06
8.222916804e-06
1.984842145e-06
4.724279533e-06
2.284467423e-07
6.086670249e-08
9.402694313e-07
2.134890127e-05
1.344127063e-06
6.325064033e-08
5.042105356e-06
0.0002284967536
2.294811317e-06
6.277634042e-05
1.336994329e-07
3.953756759e-06
9.658399826e-05
2.316822172e-06
2.819631851e-07
0.0001081195238
3.432930003e-06
3.368456504e-06
1.112864917e-08
3.295369143e-07
3.989490572e-10
1.948758233e-05
3.510758397e-07
3.210138517e-07
3.239664339e-06
8.604705327e-06
7.617690755e-10
1.878849261e-06
8.602382866e-09
2.905076666e-09
9.580640762e-08
2.761322408e-07
1.121735761e-05
8.305550955e-07
3.094512747e-06
9.008871903e-07
3.660037761e-07
4.545298754e-08
1.619265714e-09
0.01367608992
3.915696743e-06
2.025165585e-06
3.681426282e-09
4.179813437e-06
3.072778246e-08
8.896098616e-07
3.025899647e-08
1.817366927e-09
4.337353756e-10
1.271872323e-06
6.04489128e-06
9.893120074e-09
1.214081986e-05
2.958624219e-08
2.061238904e-09
1.930294996e-05
6.299188205e-07
2.643963602e-07
2.06132771e-06
4.342452252e-05
8.291567369e-07
2.498958153e-06
4.358403281e-06
9.241515571e-05
9.124195104e-05
0.000115945869
2.350206043e-06
7.845167967e-06
0.000638524343
4.723346316e-08
6.954346838e-09
5.357593527e-06
1.473700863e-07
1.053983829e-07
5.919949754e-07
1.175787989e-05
1.546190854e-07
0.0004236014522
3.617348006e-06
4.516789165e-05
3.955639799e-05
0.0001365301129
5.126068859e-05
1.077721123e-06
1.018738183e-08
2.31035007e-09
6.836601309e-07
5.827486967e-09
1.557627741e-09
6.04883445e-05
1.652825805e-07
2.956934131e-07
3.715517348e-06
1.403902522e-08
2.060672536e-07
0.001033124848
9.772923916e-07
2.935927964e-05
1.151232995e-05
7.794558457e-06
2.174319802e-05
4.056670359e-08
2.955320268e-08
8.030711247e-05
2.603164139e-07
0.0001681638933
2.93007525e-08
7.196892999e-05
1.920951539e-06
3.637534963e-09
7.893685649e-06
0.0001175194162
4.802829388e-10
4.835469293e-07
2.878769998e-05
1.898685256e-09
2.754462006e-08
1.403651125e-10
5.826442197e-07
2.43372504e-06
7.927550152e-06
1.020962274e-08
1.052564287e-07
7.737037652e-06
3.40057811e-06
2.714499296e-08
2.377749398e-09
0.0001490838408
9.612641064e-11
1.654718689e-05
1.570703582e-07
1.609375559e-06
6.94622935e-08
0.04337706424
9.860599192e-06
1.103383561e-05
0.01558774793
1.359111595e-06
2.104947148e-07
4.483240938e-08
2.52721137e-08
4.136183848e-07
1.75959296e-07
4.162277212e-07
3.660655583e-07
2.475053461e-06
2.142026381e-05
3.962155608e-07
5.075132791e-06
2.831915934e-08
8.260063969e-08
2.191422806e-08
7.392509744e-05
5.57319664e-10
1.264487785e-08
3.373106258e-05
1.515627561e-07
4.183562189e-07
1.162273072e-07
7.731661725e-06
1.077580922e-05
2.331616561e-07
4.316200243e-10
4.055086645e-09
1.282320598e-09
1.905109513e-09
2.799643084e-07
6.42991509e-08
3.336097219e-07
2.847456008e-07
2.87996438e-05
5.034716477e-08
0.0001697825431
7.357721468e-08
1.312946382e-07
1.494639153e-07
1.309657132e-09
1.452591322e-08
7.003335472e-07
2.392819852e-06
6.228178425e-06
2.32203151e-07
8.171995835e-08
4.289175929e-07
2.618956983e-07
7.684822085e-06
3.930002297e-10
1.235455443e-06
4.345617106e-07
1.019571055e-07
1.014379889e-05
1.198488003e-05
3.072332076e-05
2.471407108e-06
6.141337359e-06
2.217748858e-05
7.506667147e-07
4.156515027e-07
0.0003928577459
6.192388402e-07
9.058439765e-07
7.767728411e-06
1.506511059e-06
1.438338227e-08
1.393891846e-07
7.338259665e-05
1.013820716e-07
1.941576288e-05
1.098253136e-06
7.583333968e-07
1.924002116e-08
2.683360604e-06
0.0005387099327
4.513913649e-09
7.344897056e-06
4.863287455e-08
5.137627244e-07
1.340998671e-07
1.385068639e-09
0.001475702054
4.766804909e-09
2.788096601e-07
4.884961794e-07
4.389879349e-07
0.000605687279
2.260914346e-06
1.622018609e-07
8.044214563e-07
0.00287427997
8.545033973e-05
2.859858841e-07
3.519981239e-05
3.838329691e-05
6.701261402e-12
6.711423424e-09
1.221689939e-07
4.175349671e-09
1.276073867e-06
9.349187354e-06
5.09324804e-07
2.102018112e-05
2.546892483e-06
4.307286548e-07
0.0001405053615
1.532562441e-06
0.0043269957
1.798346185e-07
0.0002232869764
1.177464656e-07
1.784165987e-06
2.027475413e-09
3.749124642e-07
1.128514896e-07
4.135331177e-05
2.009064605e-08
8.773376156e-07
9.824605427e-10
1.631495832e-07
2.288308451e-05
2.035208047e-07
0.0001019083349
3.912211406e-06
5.328897283e-08
3.815188832e-10
1.388115787e-05
2.066863131e-08
9.173747359e-07
6.957092901e-09
2.903819076e-05
4.478705717e-09
6.770557613e-07
9.713885883e-08
1.99884353e-05
0.00179907318
1.236193169e-08
1.118353017e-05
1.144476494e-07
1.186451992e-06
4.716164372e-06
3.464387694e-07
7.532303889e-11
1.742260009e-06
1.345400054e-06
8.449377298e-08
8.954043245e-07
2.432439258e-06
0.000140064979
0.001713410238
1.075027005e-07
4.105324765e-05
3.733581957e-07
1.98491263e-06
0.001480293332
0.000115119501
0.002092834408
3.27549321e-06
2.185738747e-09
0.0004548048702
0.01476996084
3.934195111e-05
0.001243805117
5.430656178e-08
4.380768976e-07
4.297371205e-06
1.467603532e-08
1.820417519e-08
2.69946112e-09
8.420646529e-07
8.634740669e-10
1.18272243e-06
8.507977199e-09
7.184568904e-07
1.19256643e-06
1.766373439e-08
0.002934256247
1.283358802e-06
6.25493059e-06
8.703159136e-06
9.664582837e-06
2.675235499e-09
3.392583208e-08
0.0001071661518
5.832716966e-09
8.880940377e-08
0.0002166027692
1.248584038e-08
1.726725176e-06
2.996994836e-06
2.119813038e-07
7.412498676e-06
2.755943745e-08
9.960177922e-06
3.051353704e-05
1.129469998e-08
2.619012548e-08
5.532394281e-10
2.177115217e-10
0.000385188033
4.24781146e-05
1.005160549e-08
7.062926491e-08
0.0002453825231
6.569321281e-09
2.430064335e-06
1.72790267e-07
1.211706684e-05
9.947756935e-06
0.0001605725284
5.961264107e-08
1.637610222e-06
1.925017924e-07
1.895713903e-05
3.584070555e-09
2.424100642e-07
0.00100691167
1.00314273e-06
1.237459693e-06
2.006393648e-06
1.665644179e-05
8.731214124e-06
1.333844755e-07
7.301254085e-09
0.0005093424595
7.645270509e-08
2.020155895e-06
6.812955519e-08
1.108122048e-06
1.955477011e-05
2.003839607e-06
3.670861157e-07
1.728588066e-07
0.0007263652123
1.411883089e-05
1.57127984e-07
3.392219563e-06
1.377020396e-09
6.630519487e-07
5.22780019e-08
1.126988334e-06
0.0003754697012
1.865132299e-07
6.576821868e-06
3.329430704e-06
7.072429086e-06
2.881483653e-07
1.985401854e-07
6.928951207e-08
2.512210403e-06
6.217628744e-06
3.877664708e-07
1.707869127e-08
5.102229121e-08
4.572562413e-09
0.1094987466
7.448984403e-07
2.097362549e-06
9.378724955e-09
7.909888639e-08
1.223474695e-05
2.672139756e-10
2.314510579e-08
3.503273554e-13
1.60286267e-07
5.890515619e-07
1.487856468e-05
1.978663079e-05
2.552070699e-08
3.56553144e-08
0.000446743015
9.443500031e-10
1.362999756e-06
2.23574347e-07
1.136503384e-07
2.962418624e-08
9.303781273e-06
3.925004271e-10
1.017860324e-05
1.434319467e-08
4.477501165e-06
1.99694118e-09
4.66854211e-05
5.845087896e-06
1.50948917e-07
1.741914174e-06
2.215072337e-10
1.532586108e-06
1.571619745e-08
7.107335402e-05
7.234543614e-06
1.621610925e-06
2.564479119e-08
8.069230827e-09
2.035891754e-07
7.550566928e-08
8.659400091e-06
0.000432135357
3.337658317e-05
4.781987615e-05
2.323054819e-08
1.193074871e-05
0.0001350465294
2.75391678e-06
1.823090373e-08
2.411719368e-05
2.283448741e-09
3.766315025e-08
3.588272465e-05
6.787248878e-05
8.832303026e-10
7.778563208e-10
1.917323241e-08
5.849230322e-06
2.285366293e-06
6.146695971e-06
1.827323017e-05
3.231169839e-05
9.920988726e-07
3.779329381e-09
1.924111467e-05
7.69825658e-06
1.221181117e-06
1.642446586e-06
3.301149504e-07
0.0004602831218
3.293132243e-07
3.950606305e-07
2.805226741e-05
1.739173102e-08
3.772962441e-06
1.536977768e-06
2.43486381e-08
8.315791148e-05
6.888462138e-06
6.869231292e-08
4.672351113e-09
7.026257201e-05
7.305220983e-07
2.162227029e-05
0.0004328538034
1.653230976e-07
1.553380485e-05
5.165496905e-10
2.153688131e-10
8.592237275e-09
8.974803982e-09
7.52047119e-07
2.080291683e-07
1.653958549e-07
9.904164203e-09
6.542744984e-09
0.000531844742
1.144533099e-05
6.936072576e-05
7.619949025e-08
1.834407348e-08
8.099045298e-06
1.946623224e-09
4.915974524e-06
2.013182511e-05
8.510788604e-07
3.442227405e-07
1.19482958e-07
1.251208094e-07
2.847837381e-06
2.519188271e-06
3.662492908e-07
1.646238933e-07
5.373686421e-07
4.353389866e-08
8.727794552e-08
1.053005806e-05
1.083581872e-06
1.461382853e-06
1.220121571e-06
0.0009306878204
3.485645992e-07
2.207991124e-08
8.740223411e-08
0.0001005162172
5.658823113e-06
3.660387161e-07
5.210114851e-06
1.176487053e-09
1.16748696e-09
1.747200576e-07
4.652326082e-10
1.639389602e-06
9.199916708e-06
0.0001764717118
7.378923764e-09
2.542894196e-07
3.717694969e-06
4.778153607e-10
6.099042667e-07
7.389253413e-06
1.03018836e-08
5.091752513e-08
1.929282765e-06
2.953644841e-07
1.428238895e-08
3.65638878e-05
8.358046608e-06
1.269027303e-05
3.141972004e-10
1.192330682e-08
3.425122036e-07
1.238048594e-05
4.592097513e-07
4.907663521e-08
2.605751798e-07
1.290240408e-06
8.909975732e-08
1.171123948e-05
6.948191451e-06
0.0001320600695
4.467272334e-08
1.342984323e-06
4.697331637e-09
1.674844316e-07
1.158433939e-07
7.246890859e-08
0.0001328966329
1.455983069e-05
1.974761595e-05
9.916801323e-05
7.07587568e-07
7.363156082e-07
8.750552305e-07
2.098660427e-05
1.224863673e-06
6.973895334e-07
6.433262719e-08
1.65483483e-05
0.0009302696124
1.821542748e-07
6.800378641e-07
1.766120464e-07
8.004491159e-07
9.941703505e-05
2.067479421e-05
5.066437288e-11
1.506742077e-06
8.108300003e-05
7.610696189e-05
0.0002115668797
2.125674742e-09
2.910679195e-06
9.488563678e-06
1.847352953e-06
7.436406293e-05
1.286413705e-10
1.413078284e-06
4.724137688e-08
5.001760652e-07
1.140535015e-06
6.718486985e-11
8.472698691e-09
8.494761782e-08
2.498377868e-07
6.987065516e-06
4.081098152e-06
2.822241728e-07
3.985816309e-09
6.27810604e-05
9.793665521e-07
9.666541857e-07
2.708574051e-06
3.489974641e-07
3.782569786e-10
8.131668264e-08
6.664191677e-06
2.365456156e-08
1.096935914e-06
3.77600791e-07
2.990480133e-08
0.002754292096
4.047564209e-07
1.411960244e-07
3.69399348e-08
0.0002320662158
3.713311984e-07
3.835921825e-05
2.623885825e-05
6.474152282e-10
0.0003218958043
3.35594494e-07
5.603457411e-09
9.856525159e-08
2.793807664e-07
5.195272481e-06
2.632191143e-06
6.350175836e-07
3.302312082e-10
3.043465789e-08
2.372756421e-08
6.617008787e-06
6.330350038e-07
1.87009364e-06
2.338509504e-07
4.591079902e-10
5.619228611e-05
1.955919798e-08
3.044115064e-08
5.161003917e-08
0.0001378653082
-2.75065694e-09
+2.750698004e-09
2.289343402e-08
4.00103398e-09
2.246769087e-10
3.473356692e-06
7.311371105e-06
8.897919195e-08
1.481186501e-08
2.435159627e-12
1.189764093e-10
5.333618139e-06
4.606955115e-11
2.381338932e-10
1.249141019e-08
1.128608436e-07
1.550910634e-06
2.297659982e-12
5.846612314e-10
4.069273229e-10
4.063063753e-08
1.729353468e-08
5.155065374e-11
1.688263366e-11
2.774873693e-09
6.091545766e-08
2.871646004e-10
3.666363559e-10
8.258762033e-10
7.390763165e-09
1.487121411e-10
6.565655534e-12
3.128756673e-09
2.002340151e-09
7.274176809e-11
6.359622393e-09
7.94023499e-13
1.098589837e-06
7.438345866e-07
2.285104875e-10
4.499075704e-09
3.082821024e-09
1.878422387e-10
4.005506099e-09
2.752607668e-10
2.218745508e-09
5.947644788e-09
1.543301813e-10
4.907483457e-10
6.901483311e-10
2.081360469e-13
7.836953093e-10
2.032094869e-10
1.715010225e-08
4.069836853e-06
2.919452232e-08
1.565082965e-08
3.828958679e-08
1.256924112e-07
4.363104327e-09
2.287347968e-12
6.966805897e-09
5.995770348e-08
1.274222724e-10
1.695055516e-10
2.301236675e-11
9.445286023e-09
2.11536282e-08
1.662325435e-09
1.025048152e-09
8.700218431e-11
1.625672312e-11
4.133541356e-11
5.656983099e-11
5.288423408e-07
5.324133016e-07
2.221547731e-07
2.489666974e-09
2.994474704e-09
1.075604993e-09
7.782531146e-11
2.93245109e-11
3.399331328e-11
1.446376248e-08
2.187182264e-08
7.82008326e-09
2.402248677e-09
7.465791364e-10
1.280984452e-06
3.59815576e-08
1.251199722e-08
3.332237947e-11
7.61764668e-10
9.897388037e-11
4.653813812e-10
3.587118654e-06
4.446631546e-08
2.129602338e-09
4.744482097e-05
4.433870092e-09
2.26639483e-08
1.03614786e-09
1.103323529e-07
5.075149879e-08
8.363573856e-12
3.280957443e-10
3.14793974e-08
6.330314597e-10
3.863698763e-10
5.960231262e-10
4.380877852e-11
1.164085144e-07
2.505387626e-09
1.374083517e-10
9.172006774e-07
6.582144563e-08
1.427895895e-10
1.297546413e-09
8.107028639e-09
2.070936394e-06
6.936632176e-09
2.329619212e-07
2.684287588e-07
3.375453851e-12
4.7815963e-11
5.350876399e-09
9.177530263e-10
3.059207574e-08
1.113388771e-07
5.998531569e-10
3.911798644e-10
1.176672437e-08
9.191577823e-10
2.320742865e-05
2.036071812e-07
3.212921899e-08
7.129708873e-08
1.481190721e-09
9.133803854e-07
1.062371365e-10
1.481488075e-08
2.317757165e-06
3.853124468e-08
4.627532631e-12
2.494522948e-08
9.09456782e-10
1.680021386e-08
2.095780184e-08
6.85132284e-08
2.372661172e-09
1.177453156e-12
1.027899983e-11
5.17026087e-07
3.317347932e-07
1.0456376e-08
6.471580259e-07
1.86439548e-06
3.076357107e-11
4.028914166e-08
1.346803365e-11
3.296707439e-06
8.214485397e-11
3.289183668e-10
7.057835645e-14
2.679389626e-07
4.480005874e-10
2.542481361e-10
1.046857042e-09
4.950788629e-12
6.108069321e-10
5.254576286e-07
9.00341232e-11
3.637289725e-13
9.046831583e-09
1.769848077e-08
9.487515449e-07
1.288952596e-10
7.931101427e-08
2.511567859e-07
1.749569007e-09
1.157009753e-07
1.833616371e-06
4.450088162e-08
5.71525937e-10
1.226736711e-10
2.448689122e-10
6.262140377e-06
2.141665922e-08
4.258896529e-09
5.492356181e-09
9.566075733e-12
3.797467443e-11
2.705317002e-08
4.238826701e-07
1.78471189e-09
3.595342753e-08
3.222931873e-09
6.210042731e-05
2.257994023e-08
1.191490301e-13
1.51327299e-06
9.658649632e-11
5.082113798e-06
9.440035937e-10
8.436889449e-09
1.490034629e-10
4.468885132e-09
7.564217616e-12
1.360704417e-11
3.475352926e-10
8.989931067e-08
2.83363773e-11
3.075639475e-10
5.26550054e-14
4.262827162e-13
7.39316909e-13
5.385777656e-05
8.486065594e-10
1.27169587e-05
8.409442881e-11
6.150414316e-08
5.763923047e-10
1.482724362e-05
2.322767476e-10
4.431549181e-09
1.319773505e-08
1.965643997e-09
5.781996728e-06
7.525982729e-07
1.16731806e-09
1.712491312e-09
3.953724056e-07
3.136603556e-10
1.334448069e-10
2.291915354e-13
7.468688547e-11
2.139268308e-07
2.325345469e-09
7.657140694e-07
3.023146836e-07
1.974153842e-11
2.23954747e-13
1.575733249e-08
3.894628929e-07
2.527557586e-09
1.169459922e-07
7.445054153e-07
3.875934136e-06
6.204189806e-07
5.089515047e-10
2.563291505e-13
4.108999891e-09
8.53455392e-11
3.062860485e-14
5.084534294e-08
3.964359924e-09
9.560449902e-09
5.713347391e-07
8.602201637e-10
9.630300409e-09
5.575422411e-08
8.084446179e-10
5.256859943e-10
2.604559111e-08
5.185774593e-08
5.136812252e-07
4.457880285e-09
7.08739055e-06
2.866812716e-09
2.555644334e-07
1.149811348e-06
2.121457778e-11
6.538468276e-11
4.397488469e-10
1.683362541e-10
7.632025686e-10
3.194386969e-07
4.739399683e-08
5.174073262e-09
1.378366535e-06
0.0001075113033
0.0002104550971
1.259413081e-09
1.113788764e-09
2.836765427e-07
2.11937165e-05
3.257765239e-07
6.904195422e-08
1.133542583e-07
2.515686455e-07
2.868313218e-07
1.316229576e-07
1.143173676e-07
7.053036421e-10
1.261626874e-08
8.973217584e-13
8.228171868e-10
6.551051251e-09
5.085705662e-09
1.197384838e-11
7.765260495e-11
1.137626508e-10
9.984571963e-11
1.206612408e-07
9.836882175e-09
3.4882158e-08
4.035134021e-08
1.411934725e-08
1.030965664e-07
0.0001339279091
2.858703113e-12
3.300405884e-09
1.93837976e-10
3.214702272e-12
1.739126777e-09
4.024332505e-07
3.711093586e-08
3.509867475e-09
9.299169562e-11
4.059979323e-10
1.46124585e-05
9.47475461e-11
6.788839586e-14
2.422745433e-09
1.239867471e-07
2.586712504e-08
3.834753794e-09
4.980992443e-08
5.89070658e-10
2.437434325e-13
5.885757434e-12
3.850022999e-12
5.81521873e-10
2.871328231e-07
3.346893654e-09
4.474183737e-08
4.854163466e-09
6.835578263e-09
4.513572379e-10
1.936977318e-09
3.204323474e-07
2.834054461e-06
2.373861791e-10
8.958952159e-10
2.230057738e-07
2.935410608e-08
1.283627624e-05
1.768887803e-09
9.964539111e-08
1.803522538e-10
1.897974483e-06
1.555504868e-07
2.081802153e-11
1.130167125e-07
4.730203043e-10
4.419353612e-10
8.383693984e-08
1.081564925e-05
4.255910847e-09
4.961228471e-09
4.721951011e-12
2.605775896e-07
6.302152575e-12
7.028581764e-11
5.826657822e-09
3.775763541e-06
7.970550383e-06
1.07749277e-08
4.597131766e-10
6.829928693e-09
2.132369311e-10
3.596253106e-11
6.42042013e-09
2.577185048e-09
3.303168482e-09
6.187624832e-09
4.853422661e-08
4.830648803e-10
3.720002445e-06
1.356639118e-08
2.379953128e-09
7.264449115e-11
3.297503516e-09
3.435188638e-10
5.297642673e-08
4.382508996e-09
2.158635352e-08
5.682041402e-12
2.872850108e-07
1.684440833e-10
3.702296758e-08
4.342536978e-10
6.27913514e-08
2.403628029e-08
6.186423798e-09
1.906331939e-07
4.632141066e-10
2.75528781e-09
1.375301157e-08
2.872819245e-08
6.101281782e-07
3.268447178e-10
2.065018598e-07
2.51528016e-09
4.962377554e-11
1.429433599e-10
6.600434666e-08
2.769795382e-07
8.892370488e-07
4.79076441e-07
8.513265484e-09
2.371009137e-11
2.01515454e-08
1.194156708e-06
2.53608333e-06
3.538052912e-11
7.535772739e-11
3.279577835e-09
5.915460368e-09
3.02566311e-12
2.424965552e-08
1.736406389e-08
9.148306031e-11
2.773912766e-09
6.976061878e-08
1.281910328e-09
3.477331578e-11
7.026224303e-12
3.159072911e-08
2.312639776e-06
1.021070205e-11
6.311485446e-10
3.496457737e-09
5.006307921e-09
1.086232952e-09
6.417452399e-07
3.403922546e-08
1.040513338e-07
2.01520989e-08
2.367352254e-08
1.453582854e-08
4.519467641e-07
6.134987707e-08
1.681912432e-08
3.378933948e-10
1.788912469e-11
8.65895426e-08
1.581698325e-10
3.168928376e-08
1.269427393e-09
3.307594478e-11
6.742913842e-11
4.303451928e-10
1.049209168e-07
7.72679499e-12
5.740213241e-07
7.011858745e-07
5.337254201e-09
6.364098195e-09
6.365181097e-10
7.724336147e-10
2.368368232e-08
1.761684159e-07
3.074352635e-08
5.930269563e-06
1.348260163e-12
7.971765332e-11
1.851907721e-07
1.679589918e-10
6.480753213e-06
9.759320168e-09
3.448700092e-09
7.000832649e-11
2.060702385e-09
5.016091717e-09
1.481166798e-09
5.885608742e-10
6.092754546e-10
1.149674113e-08
8.39211811e-11
2.869929102e-13
9.120025443e-09
1.023019663e-06
1.843312565e-11
1.004917355e-06
1.437937821e-09
3.104223957e-06
3.736773683e-10
9.734184114e-11
1.231999544e-11
3.535861091e-07
7.131073672e-08
2.337280167e-10
2.824053112e-07
3.138209131e-08
5.494664929e-12
3.070138908e-12
2.011205543e-10
1.461424596e-11
3.081987435e-09
1.46093963e-10
3.186536889e-07
8.613722973e-10
4.227278814e-07
1.272955977e-08
1.38585023e-08
2.535073827e-09
5.647503952e-08
1.712170102e-10
9.68893623e-07
2.23301625e-09
2.607896714e-11
1.678599102e-13
2.258614162e-06
1.695694084e-09
1.706056963e-11
4.927634237e-12
1.42650763e-05
1.449408916e-07
4.244679669e-08
3.617422863e-10
8.345513919e-12
9.320955516e-13
8.342751299e-09
3.323822411e-09
1.436276816e-08
5.885513695e-10
3.584811231e-08
9.154788215e-09
9.640669763e-11
1.908867792e-09
2.558486e-10
5.264752426e-10
1.261283666e-09
2.53541926e-10
1.402713106e-06
9.404881776e-10
1.920817629e-07
4.404107301e-11
5.005042888e-10
1.427876454e-07
5.13977659e-10
4.146766934e-07
1.166277032e-07
7.051068181e-09
4.11849211e-09
9.268930022e-07
1.226392807e-10
1.045597589e-06
4.947490967e-07
6.160106382e-12
1.351622658e-08
2.784201053e-09
1.934479305e-08
5.745748437e-13
1.357138484e-09
4.977337067e-12
6.477023506e-11
4.234494139e-08
4.918197299e-13
3.683378601e-08
3.727148631e-09
1.040223561e-09
6.604422035e-10
3.806876455e-09
1.097067221e-08
1.705648579e-08
6.319765378e-07
1.574099802e-09
3.013722918e-08
4.082111113e-09
3.679411129e-09
8.938433543e-10
4.833608226e-09
1.526922169e-05
3.497170215e-10
9.33335016e-08
1.359274e-08
6.404241971e-11
1.5355979e-11
2.323665288e-09
4.662114978e-09
1.022032175e-10
1.680116769e-07
1.619039888e-12
1.507964967e-08
2.056097551e-09
1.984726853e-05
5.313179169e-09
1.04747219e-08
1.156429507e-07
1.201419215e-06
1.584080335e-09
2.098542248e-11
3.445376555e-09
4.164592916e-06
4.092695408e-12
9.944139167e-10
5.275904704e-09
3.628906351e-11
2.196977166e-10
1.052794797e-08
1.189717354e-08
3.768844586e-11
1.364002007e-07
3.40520559e-10
5.381516624e-08
6.729189153e-12
1.207745127e-07
4.017462281e-09
1.278030933e-08
7.005379077e-11
1.491096858e-07
7.198140186e-10
6.627930817e-08
1.215887355e-10
6.230801001e-09
0.0001012323593
5.309321999e-10
5.741270954e-10
2.429844937e-08
7.32585206e-12
3.533917464e-09
7.92607808e-10
8.291024594e-11
3.746228723e-11
5.242197859e-11
8.718274092e-08
6.175927737e-11
5.510569329e-09
2.644291342e-09
2.306855272e-08
7.426139063e-10
1.884007562e-06
3.755557449e-08
1.659745103e-08
5.494894979e-08
3.546656967e-09
2.674192614e-07
1.028231645e-09
9.289065159e-06
7.959800508e-08
1.068059708e-11
5.674230217e-07
6.214339587e-09
1.601035056e-06
1.977699126e-06
3.462717682e-07
7.373550469e-09
5.735187951e-12
diff --git a/t/ME_data/ME_Wp.dat b/t/ME_data/ME_Wp.dat
index 99d9469..74eec16 100644
--- a/t/ME_data/ME_Wp.dat
+++ b/t/ME_data/ME_Wp.dat
@@ -1,1289 +1,1289 @@
0.0001355692831
3.608650493e-05
4.306256176e-07
9.75530888e-07
7.965371874e-08
1.158718086e-06
2.661789326e-06
2.304777927e-09
2.352088884e-07
1.82969808e-06
7.663185819e-10
1.356585238e-07
3.025103442e-07
9.674650592e-08
1.085748789e-08
1.140587369e-07
9.611216759e-06
9.925878133e-06
0.000239556953
3.934982623e-07
1.135265775e-06
3.820124311e-06
9.70740276e-05
1.9898323e-08
4.670871268e-06
6.369333758e-08
5.769151072e-06
6.100650415e-07
4.177523131e-07
5.131047905e-05
2.623611134e-08
9.921070288e-07
2.818568347e-05
2.357681373e-06
1.529378671e-07
3.81204264e-09
4.089810134e-06
5.519469061e-08
3.10840544e-05
6.274958158e-07
0.0004711537356
1.064824795e-10
3.599835498e-07
0.0005693323828
5.34415283e-06
1.369492275e-07
4.610488411e-08
2.759081421e-08
0.0002011622772
8.792498535e-06
2.703844037e-06
1.87129965e-09
6.645074467e-07
0.01345885061
8.025636853e-06
1.131164601e-07
2.178426779e-07
3.510820573e-05
1.899806627e-07
1.04510564e-06
1.869058199e-06
2.768491575e-06
8.55084388e-08
7.773245411e-07
4.531353015e-06
1.601146181e-06
1.364138406e-07
5.968172458e-08
1.972845581e-07
7.86996963e-08
1.027062634e-09
1.687440419e-07
6.311963568e-10
3.39176791e-07
4.915799802e-06
1.471280469e-08
5.806849675e-07
1.107334043e-06
2.913893129e-05
2.123996588e-06
2.129863733e-07
1.356401739e-06
2.199181613e-08
9.711164556e-08
0.0001385546466
6.515093877e-07
2.535248387e-05
9.498630852e-06
7.27526662e-08
3.660448579e-08
1.057431404e-07
6.036950822e-08
1.535358058e-05
8.248221444e-07
3.681365448e-07
6.97866049e-07
0.0001134600947
0.003172309993
0.0002184096786
1.851539911e-06
2.38713427e-07
0.0002408689592
2.563539598e-07
1.000618893e-09
4.853248021e-06
1.037533394e-07
7.119446461e-06
2.561842841e-06
6.552654227e-08
2.501602017e-08
3.220051642e-08
4.669139132e-07
1.713283931e-05
1.245230123e-09
7.470021795e-06
1.121772404e-10
2.566755035e-08
0.0006902390385
1.619658523e-08
0.007842760801
3.718504779e-05
0.0007152504841
3.366806183e-05
3.590209047e-05
0.0001083332252
1.746018239e-07
2.714420914e-05
5.495029856e-05
6.534951941e-05
5.041589637e-07
3.546022348e-06
1.201771525e-08
2.615134233e-08
0.005544395296
2.66910055e-08
1.124390825e-07
2.09763416e-07
0.0004496933581
3.522832045e-06
3.609724609e-10
2.381496829e-06
0.0001723801775
8.004950126e-06
2.558924233e-05
1.661271458e-05
7.829918874e-06
8.688950059e-05
1.604189836e-05
1.563248173e-07
7.475079052e-08
3.375214615e-08
7.208326593e-08
5.524045248e-10
5.410264303e-09
9.793248952e-05
1.210164578e-06
5.403129487e-06
9.39585062e-07
9.218507242e-07
1.916781375e-05
2.177827067e-06
0.0003020837432
6.197089683e-08
7.288265251e-08
1.364167069e-05
1.038814966e-05
1.181869351e-07
5.0050635e-06
3.611153681e-08
1.221950242e-08
1.942721612e-05
3.734234992e-06
0.0001093973514
2.068010524e-06
8.857332348e-06
4.17303507e-06
6.95781085e-09
1.448248179e-06
4.551963277e-07
1.371409986e-06
8.524736448e-07
1.659108559e-09
1.245838323e-07
0.0002242593113
3.699002275e-06
2.084603102e-10
2.212804342e-08
8.425247645e-08
2.259040476e-07
0.00725976809
4.467513931e-06
3.025291035e-06
1.224916912e-08
1.156098517e-07
1.910457772e-06
6.956962769e-08
3.814058724e-08
4.664718224e-07
1.0517947e-05
8.361308985e-09
6.131166886e-07
6.896484091e-06
2.90234056e-05
7.595631192e-09
1.528964392e-08
5.225436587e-09
3.429390849e-05
2.622200403e-09
0.0005681950565
7.453650742e-07
1.807405542e-05
9.740788873e-07
2.028438792e-06
9.039473405e-08
2.480393012e-06
0.0004804465528
1.082382535e-08
2.792523478e-06
0.0002121878745
9.812357045e-05
7.82823522e-07
0.00145627451
1.943508756e-08
8.760399204e-08
6.867164243e-10
1.796199359e-07
3.651622409e-07
1.227392702e-07
6.646464152e-06
8.406802706e-05
0.001885252898
0.0001570299162
2.314662575e-06
0.00029826679
0.003229544825
2.42751347e-10
1.896000711e-06
1.190921017e-08
4.219862349e-07
7.271914713e-08
1.591930361e-06
2.437047488e-07
2.465409596e-10
3.867330857e-06
7.355187139e-09
3.882716359e-06
0.003249883727
1.636214738e-05
0.0003223729516
2.270719493e-05
1.15818731e-07
4.803334136e-06
1.938462521e-06
1.000097405e-06
7.966369258e-08
1.97712505e-05
3.697772904e-06
0.0001118251392
8.541955637e-09
3.548947328e-05
4.419488231e-08
0.0004157153041
1.269109967e-06
1.336059596e-09
0.0001810730797
5.607415682e-08
7.073695965e-08
6.226253847e-05
3.84284111e-07
6.352157121e-05
1.767499503e-08
0.002554200991
1.646506323e-08
1.960606455e-08
0.0002646754042
1.199891548e-05
2.978452263e-08
0.0001169653935
1.112477467e-08
4.229628711e-07
4.906912054e-08
1.380144673e-06
1.086541811e-05
4.639273361e-08
5.97799329e-07
7.468444459e-08
2.128065224e-08
3.477463012e-06
1.453685343e-09
3.088656052e-05
3.041112154e-05
2.604146292e-07
1.592414079e-06
1.544402089e-07
3.372312274e-06
0.0001468075351
5.574656824e-07
4.759033098e-10
0.0009306534246
5.974128962e-06
0.0002033110946
2.248575165e-06
1.737444074e-05
3.508096174e-06
2.435450701e-09
6.368877015e-06
9.703535632e-09
1.740674328e-06
3.283241228e-07
1.701391237e-08
4.192832836e-08
4.498754436e-05
3.833969139e-05
1.152598971e-09
1.277568685e-05
0.0002811181233
1.725034388e-07
1.077479686e-06
0.004189433812
1.659673188e-07
1.343332109e-07
8.327404712e-07
2.34602978e-05
5.667205166e-07
1.495413468e-09
6.283003988e-07
7.691015317e-09
1.613339674e-08
4.367646433e-08
1.328388611e-09
1.452884787e-06
2.02551508e-07
1.293118486e-06
3.674422254e-07
4.037613848e-06
2.08608188e-08
1.579955055e-06
7.231771617e-10
3.1569536e-07
0.0003659964292
0.0008948876283
5.638061947e-08
3.288523078e-06
5.432396316e-08
5.180656108e-10
5.121101877e-08
1.939412345e-08
1.979103443e-05
2.528952084e-08
0.0002538232488
3.803903342e-05
7.172183729e-06
2.75905474e-06
3.708702768e-06
4.712896021e-08
0.000108439689
5.277716937e-08
4.363976525e-08
1.209873957e-07
3.136828403e-08
1.044074244e-07
6.962817542e-08
6.454218798e-06
2.031013903e-09
4.064683174e-07
1.125182884e-06
3.916618418e-06
-5.611673539e-07
+5.611471592e-07
1.84857491e-07
2.582351752e-07
8.839694973e-05
2.850270239e-06
2.108250587e-07
1.232780491e-08
9.539060161e-08
0.006539743101
2.250901191e-05
1.893445914e-07
2.619125419e-11
3.531362021e-10
1.155012345e-07
1.892426293e-07
1.339882852e-08
5.529288677e-08
1.151877164e-05
0.00588960458
8.216093177e-05
1.557089087e-07
-9.089131199e-08
+9.089307108e-08
4.410102007e-07
4.125576496e-08
1.888774768e-07
5.726792806e-05
1.874521698e-09
9.008279257e-12
1.709010837e-05
1.679280681e-07
8.335400164e-06
1.141575757e-08
0.002680403832
1.068376798e-08
1.957144802e-06
2.244322871e-05
1.231886388e-05
4.629843546e-05
4.576136028e-06
1.980948305e-08
3.939031957e-07
2.65004331e-05
8.453580686e-06
3.871900134e-08
0.002862634137
4.135209222e-08
0.0002240265876
1.625991922e-07
1.446668925e-07
1.815685575e-06
3.823360471e-07
1.519891535e-06
7.180598153e-10
3.485367321e-05
7.35805394e-08
7.791648477e-06
7.339952614e-07
1.648796899e-05
5.701264433e-07
1.101790041e-07
2.180070124e-08
0.005048514103
1.410925493e-08
1.79659364e-05
3.537515255e-05
2.034970872e-05
1.003874866e-08
0.0007062917902
8.670420835e-05
6.826257715e-08
5.900590474e-06
1.447324088e-08
1.072029142e-05
2.460261981e-05
3.389879349e-08
3.777393397e-07
7.004751459e-07
3.633234866e-05
4.394822707e-07
7.401507266e-07
0.004667610557
4.822447961e-07
1.289578788e-08
0.0001296610888
1.387004875e-06
2.915918021e-05
1.716180445e-07
2.901178909e-08
7.505099919e-10
7.294909916e-08
5.765709885e-07
3.265707371e-08
1.700562686e-06
3.410395422e-08
7.817695837e-06
0.0001713560836
4.092858435e-08
9.55494033e-08
2.412147443e-06
3.18345347e-09
0.0006930357429
1.227979147e-07
8.804484246e-05
2.135133831e-06
4.692924121e-06
5.984249476e-08
2.487695226e-06
1.77617142e-07
6.533554512e-07
1.608219361e-06
4.482959485e-06
1.390185368e-06
2.463223456e-09
0.0001373151805
4.692426131e-10
6.223525821e-09
3.026838853e-06
5.540489257e-08
5.422757022e-08
8.715186702e-08
3.005202198e-06
4.491427703e-06
1.296306315e-08
8.033990981e-07
8.445507163e-07
1.3313495e-08
1.8082169e-06
4.811729624e-06
3.356569741e-06
3.718129986e-07
1.633178145e-06
1.145725674e-06
1.082062104e-06
6.645371526e-07
7.689911527e-06
9.446467571e-07
0.0001628133531
0.0004161446517
2.733127802e-06
4.797886191e-06
0.0001243167241
0.000448968545
4.018914394e-08
6.294006476e-05
0.000630081244
3.276911502e-07
5.253342347e-07
6.01350357e-09
6.133257751e-07
7.49793288e-08
9.707663258e-08
2.891962755e-07
0.0008858323235
5.834108655e-12
2.657205251e-05
1.937012009e-08
1.014187607e-08
2.078375068e-09
9.092481277e-09
1.171047339e-06
6.185885599e-07
7.316480548e-08
1.137906585e-05
3.92613425e-07
9.411159935e-08
2.358834566e-05
4.56232593e-07
8.746105152e-07
2.158934095e-07
7.682363523e-08
4.623384694e-08
1.822606698e-09
5.649524127e-07
8.860919132e-07
3.66172606e-06
0.005151181886
1.333234926e-08
0.01070056746
3.557241357e-07
7.955939989e-07
1.102004449e-05
9.09890981e-08
9.760274691e-07
4.02599664e-07
3.148655627e-05
3.652062557e-08
5.832473901e-08
0.0002758547993
5.979707005e-06
2.336266799e-06
1.214254103e-06
1.256993645e-05
4.380008417e-07
8.014538804e-07
0.003169517808
9.221440216e-05
4.690667282e-07
7.647898282e-07
1.418662628e-05
2.220757827e-07
3.141955147e-07
1.911552946e-07
2.452935386e-08
0.0001445565639
2.061452929e-07
1.315418046e-09
3.270353267e-05
2.302094503e-07
4.002015452e-06
1.609746316e-06
1.156175809e-08
4.490179986e-08
0.0002105147887
1.502852601e-10
0.0003938429791
2.359475691e-06
5.766404671e-05
4.012664204e-07
1.060115376e-07
0.0002141724857
2.679419004e-06
1.069242724e-05
2.484273569e-09
0.0001019811318
1.193107128e-06
1.086686152e-06
1.985318786e-07
1.172745202e-10
4.982354956e-07
2.574808657e-05
0.001250228135
2.058018342e-09
1.333953989e-07
9.818589429e-09
5.153985583e-07
0.0001083911806
1.428530264e-06
1.746216459e-05
2.151516381e-05
0.01028483602
7.323846337e-07
2.685813988e-08
4.250429571e-05
0.0004356709049
3.28939188e-07
0.003236795043
3.007501797e-08
2.793831556e-06
1.897873992e-05
1.619696108e-07
8.953551025e-05
4.80193982e-06
7.858777009e-06
7.610679109e-07
2.314240509e-05
4.939677487e-09
4.926048471e-09
8.469201572e-07
3.809176042e-05
8.401919445e-09
4.426718159e-07
2.426495876e-06
3.814034779e-08
1.746205254e-06
1.524431761e-07
7.921023056e-05
2.101316821e-07
4.839140599e-07
5.351990236e-07
4.622182851e-05
0.001900102798
1.924053736e-06
0.001430517818
-4.980374668e-08
+4.980174317e-08
0.001646472474
1.600583774e-05
2.692261472e-06
4.189059841e-06
9.98835345e-08
3.325344145e-05
1.063350006e-06
4.230188624e-06
3.988458639e-06
0.0001050664497
5.985955822e-08
2.025523764e-06
4.345105365e-07
8.910417748e-05
2.543223041e-06
0.0002129901487
8.738312346e-06
4.913864555e-09
3.35325571e-08
3.330566284e-05
1.951472299e-06
1.384757071e-05
2.148294489e-07
1.685285196e-07
1.114823446e-07
6.020716679e-05
1.280863842e-07
4.215213019e-07
3.062370712e-07
6.986571467e-07
1.615302443e-07
1.317636931e-07
0.0008030051004
2.015887208e-05
5.191764639e-06
5.499934662e-07
7.154064268e-08
0.0002157242896
9.439855684e-07
6.512827947e-08
1.039015038e-08
0.0001795882702
1.962855737e-05
1.739481498e-07
2.125722499e-08
1.743456005e-05
1.541385725e-06
3.391912637e-06
8.358930459e-06
1.415585549e-09
7.887403159e-07
1.745047173e-07
1.932691352e-07
9.827865394e-07
1.29143604e-07
9.277986603e-05
4.26401939e-06
3.754381337e-07
3.452543569e-05
8.101147971e-06
4.056176049e-09
2.2530047e-07
1.172500327e-05
4.796353912e-07
3.78522037e-09
2.390433748e-05
5.308730798e-08
4.670503157e-05
6.758043071e-05
3.147573094e-09
1.661535227e-05
5.89344189e-08
3.900862268e-05
3.504192369e-08
8.675791814e-05
0.0001335963742
1.30228891e-06
1.14637156e-06
6.394081592e-06
2.519446964e-05
0.0001367023723
4.154471425e-09
2.099781601e-06
2.146405871e-05
1.233838555e-09
0.001728366589
4.092072252e-05
0.007151322001
4.760094246e-09
0.0007525543832
1.763178471e-08
2.706266303e-07
4.168856448e-06
6.813637481e-07
6.199556862e-06
1.450329435e-09
0.001178920548
3.552454118e-07
0.001006319382
3.684787183e-07
2.167631697e-07
2.537791672e-06
7.65843049e-07
2.591360291e-07
4.604264858e-08
4.168291081e-07
7.23863419e-09
1.170104498e-06
4.393494707e-07
0.0007919172208
1.02875714e-06
6.918203902e-08
9.436152555e-08
1.591152279e-06
2.089437034e-06
1.613336951e-07
3.637029901e-07
9.604787084e-07
9.357789457e-07
2.887457829e-06
6.126125816e-06
1.015640357e-05
5.754930623e-06
0.002323638092
8.012837173e-07
5.265094095e-07
0.0002032528678
1.451002447e-07
5.476912535e-07
7.46876402e-07
1.172249495e-05
5.000091719e-06
1.5728585e-06
6.441937094e-08
2.908823657e-05
4.56908538e-05
0.0003671897961
1.009463738e-07
1.696629017e-09
1.138323555e-08
7.716768592e-07
3.529182645e-06
6.362980666e-07
1.932763096e-07
5.985715811e-07
3.856195618e-06
8.483660513e-07
7.519020755e-09
6.046415155e-06
8.331849394e-05
0.0004929171874
1.0663275e-06
6.354542182e-05
2.325133372e-05
3.330055739e-07
7.906912957e-09
3.514921427e-06
6.695060974e-06
7.791800029e-08
3.007347254e-09
1.687228718e-06
2.965956639e-06
2.551834128e-08
0.02756992212
1.902086655e-06
1.858826374e-06
1.174566239e-07
6.264996388e-07
0.0001498026534
2.510027485e-07
0.003704484331
7.864414886e-05
1.776329747e-08
0.0002011665187
4.239768745e-06
1.056552608e-06
0.0008521357496
4.768399496e-07
5.244443485e-06
0.0001425007391
4.507198467e-07
3.092477214e-07
0.007681380868
1.571782134e-07
3.341334167e-08
9.57706182e-07
1.135332842e-09
1.200841644e-07
4.167709339e-06
7.82760817e-06
1.444096635e-06
2.566989594e-05
3.703035915e-07
2.683318538e-05
2.223130286e-07
1.771561768e-05
1.154453031e-07
4.676014437e-07
1.045571408e-08
5.174119793e-07
8.734308333e-07
1.146383843e-08
8.783876977e-07
7.500376201e-10
1.921618268e-05
5.899024727e-08
7.304615354e-05
1.070796034e-08
2.432284524e-05
2.466254787e-09
1.587769497e-05
5.253475459e-08
1.186364484e-06
2.513004088e-08
2.87819867e-05
1.268569695e-09
4.036264205e-07
1.480215602e-05
0.007764494839
1.266601317e-05
2.673263662e-07
2.052245037e-07
3.026351381e-09
3.958705306e-09
5.153048282e-10
2.30485786e-10
4.00164636e-09
4.681584535e-12
1.079720947e-08
7.820293168e-07
2.049257559e-13
2.452534431e-10
6.723860105e-09
1.881353189e-09
8.769041652e-09
1.405579655e-07
8.79667941e-10
8.5986373e-09
1.838834465e-09
1.411765524e-09
3.18067218e-09
1.657130925e-09
6.291536437e-07
1.407357562e-08
2.446358306e-10
3.694611522e-06
7.168161371e-08
1.515101018e-07
5.976714028e-07
9.943558916e-11
1.770298093e-10
6.923386378e-10
3.490325032e-09
1.003706187e-11
2.246634535e-10
2.180960516e-08
1.258778496e-08
1.510668332e-09
1.860249337e-08
7.435896516e-11
8.656495714e-11
8.343559836e-07
8.213468567e-10
7.372052076e-09
1.737913173e-08
4.796532688e-10
1.006823211e-11
2.297476668e-07
2.964280788e-09
4.99821361e-06
1.217074468e-09
1.247831141e-10
2.270410706e-10
6.846887143e-08
1.126062262e-08
9.36618524e-09
8.815178295e-13
1.248917793e-09
1.287973388e-06
1.86196907e-05
5.985896019e-09
3.576308021e-09
1.396180119e-07
5.597665378e-09
1.227844729e-07
3.053732271e-10
1.281522958e-09
3.785294808e-08
8.065221219e-10
1.797290698e-08
4.870819508e-09
1.481392231e-08
6.148531319e-13
2.607634113e-08
7.503249808e-09
7.245693233e-10
3.917501637e-06
2.455952184e-09
1.635840103e-09
1.504256884e-07
4.625984077e-12
2.544655175e-10
2.326869545e-09
1.209441139e-07
9.766041156e-09
2.844111981e-07
1.20119614e-06
7.294984792e-08
2.673419175e-08
3.217599397e-11
3.051042541e-08
8.284597023e-08
2.04727224e-10
2.829564576e-09
6.317615443e-05
5.647542344e-08
3.939768099e-09
6.94521861e-09
2.688795896e-09
1.246674677e-07
1.659871523e-07
2.236153688e-05
9.480324188e-09
9.50066876e-07
7.270988017e-08
4.028863152e-07
7.587726441e-10
1.151204932e-07
3.256182092e-08
1.170362629e-12
5.503648305e-08
1.043808855e-11
1.203531525e-08
1.075547057e-10
8.107924721e-07
4.847758787e-11
7.343486683e-06
1.187542764e-11
3.503963419e-08
7.017971335e-12
2.463602443e-10
2.330168221e-10
1.513084866e-09
8.01852224e-09
3.929366966e-06
2.084034312e-07
5.951098049e-10
3.506971351e-06
5.267869558e-12
3.372542047e-06
1.801044954e-11
1.866988743e-10
2.684157842e-10
3.831407595e-10
6.969326172e-05
-1.534908483e-08
+1.534892651e-08
7.039003686e-12
1.388232647e-10
5.247749174e-08
4.682028909e-08
1.641598156e-10
2.193557269e-08
1.553195144e-09
2.996101226e-08
2.723666317e-11
1.942276011e-09
1.292931245e-11
4.121903263e-12
1.610129287e-08
1.68730662e-06
1.925468729e-09
4.331103214e-13
6.454416036e-07
6.024026922e-07
2.580228274e-10
1.225798734e-05
1.983091967e-11
4.860872791e-11
2.362341657e-08
1.373652589e-08
3.790341756e-12
1.087152141e-08
3.129247636e-07
1.364371415e-09
1.71528796e-05
3.172915696e-08
2.25476832e-08
3.114465198e-09
3.261406498e-10
2.126159374e-05
5.76634479e-08
9.666079594e-10
7.465172199e-10
1.354594056e-10
2.460812063e-08
1.087490927e-09
2.736632646e-07
5.112932147e-10
6.30608916e-07
3.16614122e-09
2.044783566e-08
6.895475483e-10
6.628499186e-11
5.398066665e-09
1.01564995e-07
5.845619561e-10
4.836690216e-08
5.100722965e-10
1.572957817e-08
1.787346401e-10
1.632619667e-08
7.017140249e-08
7.320144324e-12
1.96757824e-11
1.409100801e-12
1.238542063e-06
8.502558829e-08
3.438481899e-08
4.863624248e-07
5.154336834e-07
2.080068277e-09
3.431232831e-10
1.781888545e-08
1.544255738e-09
4.308591117e-09
2.652013597e-09
1.604288279e-11
1.019669685e-05
1.12327842e-09
2.815515241e-12
1.078735039e-10
6.773434559e-09
2.399387849e-09
1.75474283e-07
6.117989542e-11
8.592612784e-06
6.887595232e-11
1.79993583e-09
2.674398001e-10
1.267708042e-08
1.616751919e-09
2.418384582e-10
2.161734611e-07
3.594695838e-06
1.314781583e-05
3.035129573e-06
6.859937833e-12
8.125025163e-07
2.904326031e-07
8.278771155e-08
1.491230091e-12
3.232705305e-10
3.014800965e-07
3.869271268e-08
4.060064314e-09
2.722533999e-09
1.252200287e-06
8.533922342e-10
1.085340928e-13
1.013028025e-07
1.173657504e-10
1.230854977e-06
3.602985461e-07
2.340450846e-10
1.622721489e-07
1.855075276e-10
1.873406688e-13
1.477955216e-10
5.844396086e-09
5.502790891e-08
6.238719977e-09
2.281603193e-08
2.708614274e-09
8.942145684e-07
8.600069274e-08
6.474837678e-08
1.15030578e-06
4.235841586e-10
2.535746467e-12
1.472699284e-07
2.05450342e-10
2.45239943e-06
2.149773698e-09
1.302041059e-08
1.354492219e-05
1.493744046e-10
4.238601095e-09
1.975404554e-09
7.883672681e-09
3.016978986e-08
6.106630141e-07
7.432899837e-09
3.802000369e-08
1.356515625e-08
7.945785661e-11
6.146806632e-09
1.061252645e-06
3.425558912e-09
4.724377322e-09
2.264944373e-06
5.206889333e-12
1.903801687e-11
5.371716621e-09
5.143287236e-12
8.796169898e-08
2.01831747e-09
7.391246947e-07
6.374455413e-11
5.084300761e-08
4.70595383e-10
7.578932433e-10
1.02976571e-10
1.065285301e-10
2.629754276e-10
7.25335847e-09
3.644329828e-09
3.875224101e-06
1.911472842e-10
7.705710398e-09
9.497312756e-06
1.776151243e-07
7.671908505e-12
1.053728533e-07
4.52842088e-07
2.689018224e-11
2.172694222e-09
1.388552825e-10
8.063784626e-10
5.66398293e-08
1.744770454e-09
2.160682805e-11
2.110818333e-11
4.471876064e-07
1.547841859e-09
1.123968452e-08
1.220442239e-08
7.869244588e-07
1.643473453e-10
1.008574226e-08
3.353923888e-07
6.498701544e-10
1.631294723e-10
3.805045733e-07
8.427141619e-07
1.498045908e-09
1.039453208e-10
1.111541375e-08
1.747478719e-09
2.355118422e-09
1.678201094e-08
2.202184418e-10
5.649709141e-08
7.45087756e-08
3.331581749e-11
3.654180177e-11
-2.135416289e-11
+2.135443335e-11
4.454404883e-09
1.468811182e-11
6.57114871e-09
5.984556178e-10
1.45099015e-07
1.39419768e-10
1.042111653e-11
3.257323655e-11
9.125783554e-09
1.15776364e-08
1.347035226e-10
1.772518894e-07
2.065966085e-07
1.18335761e-10
4.340575369e-11
1.724054892e-07
1.017507229e-09
4.448105514e-08
9.297852627e-10
1.457953015e-08
1.344680446e-06
3.942919772e-10
6.431630422e-10
2.598292431e-09
9.723217413e-11
5.887976816e-07
1.250935985e-08
1.91892714e-08
1.266147267e-09
2.358514388e-08
9.312371762e-07
4.777005653e-09
1.101297338e-05
4.464552646e-10
7.607715154e-08
-4.820145421e-10
+4.82066931e-10
2.264473432e-08
9.048298519e-09
8.615187911e-07
1.40027638e-06
2.971161233e-07
1.70303243e-09
9.085713716e-08
7.676704887e-09
3.495574146e-07
5.348513441e-10
8.550322518e-10
5.35951439e-10
3.270519121e-08
2.199299151e-11
3.481764536e-07
6.965854851e-10
1.079164489e-09
2.017281663e-11
1.225752127e-08
1.755199899e-09
3.79819579e-08
3.251519125e-06
1.339537878e-08
1.067186532e-10
2.907350796e-09
2.121014228e-09
1.792310313e-10
5.485037209e-10
2.677594237e-10
3.400494514e-09
2.722261269e-07
1.221064105e-07
4.918083741e-10
1.448712678e-08
1.186087559e-10
2.383805467e-09
6.366306819e-10
2.884744386e-08
3.233589437e-09
0.0001662944389
2.948950879e-10
2.818085744e-08
2.149521422e-09
2.878511074e-11
1.429955578e-11
2.151447291e-11
4.991716189e-07
7.608410036e-11
4.322111296e-08
8.515575318e-10
4.52193765e-09
4.65887361e-08
4.872378255e-10
1.225930736e-08
9.732262711e-09
5.360870618e-08