diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 2314ca5..637c3a0 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1761 +1,1743 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <unordered_map>
#include <utility>
#include "CLHEP/Vector/LorentzVector.h"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/Constants.hh"
#include "HEJ/currents.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/utility.hh"
namespace HEJ{
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const {
const double lambda = param_.regulator_lambda;
const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
// use alpha_s(sqrt(q_j*lambda)), evolved to mur
return (
1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(
Event const & event
) const {
return tree(event)*virtual_corrections(event);
}
Weights MatrixElement::tree(
Event const & event
) const {
return tree_param(event)*tree_kin(event);
}
Weights MatrixElement::tree_param(
Event const & event
) const {
if(! is_HEJ(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
Weights MatrixElement::virtual_corrections(
Event const & event
) const {
if(! is_HEJ(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = virtual_corrections(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = virtual_corrections(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
double mur,
Particle const & WBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
bool wc = true;
bool wqq = false;
// With extremal qqx or unordered gluon outside the extremal
// partons then it is not part of the FKL ladder and does not
// contribute to the virtual corrections. W emitted from the
// most backward leg must be taken into account in t-channel
if (event.type() == event_type::FKL) {
if (in[0].type != partons[0].type ){
q -= WBoson.p;
wc = false;
}
}
else if (event.type() == event_type::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
wc = false;
}
}
else if (event.type() == event_type::qqxexb) {
q -= partons[1].p;
++first_idx;
if (abs(partons[0].type) != abs(partons[1].type)){
q -= WBoson.p;
wc = false;
}
}
if(event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(abs(backquark->type) != abs((backquark+1)->type)) {
wqq=true;
wc=false;
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -=partons[j+1].p;
} // End Loop one
if (last_idx != first_idx_qqx) q -= partons[last_idx+1].p;
if (wqq) q -= WBoson.p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -= partons[j+1].p;
}
if (wc) q -= WBoson.p;
assert(
nearby(q, -1*pb, norm)
|| is_AWZH_boson(partons.back().type)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
double MatrixElement::virtual_corrections(
Event const & event,
double mur
) const{
auto const & in = event.incoming();
auto const & out = event.outgoing();
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
const auto AWZH_boson = std::find_if(
begin(out), end(out),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson != end(out) && abs(AWZH_boson->type) == pid::Wp){
return virtual_corrections_W(event, mur, *AWZH_boson);
}
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
size_t first_idx = 0;
size_t last_idx = out.size() - 1;
// if there is a Higgs boson, extremal qqx or unordered gluon
// outside the extremal partons then it is not part of the FKL
// ladder and does not contribute to the virtual corrections
if((out.front().type == pid::Higgs)
|| event.type() == event_type::unob
|| event.type() == event_type::qqxexb){
q -= out[1].p;
++first_idx;
}
if((out.back().type == pid::Higgs)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(out) + 1, end(out) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
);
if(backquark == end(out) || (backquark+1)->type==pid::gluon) return 0;
last_idx = std::distance(begin(out), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
if (last_idx != first_idx_qqx) q -= out[last_idx+1].p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| out.back().type == pid::Higgs
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
} // namespace HEJ
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(
CLHEP::HepLorentzVector qav,
CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1,
CLHEP::HepLorentzVector p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
return Cls-4./(kperp*kperp);
}
}
//! Lipatov vertex
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip,
CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(
CLHEP::HepLorentzVector qav,
CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pa,
CLHEP::HepLorentzVector pb,
CLHEP::HepLorentzVector p1,
CLHEP::HepLorentzVector p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
double temp=Cls-4./(kperp*kperp);
return temp;
}
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
if (aptype==21&&bptype==21) {
return jM2gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2qg(pn,pb,p1,pa);
else
return jM2qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return jM2qg(p1,pa,pn,pb);
else
return jM2qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2qQ(pn,pb,p1,pa);
else
return jM2qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return jM2qQbar(p1,pa,pn,pb);
else
return jM2qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_W_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// We know it cannot be gg incoming.
assert(!(aptype==21 && bptype==21));
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jMWqg(pn,plbar,pl,pb,p1,pa);
else
return jMWqbarg(pn,plbar,pl,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return jMWqg(p1,plbar,pl,pa,pn,pb);
else
return jMWqbarg(p1,plbar,pl,pa,pn,pb);
}
else { // they are both quark
if (wc==true){ // emission off b, (first argument pbout)
if (bptype>0) {
if (aptype>0)
return jMWqQ(pn,plbar,pl,pb,p1,pa);
else
return jMWqQbar(pn,plbar,pl,pb,p1,pa);
}
else {
if (aptype>0)
return jMWqbarQ(pn,plbar,pl,pb,p1,pa);
else
return jMWqbarQbar(pn,plbar,pl,pb,p1,pa);
}
}
else{ // emission off a, (first argument paout)
if (aptype > 0) {
if (bptype > 0)
return jMWqQ(p1,plbar,pl,pa,pn,pb);
else
return jMWqQbar(p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0)
return jMWqbarQ(p1,plbar,pl,pa,pn,pb);
else
return jMWqbarQbar(p1,plbar,pl,pa,pn,pb);
}
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*/
double ME_W_unob_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// we know they are not both gluons
if (bptype == 21 && aptype != 21) { // b gluon => W emission off a
if (aptype > 0)
return jM2Wunogqg(pg,p1,plbar,pl,pa,pn,pb);
else
return jM2Wunogqbarg(pg,p1,plbar,pl,pa,pn,pb);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return junobMWqQg(pn,plbar,pl,pb,p1,pa,pg);
else
return junobMWqQbarg(pn,plbar,pl,pb,p1,pa,pg);
}
else{
if (aptype>0)
return junobMWqbarQg(pn,plbar,pl,pb,p1,pa,pg);
else
return junobMWqbarQbarg(pn,plbar,pl,pb,p1,pa,pg);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return jM2WunogqQ(pg,p1,plbar,pl,pa,pn,pb);
else //qqbar
return jM2WunogqQbar(pg,p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return jM2WunogqbarQ(pg,p1,plbar,pl,pa,pn,pb);
else //qbarqbar
return jM2WunogqbarQbar(pg,p1,plbar,pl,pa,pn,pb);
}
}
}
}
/** Matrix element squared for uno forward tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unof Tree-Level Current-Current Scattering
*/
double ME_W_unof_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// we know they are not both gluons
if (aptype==21 && bptype!=21) {//a gluon => W emission off b
if (bptype > 0)
return jM2Wunogqg(pg, pn,plbar, pl, pb, p1, pa);
else
return jM2Wunogqbarg(pg, pn,plbar, pl, pb, p1, pa);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return jM2WunogqQ(pg,pn,plbar,pl,pb,p1,pa);
else
return jM2WunogqQbar(pg,pn,plbar,pl,pb,p1,pa);
}
else{
if (aptype>0)
return jM2WunogqbarQ(pg,pn,plbar,pl,pb,p1,pa);
else
return jM2WunogqbarQbar(pg,pn,plbar,pl,pb,p1,pa);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return junofMWgqQ(pg,pn,pb,p1,plbar,pl,pa);
else //qqbar
return junofMWgqQbar(pg,pn,pb,p1,plbar,pl,pa);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return junofMWgqbarQ(pg,pn,pb,p1,plbar,pl,pa);
else //qbarqbar
return junofMWgqbarQbar(pg,pn,pb,p1,plbar,pl,pa);
}
}
}
}
/** \brief Matrix element squared for backward qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxb Tree-Level Current-Current Scattering
*/
double ME_W_qqxb_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFbackward;
if (pqbar.rapidity() > pq.rapidity()){
swapQuarkAntiquark=true;
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pq.minus())+(pq.minus())/pa.minus())+1./3.)*3./4.;
}
else{
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pqbar.minus())+(pqbar.minus())/pa.minus())+1./3.)*3./4.;
}
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swapQuarkAntiquark){
return jM2Wggtoqqbarg(pa, pqbar, plbar, pl, pq, pn,pb)*CFbackward;}
else {
return jM2Wggtoqbarqg(pa, pq, plbar, pl, pqbar, pn,pb)*CFbackward;}
}
else if (aptype==21&&bptype!=21 ) {//a gluon => W emission off b leg or qqx
if (wc!=1){ // W Emitted from backwards qqx
if (swapQuarkAntiquark){
return jM2WgQtoqqbarQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
else{
return jM2WgQtoqbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
}
else { // W Must be emitted from forwards leg.
if(bptype > 0){
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pb, pa, pn, pqbar, pq, plbar, pl, false)*CFbackward;}
else{
return jM2WgqtoQQqW(pb, pa, pn, pq, pqbar, plbar, pl, false)*CFbackward;}
} else {
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pb, pa, pn, pqbar, pq, plbar, pl, true)*CFbackward;}
else{
return jM2WgqtoQQqW(pb, pa, pn, pq, pqbar, plbar, pl, true)*CFbackward;}
}
}
}
else{
throw std::logic_error("Incompatible incoming particle types with qqxb");
}
}
/* \brief Matrix element squared for forward qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param p1 Final state 1 Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxf Tree-Level Current-Current Scattering
*/
double ME_W_qqxf_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFforward;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pq.plus())+(pq.plus())/pb.plus())+1./3.)*3./4.;
}
else{
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pqbar.plus())+(pqbar.plus())/pb.plus())+1./3.)*3./4.;
}
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swapQuarkAntiquark){
return jM2Wggtoqqbarg(pb, pqbar, plbar, pl, pq, p1,pa)*CFforward;}
else {
return jM2Wggtoqbarqg(pb, pq, plbar, pl, pqbar, p1,pa)*CFforward;}
}
else if (bptype==21&&aptype!=21) {// b gluon => W emission off a or qqx
if (wc==1){ // W Emitted from forwards qqx
if (swapQuarkAntiquark){
return jM2WgQtoqbarqQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
else {
return jM2WgQtoqqbarQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
}
// W Must be emitted from backwards leg.
if (aptype > 0){
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pa,pb, p1, pqbar, pq, plbar, pl, false)*CFforward;}
else{
return jM2WgqtoQQqW(pa,pb, p1, pq, pqbar, plbar, pl, false)*CFforward;}
} else
{
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pa,pb, p1, pqbar, pq, plbar, pl, true)*CFforward;}
else{
return jM2WgqtoQQqW(pa,pb, p1, pq, pqbar, plbar, pl, true)*CFforward;}
}
}
else{
throw std::logic_error("Incompatible incoming particle types with qqxf");
}
}
/* \brief Matrix element squared for central qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqxpair
* @param nbelow Number of gluons emitted after central qqxpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum\
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wqq Boolean. True siginfies W boson is emitted from Central qqx
* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
*/
double ME_W_qqxmid_current(
int aptype, int bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<HLV> partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
bool swapQuarkAntiquark=false;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
}
double CFforward = (0.5*(3.-1./3.)*( pb.plus()/(partons[partons.size()-1].plus())
+ (partons[partons.size()-1].plus())/pb.plus() )+1./3.)*3./4.;
double CFbackward = (0.5*(3.-1./3.)*( pa.minus()/(partons[0].minus())
+ (partons[0].minus())/pa.minus() )+1./3.)*3./4.;
double wt=1.;
if (aptype==21) wt*=CFbackward;
if (bptype==21) wt*=CFforward;
if (aptype <=0 && bptype <=0){ // Both External AntiQuark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons,true,true, swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true,true, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true,true, swapQuarkAntiquark, nabove, nbelow, false);
}
} // end both antiquark
else if (aptype<=0){ // a is antiquark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons, false, true, swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons,false,true, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, true, swapQuarkAntiquark, nabove, nbelow, false);
}
} // end a is antiquark
else if (bptype<=0){ // b is antiquark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow, false);
}
} //end b is antiquark
else{ //Both Quark or gluon
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl, plbar, partons, false, false, swapQuarkAntiquark, nabove);}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow, false);
}
}
}
/** \brief Matrix element squared for tree-level current-current scattering with Higgs
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
*/
double ME_Higgs_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype==21) // gg initial state
return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered forward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param punof Unordered Particle Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered forward emission
*/
double ME_Higgs_current_unof(
int aptype, int bptype,
CLHEP::HepLorentzVector const & punof,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (aptype>0)
return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param punob Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*/
double ME_Higgs_current_unob(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & punob,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (bptype==21&&aptype!=21) {
if (aptype > 0)
return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (bptype>0)
return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
void validate(HEJ::MatrixElementConfig const & config) {
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(!config.Higgs_coupling.use_impact_factors) {
throw std::invalid_argument{
"Invalid Higgs coupling settings.\n"
"HEJ without QCDloop support can only use impact factors.\n"
"Set use_impact_factors to true or recompile HEJ.\n"
};
}
#endif
if(config.Higgs_coupling.use_impact_factors
&& config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) {
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
} // namespace anonymous
namespace HEJ{
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
double MatrixElement::tree_kin(
Event const & ev
) const {
if(! is_HEJ(ev.type())) return 0.;
auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing()))
return tree_kin_jets(ev);
switch(AWZH_boson->type){
case pid::Higgs:
return tree_kin_Higgs(ev);
case pid::Wp:
case pid::Wm:
return tree_kin_W(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
namespace{
constexpr int extremal_jet_idx = 1;
constexpr int no_extremal_jet_idx = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == extremal_jet_idx;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
double lambda
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
Event const & ev
) const{
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
// TODO: avoid reclustering
fastjet::ClusterSequence cs(to_PseudoJet(out_partons), ev.jet_def());
const auto jets = sorted_by_rapidity(cs.inclusive_jets(ev.min_jet_pt()));
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission or qqx
if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){
assert(jets.size() >= 3);
++most_backward;
}
else if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
assert(jets.size() >= 3);
--most_forward;
}
const auto extremal_jet_indices = cs.particle_jet_indices(
{*most_backward, *most_forward}
);
assert(extremal_jet_indices.size() == out_partons.size());
for(size_t i = 0; i < out_partons.size(); ++i){
assert(HEJ::is_parton(out_partons[i]));
const int idx = (extremal_jet_indices[i]>=0)?
extremal_jet_idx:
no_extremal_jet_idx;
out_partons[i].p.set_user_index(idx);
}
return out_partons;
}
double MatrixElement::tree_kin_jets(
Event const & ev
) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
if(is_uno(ev.type())){
throw not_implemented("unordered emission not implemented for pure jets");
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
pa - p1, pa, pb, p1, pn,
param_.regulator_lambda
);
}
namespace{
double tree_kin_W_FKL(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
double tree_kin_W_unob(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
auto pg = to_HepLorentzVector(partons[0]);
auto p1 = to_HepLorentzVector(partons[1]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - p1 -pg;
if (aptype!=partons[1].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_unob_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_unof(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 2]);
auto pg = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_unof_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxb(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
HLV pq,pqbar;
if(is_quark(partons[0])){
pq = to_HepLorentzVector(partons[0]);
pqbar = to_HepLorentzVector(partons[1]);
}
else{
pq = to_HepLorentzVector(partons[1]);
pqbar = to_HepLorentzVector(partons[0]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - pq - pqbar;
if (partons[1].type!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_qqxb_current(
aptype, bptype, pa, pb,
pq, pqbar, pn, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxf(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
HLV pq,pqbar;
if(is_quark(partons[partons.size() - 1])){
pq = to_HepLorentzVector(partons[partons.size() - 1]);
pqbar = to_HepLorentzVector(partons[partons.size() - 2]);
}
else{
pq = to_HepLorentzVector(partons[partons.size() - 2]);
pqbar = to_HepLorentzVector(partons[partons.size() - 1]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_qqxf_current(
aptype, bptype, pa, pb,
pq, pqbar, p1, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxmid(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
bool wc, wqq;
if (backmidquark->type == -(backmidquark+1)->type){ // Central qqx does not emit
wqq=false;
if (aptype==partons[0].type) {
wc = true;
}
else{
wc = false;
q0-=pl+plbar;
}
}
else{
wqq = true;
wc = false;
qqxt-=pl+plbar;
}
auto begin_ladder = begin(partons) + 1;
auto end_ladder_1 = (backmidquark);
auto begin_ladder_2 = (backmidquark+2);
auto end_ladder = end(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
int nabove = std::distance(begin_ladder, backmidquark);
int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqxt, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
} // namespace anonymous
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
auto const & decays(ev.decays());
HLV plbar, pl;
for (auto& x: decays) {
if (x.second.at(0).type < 0){
plbar = to_HepLorentzVector(x.second.at(0));
pl = to_HepLorentzVector(x.second.at(1));
}
else{
pl = to_HepLorentzVector(x.second.at(0));
plbar = to_HepLorentzVector(x.second.at(1));
}
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
if(ev.type() == unordered_backward){
return tree_kin_W_unob(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == unordered_forward){
return tree_kin_W_unof(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == extremal_qqxb){
return tree_kin_W_qqxb(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == extremal_qqxf){
return tree_kin_W_qqxf(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == central_qqx){
return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
double MatrixElement::tree_kin_Higgs(
Event const & ev
) const {
if(is_uno(ev.type())){
return tree_kin_Higgs_between(ev);
}
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev);
}
return tree_kin_Higgs_between(ev);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
#ifdef HEJ_BUILD_WITH_QCDLOOP
// TODO: code duplication with currents.cc
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(type == ParticleID::gluon) return K_g(pout, pin);
return C_F;
}
#endif
// Colour factor in strict MRK limit
double K_MRK(ParticleID type) {
return (type == ParticleID::gluon)?C_A:C_F;
}
}
double MatrixElement::MH2_forwardH(
CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
ParticleID type2,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
// gluon case
#ifdef HEJ_BUILD_WITH_QCDLOOP
if(!param_.Higgs_coupling.use_impact_factors){
return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*MH2gq_outsideH(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(
Event const & ev
) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(outgoing[1].type != pid::gluon) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
const auto q0 = pa - p1 - pH;
const double t1 = q0.m2();
const double t2 = (pn - pb).m2();
return MH2_forwardH(
p1, pa, incoming[1].type, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_last(
Event const & ev
) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.back().type == pid::Higgs);
if(outgoing[outgoing.size()-2].type != pid::gluon) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
const double t1 = q0.m2();
const double t2 = (pn + pH - pb).m2();
return MH2_forwardH(
pn, pb, incoming[0].type, p1, pa, pH,
t2, t1
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_between(
Event const & ev
) const {
using namespace event_type;
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[(ev.type() == unob)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - ((ev.type() == unof)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
(ev.type() == unob)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
(ev.type() == unof)
|| partons.front().type != pid::gluon
))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor;
if(ev.type() == unob){
current_factor = C_A*C_A/2.*ME_Higgs_current_unob( // 1/2 = "K_uno"
incoming[0].type, incoming[1].type,
pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
current_factor = C_A*C_A/2.*ME_Higgs_current_unof( // 1/2 = "K_uno"
incoming[0].type, incoming[1].type,
to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn,
param_.regulator_lambda
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn,
param_.regulator_lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
- double MatrixElement::tree_param_partons(
- double alpha_s, double mur,
- std::vector<Particle> const & partons
- ) const{
- const double gs2 = 4.*M_PI*alpha_s;
- double wt = std::pow(gs2, partons.size());
- if(param_.log_correction){
- // use alpha_s(q_perp), evolved to mur
- assert(partons.size() >= 2);
- for(size_t i = 1; i < partons.size()-1; ++i){
- wt *= 1. + alpha_s/(2.*M_PI)*beta0*log(mur/partons[i].p.perp());
- }
- }
- return wt;
- }
-
namespace {
double get_AWZH_coupling(Event const & ev, double alpha_s) {
const auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](auto const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing())) return 1.;
switch(AWZH_boson->type){
case pid::Higgs:
return alpha_s*alpha_s;
case pid::Wp:
case pid::Wm:
return gw*gw*gw*gw/4.;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
}
double MatrixElement::tree_param(
Event const & ev,
double mur
) const{
assert(is_HEJ(ev.type()));
- const auto & out = ev.outgoing();
+ const auto partons = filter_partons(ev.outgoing());
const double alpha_s = alpha_s_(mur);
- const double AWZH_coupling = get_AWZH_coupling(ev, alpha_s);
- if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){
- return AWZH_coupling*4.*M_PI*alpha_s*tree_param_partons(
- alpha_s, mur, filter_partons({begin(out) + 1, end(out)})
- );
- }
- if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
- return AWZH_coupling*4.*M_PI*alpha_s*tree_param_partons(
- alpha_s, mur, filter_partons({begin(out), end(out) - 1})
- );
+ const double gs2 = 4.*M_PI*alpha_s;
+ double res = std::pow(gs2, partons.size());
+ if(param_.log_correction){
+ // use alpha_s(q_perp), evolved to mur
+ assert(partons.size() >= 2);
+ for(size_t i = 1; i < partons.size()-1; ++i){
+ res *= 1. + alpha_s/(2.*M_PI)*beta0*log(mur/partons[i].p.perp());
+ }
}
- return AWZH_coupling*tree_param_partons(alpha_s, mur, filter_partons(out));
+ return get_AWZH_coupling(ev, alpha_s)*res;
}
} // namespace HEJ
diff --git a/t/ME_data/ME_Wm.dat b/t/ME_data/ME_Wm.dat
index 5c37191..95149c7 100644
--- a/t/ME_data/ME_Wm.dat
+++ b/t/ME_data/ME_Wm.dat
@@ -1,1308 +1,1308 @@
2.928675692e-06
1.647116461e-06
6.792853649e-07
8.309213874e-05
1.188586835e-09
7.418850985e-09
4.641858351e-07
4.469299635e-07
6.109075964e-06
2.937428029e-06
4.834759305e-06
8.223080997e-06
1.984882889e-06
4.724373663e-06
2.284512955e-07
6.086791509e-08
9.402879575e-07
2.134932663e-05
1.344153838e-06
6.325190057e-08
5.042205822e-06
0.0002285013435
2.294857015e-06
6.277759204e-05
1.337020968e-07
3.953835541e-06
9.658506379e-05
2.316868337e-06
2.819688016e-07
0.0001081216619
3.432998402e-06
3.368523741e-06
1.11288709e-08
3.29543466e-07
3.989570061e-10
1.948796723e-05
1.282013086e-10
3.210202477e-07
3.239728865e-06
8.604875852e-06
7.617842534e-10
1.878886612e-06
8.602552313e-09
2.905134545e-09
9.580831652e-08
2.761377426e-07
1.121758088e-05
8.305715929e-07
3.094575366e-06
9.009051399e-07
3.660116555e-07
4.545389317e-08
1.619299497e-09
0.01367636241
3.915774786e-06
2.025205952e-06
3.6814996e-09
4.179895596e-06
3.07283947e-08
8.89627591e-07
-2.048586989e-08
+1.381537057737695e-08
1.817403251e-09
4.337439968e-10
1.271896091e-06
6.045011555e-06
9.89332115e-09
1.214106087e-05
2.958683161e-08
2.061281825e-09
1.93033348e-05
6.299313702e-07
2.644016282e-07
2.061367847e-06
4.34253878e-05
1.164199395e-06
2.499007952e-06
4.358489891e-06
9.241699824e-05
9.12437674e-05
0.0001159481718
2.350252878e-06
7.845324212e-06
0.000638536879
4.723440426e-08
6.954486577e-09
5.357700283e-06
1.473730225e-07
1.054004829e-07
5.920067571e-07
1.175811475e-05
-5.246498884e-08
+3.538161996944155e-08
0.0004236100189
3.617420289e-06
4.516879429e-05
3.955718614e-05
0.0001365328248
5.126170478e-05
2.005036424e-07
1.018758481e-08
-9.831902939e-09
+6.630491325177325e-09
6.836738619e-07
5.827603078e-09
1.557658774e-09
6.048954972e-05
1.652858793e-07
2.956993046e-07
3.715591231e-06
1.403931806e-08
2.060713606e-07
0.001033146046
9.773118638e-07
2.935985027e-05
1.151256054e-05
7.794713809e-06
2.17436307e-05
4.056751186e-08
2.955377882e-08
8.030872004e-05
2.603215792e-07
0.0001681672442
2.930133644e-08
7.197036394e-05
1.920989813e-06
3.637607439e-09
7.893842383e-06
0.0001175217074
4.802925073e-10
4.83556567e-07
2.878827354e-05
1.898723086e-09
2.754516812e-08
1.403679092e-10
5.826558306e-07
2.433773671e-06
7.927708105e-06
1.020983626e-08
1.052585256e-07
7.737192672e-06
3.400645867e-06
2.714553382e-08
6.774637317e-09
0.0001490868071
9.612832588e-11
1.654742249e-05
1.57073488e-07
1.609407803e-06
2.782539582e-09
0.04337790658
9.860795665e-06
1.103405546e-05
0.01558805851
1.359140493e-06
2.104988229e-07
4.483328013e-08
2.527261731e-08
1.967859428e-09
1.75962763e-07
4.16236019e-07
3.66072855e-07
2.475097776e-06
2.142069091e-05
3.962234552e-07
5.0752339e-06
2.831972358e-08
8.260228591e-08
2.191466894e-08
7.392657036e-05
1.244008807e-09
1.264513056e-08
3.373173459e-05
1.515657759e-07
4.183645033e-07
1.162296215e-07
7.731815775e-06
1.077601875e-05
2.331663061e-07
4.316286429e-10
4.05516744e-09
1.282346148e-09
-2.825013126e-09
+1.905147471814478e-09
2.799718364e-07
6.430044079e-08
3.336162544e-07
2.847516226e-07
2.8800296e-05
5.034816791e-08
0.0001697859258
7.357868068e-08
-1.946917349e-07
+1.312972541815244e-07
1.494668933e-07
1.309682993e-09
1.452620265e-08
7.003475182e-07
2.392867531e-06
6.228302547e-06
2.32207831e-07
8.172158542e-08
4.289261389e-07
2.619009165e-07
7.684975193e-06
3.149447044e-10
1.235480059e-06
4.345702611e-07
1.019591398e-07
1.014400081e-05
1.198511446e-05
3.072393212e-05
2.47145635e-06
6.141459748e-06
2.217791841e-05
7.506816715e-07
4.156597855e-07
0.0003928655193
6.192512495e-07
9.058620171e-07
7.767884384e-06
1.506541078e-06
1.438366886e-08
2.856106448e-08
-4.815309756e-05
+3.247374364864003e-05
1.013840918e-07
1.941618603e-05
1.098271271e-06
7.583485063e-07
1.924040495e-08
2.683414048e-06
0.0005387206665
4.514003014e-09
7.345041607e-06
4.863384343e-08
5.137710227e-07
1.341025387e-07
1.385096075e-09
0.001475731456
4.766899706e-09
2.788152153e-07
4.885059085e-07
4.389966638e-07
0.0006056993831
2.260959244e-06
1.622050932e-07
8.044376305e-07
0.002874337098
8.545204208e-05
2.859915875e-07
3.520051428e-05
3.83840617e-05
-6.693256832e-11
+4.513834365350505e-11
6.711555934e-09
1.221714284e-07
4.175432874e-09
1.276099272e-06
9.349373741e-06
5.093349519e-07
2.102059997e-05
2.546943624e-06
4.307372366e-07
0.0001405081607
1.532593492e-06
0.004327081818
1.798381994e-07
0.0002232914245
1.177488117e-07
1.784201554e-06
2.027515818e-09
3.74919864e-07
1.128537898e-07
4.135413677e-05
1.320654261e-08
-1.300969977e-06
+8.773550961839615e-07
9.824802351e-10
-2.419281993e-07
+1.631528338714485e-07
2.288354101e-05
2.035248598e-07
0.0001019103655
3.912289372e-06
5.329003459e-08
3.815264848e-10
1.388142926e-05
-3.064871302e-08
+2.066904312077937e-08
9.173930152e-07
6.957230348e-09
2.903877008e-05
4.478794954e-09
6.770693006e-07
9.714079428e-08
1.998883986e-05
0.001799109057
1.236217788e-08
1.118375624e-05
1.144499288e-07
1.1864756e-06
4.716258378e-06
3.464456734e-07
7.532453996e-11
1.742294659e-06
4.696151226e-09
8.449545648e-08
8.954221651e-07
2.432487759e-06
0.0001400677698
0.001713444386
1.075048353e-07
4.105406541e-05
3.733656311e-07
1.984952182e-06
0.001480322688
0.0001151217974
0.00209287627
3.27555809e-06
2.185782336e-09
0.0004548139237
0.01477027961
3.934273672e-05
0.001243829892
5.430764382e-08
4.380856265e-07
4.297458002e-06
1.467632945e-08
1.820453663e-08
2.699514835e-09
8.420814307e-07
8.634913084e-10
1.182745993e-06
8.508146264e-09
7.184712053e-07
1.192590185e-06
1.766408443e-08
0.002934313619
1.283384372e-06
6.255055172e-06
8.703331507e-06
9.664774229e-06
2.675288802e-09
3.392650797e-08
0.0001071682868
5.832833131e-09
8.881117172e-08
0.0002166070842
1.248608916e-08
1.726759584e-06
2.99705455e-06
2.119854874e-07
7.412646571e-06
2.755998704e-08
9.960376424e-06
3.051416349e-05
1.129492506e-08
-3.883632292e-08
+2.619064730671114e-08
5.532504512e-10
2.177158595e-10
0.0003851956976
4.247896087e-05
1.005180535e-08
7.063067217e-08
0.0002453874122
6.569452237e-09
2.430119913e-06
1.727937098e-07
1.211730869e-05
9.94795514e-06
0.0001605757277
5.961378248e-08
1.63764285e-06
1.925056588e-07
1.895751671e-05
3.584139499e-09
2.424148932e-07
0.001006931825
6.133075998e-07
1.237484464e-06
2.006432918e-06
1.665677334e-05
8.731387573e-06
1.333871332e-07
7.301399559e-09
0.0005093528313
7.645422963e-08
2.020202554e-06
6.813091876e-08
1.10814416e-06
1.95551597e-05
2.003879505e-06
3.670933614e-07
1.728622505e-07
0.0007263797041
1.411925493e-05
1.571311147e-07
3.392287113e-06
1.377047714e-09
6.630650805e-07
5.227904356e-08
1.127011492e-06
0.0003754771874
1.865169444e-07
6.576952888e-06
3.329497044e-06
7.07257004e-06
2.88154104e-07
-2.944075529e-07
+1.98544141222215e-07
6.929088895e-08
2.565897108e-07
6.217791051e-06
3.877741985e-07
1.707903155e-08
5.10233078e-08
4.572653519e-09
0.1095009268
7.449132821e-07
2.09740422e-06
9.378911822e-09
7.910046337e-08
1.223499404e-05
2.672193111e-10
1.02969936e-09
3.503343355e-13
1.602894349e-07
5.890632996e-07
1.487885964e-05
1.978704988e-05
2.552121299e-08
3.565602472e-08
0.0004467518296
1.005865545e-09
1.554902214e-07
2.235788016e-07
1.136525751e-07
2.962478329e-08
9.303967807e-06
3.925082515e-10
1.017880604e-05
1.434347997e-08
4.47759051e-06
1.996980966e-09
4.66863509e-05
5.845204451e-06
-2.238363039e-07
+1.509519246185828e-07
1.741949098e-06
2.215116471e-10
1.532616641e-06
1.571651059e-08
7.107477006e-05
7.234687759e-06
1.621643207e-06
2.56453036e-08
8.069391642e-09
2.035933164e-07
7.550817486e-08
8.659572881e-06
0.000432153513
3.337724504e-05
4.782082889e-05
2.323101105e-08
1.19309915e-05
0.0001350492201
2.75397447e-06
1.823126697e-08
2.411767675e-05
2.283494193e-09
3.766390067e-08
3.588343942e-05
6.787383948e-05
8.832490704e-10
5.621457713e-10
1.917361443e-08
5.849351718e-06
2.28541183e-06
6.146818441e-06
1.827359426e-05
3.231234221e-05
9.921186251e-07
3.779404373e-09
6.904752827e-06
7.698409214e-06
-1.648217403e-06
+1.111533673613149e-06
1.642479046e-06
3.301215266e-07
0.0004602920172
3.293196018e-07
3.950685019e-07
2.805273751e-05
1.739207752e-08
3.773037585e-06
1.537008391e-06
2.434913209e-08
8.315956868e-05
6.88859927e-06
-1.018611709e-07
+6.869368158598052e-08
4.672444205e-09
7.02639721e-05
7.305366536e-07
2.162269825e-05
0.0004328624246
1.653264992e-07
1.553411281e-05
5.165593047e-10
-3.193620733e-10
-1.548539649e-08
+2.153731042127394e-10
+1.044312456634803e-08
8.974979512e-09
7.520620911e-07
6.416418106e-08
1.653991504e-07
9.9043617e-09
6.542874961e-09
0.0005318553158
1.144555776e-05
8.760870822e-10
7.62010085e-08
3.426636308e-07
8.099206646e-06
1.154675564e-09
4.916072586e-06
2.013222915e-05
8.510967278e-07
3.442296018e-07
1.194853419e-07
1.251233575e-07
2.847893723e-06
2.51170281e-08
3.662564708e-07
1.646271735e-07
5.373793619e-07
4.353476735e-08
8.727968448e-08
1.053026857e-05
8.908419895e-07
1.461412043e-06
1.220147207e-06
0.0009307060856
3.48570189e-07
2.208034886e-08
8.740398335e-08
0.0001005182181
5.658935843e-06
3.660460092e-07
5.210218655e-06
1.176510494e-09
1.167510222e-09
1.74723538e-07
-6.898754196e-10
-2.43098736e-06
+4.652418777294805e-10
+1.639422266502427e-06
9.200100093e-06
0.0001764752281
7.379070775e-09
2.542944745e-07
3.71783565e-06
4.778249403e-10
6.09916388e-07
7.389397856e-06
1.030208854e-08
5.091853226e-08
1.929321992e-06
2.953703561e-07
1.428266933e-08
3.65646162e-05
8.358213962e-06
1.269052515e-05
3.142034606e-10
1.192354439e-08
3.42519028e-07
1.238073279e-05
4.592188983e-07
4.907760576e-08
2.605803561e-07
1.290266109e-06
8.91015308e-08
1.171147313e-05
6.948329894e-06
0.0001320627046
4.467361346e-08
1.343011103e-06
2.786658719e-09
1.674877668e-07
1.158457021e-07
7.247035325e-08
0.0001328992567
1.456012145e-05
1.974800428e-05
9.916997919e-05
7.076016655e-07
7.36330363e-07
8.750739052e-07
2.098702267e-05
1.200784083e-09
6.974041336e-07
6.433390899e-08
1.654868124e-05
0.0009302881479
1.821579155e-07
6.800514105e-07
1.766155755e-07
8.004650647e-07
9.941901599e-05
2.06752057e-05
1.488709742e-10
1.506772094e-06
8.10845767e-05
7.610845142e-05
0.0002115710954
2.12571763e-09
2.910738954e-06
9.48875277e-06
1.847390948e-06
7.436554496e-05
1.286439322e-10
1.413106449e-06
4.724231895e-08
5.001854028e-07
1.14055774e-06
6.718621271e-11
8.472867506e-09
9.810486272e-10
2.498427647e-07
6.987204698e-06
4.08117914e-06
2.82229796e-07
3.98589612e-09
6.278231145e-05
9.793860499e-07
9.666736058e-07
1.367103841e-08
3.490044251e-07
3.782645161e-10
8.131830282e-08
6.664325458e-06
2.365503287e-08
1.096958037e-06
3.776083109e-07
2.990539686e-08
0.002754346975
4.047644861e-07
1.411988389e-07
3.694067085e-08
0.0002320708655
3.71338597e-07
3.835998288e-05
2.623922776e-05
6.474281304e-10
0.0003219022096
-4.976400801e-07
+3.356011806202945e-07
4.178663868e-10
9.856721257e-08
2.793861685e-07
5.195375958e-06
2.632243695e-06
6.350302361e-07
-4.896870712e-10
+3.302377879428893e-10
3.043526423e-08
2.372803041e-08
6.617140209e-06
6.330476168e-07
1.870130608e-06
2.338555597e-07
4.591171378e-10
5.61933523e-05
1.955958771e-08
3.044175709e-08
5.161106748e-08
0.0001378680549
2.750711746e-09
2.289389014e-08
4.001113699e-09
2.236415233e-10
3.473425904e-06
7.311516826e-06
8.898096465e-08
1.481216022e-08
2.435208033e-12
1.189787819e-10
5.333724185e-06
4.664375968e-11
2.381386364e-10
1.249165957e-08
1.128630934e-07
1.550941521e-06
2.360856297e-14
5.846728434e-10
4.069354108e-10
4.063144762e-08
1.749554102e-08
5.155167146e-11
1.688297139e-11
2.774930205e-09
6.091667068e-08
2.87170322e-10
3.66643664e-10
8.258926512e-10
7.390910771e-09
1.48715158e-10
6.565786352e-12
3.475089984e-09
2.002380042e-09
7.274321744e-11
6.359749104e-09
7.940393308e-13
1.098611945e-06
7.438494044e-07
2.339207348e-10
4.499165346e-09
3.082882191e-09
1.878464872e-10
4.005585989e-09
2.63828057e-10
2.218789497e-09
5.827641219e-09
1.543332934e-10
4.907581247e-10
6.901620609e-10
2.088068542e-13
7.837109261e-10
1.937213353e-10
1.715044299e-08
4.069918057e-06
2.91951045e-08
1.565114165e-08
3.443423804e-08
1.256948998e-07
4.363190898e-09
2.287393528e-12
6.966944807e-09
5.995890076e-08
1.289835029e-10
1.695089305e-10
2.301282525e-11
9.445474271e-09
2.115405965e-08
1.662358545e-09
1.025068565e-09
-1.441219578e-10
-2.569402823e-11
+9.719373728155117e-11
+1.732767627245824e-11
4.41946888e-11
5.657095809e-11
5.288646625e-07
5.324238867e-07
2.221589565e-07
2.48971658e-09
1.285025395e-12
1.075626423e-09
7.782686578e-11
2.946128307e-11
3.399399058e-11
1.446405065e-08
2.187225844e-08
7.820238975e-09
2.40229777e-09
7.46594076e-10
1.281010307e-06
3.598219943e-08
2.214721626e-12
3.332304338e-11
7.617797708e-10
9.909367591e-11
4.653906139e-10
3.587191007e-06
4.446718734e-08
1.894297187e-09
4.744577059e-05
4.433959517e-09
2.210894416e-08
1.036168503e-09
1.103345617e-07
-7.407100956e-08
+4.995240386235438e-08
8.363740584e-12
3.281022886e-10
3.148002148e-08
6.330438798e-10
3.863775656e-10
5.96035067e-10
4.255627067e-11
1.165981634e-07
2.50543766e-09
1.33188621e-10
9.172189534e-07
6.582275731e-08
-2.287354562e-10
+1.542558411660617e-10
1.297572348e-09
8.107196167e-09
2.070977659e-06
6.93677041e-09
2.329665548e-07
2.684341075e-07
3.375521064e-12
4.781691571e-11
5.350983865e-09
9.177711402e-10
3.059268527e-08
1.113405482e-07
5.998651665e-10
3.911876462e-10
1.176693635e-08
8.693402036e-10
2.320788866e-05
2.036112191e-07
3.21298594e-08
7.836398482e-08
1.481220144e-09
9.13398584e-07
9.416502413e-11
1.481517933e-08
2.317808347e-06
2.429491433e-08
4.627625989e-12
2.494572686e-08
9.094749026e-10
1.806091619e-08
2.095822035e-08
7.215954213e-08
2.372708431e-09
9.374174236e-13
1.029823717e-11
5.170443463e-07
3.317414045e-07
1.109713675e-08
6.47170925e-07
1.864432618e-06
3.397272389e-11
4.168205664e-08
1.346830208e-11
3.171866066e-06
7.866352681e-11
3.289251567e-10
7.057976269e-14
2.67944287e-07
4.480095364e-10
2.542532019e-10
1.046877898e-09
4.692932633e-12
5.841384928e-10
5.254681082e-07
9.003591534e-11
3.637362188e-13
9.047013259e-09
1.769883631e-08
9.48770427e-07
2.457011926e-10
7.931261407e-08
2.511617984e-07
1.749603857e-09
1.158718175e-07
1.833652906e-06
4.450176914e-08
5.71537314e-10
1.070556403e-10
2.448737924e-10
6.262272951e-06
2.141707476e-08
4.258981326e-09
5.49246536e-09
9.566266675e-12
3.797543387e-11
2.705369876e-08
4.238910221e-07
1.784747449e-09
3.595414385e-08
3.222996111e-09
6.210164067e-05
2.258039015e-08
1.191514021e-13
1.513303142e-06
9.658853509e-11
5.082215114e-06
9.44022394e-10
8.437053744e-09
1.463048848e-10
-5.108251604e-09
+3.444930056306798e-09
7.868562747e-12
1.272388225e-11
3.427574116e-10
8.990110201e-08
2.83369365e-11
3.0757007e-10
5.265605819e-14
4.26291208e-13
7.393316606e-13
5.385884913e-05
8.367399832e-10
1.271720599e-05
8.40960929e-11
6.150536572e-08
5.764038353e-10
1.4827533e-05
2.322813811e-10
4.431637476e-09
1.319799794e-08
1.965682244e-09
5.782112632e-06
7.526132679e-07
-1.594961722e-09
+1.0756188221237e-09
8.186541654e-13
3.952010247e-07
3.136663182e-10
1.366784094e-10
2.291961377e-13
7.468862634e-11
2.13931093e-07
2.325391866e-09
7.657293509e-07
3.023206769e-07
1.974192114e-11
2.239591754e-13
1.575764712e-08
3.894663093e-07
2.463427596e-09
1.169483637e-07
7.445202399e-07
3.876016309e-06
6.204312046e-07
5.089616453e-10
2.563342578e-13
4.108963768e-09
8.534723944e-11
3.062921512e-14
5.084635637e-08
3.964438906e-09
9.56064042e-09
5.713458485e-07
8.602373058e-10
9.674381801e-09
5.575533604e-08
-1.147257577e-09
+7.736937051663231e-10
5.25696461e-10
2.604611015e-08
5.185877553e-08
5.136914596e-07
4.457967612e-09
7.087532094e-06
2.866869842e-09
2.555695269e-07
1.14983411e-06
2.121499808e-11
6.538598731e-11
4.397576097e-10
1.68339254e-10
7.632177773e-10
3.194450213e-07
4.743002807e-08
5.174172592e-09
1.378394044e-06
0.0001075134541
0.0002104593459
1.291514926e-09
1.113812317e-09
2.836821935e-07
2.119414007e-05
3.257830135e-07
6.904331208e-08
1.133565101e-07
2.515736584e-07
2.868370103e-07
1.316255824e-07
1.143196584e-07
6.920318436e-10
1.261652011e-08
7.917924522e-13
7.954494796e-10
6.551135771e-09
5.085806996e-09
1.197408695e-11
7.765435603e-11
1.137649191e-10
9.984771335e-11
1.206636661e-07
9.837079112e-09
3.488285266e-08
4.035219127e-08
1.411962856e-08
1.030986206e-07
0.0001339305599
2.858760077e-12
3.300471596e-09
1.938419056e-10
3.176459434e-12
1.739161428e-09
4.024412744e-07
3.711166771e-08
2.706304314e-09
9.299354827e-11
4.060060302e-10
1.4612779e-05
9.474942894e-11
6.71405306e-14
2.422793691e-09
1.239892236e-07
2.586764905e-08
3.834829456e-09
4.981089193e-08
4.442178086e-10
-2.539439066e-13
+1.712560508754375e-13
5.885874705e-12
3.850099718e-12
5.815334596e-10
2.871382336e-07
3.360539652e-09
4.474272729e-08
4.854260426e-09
6.835712438e-09
4.513662362e-10
1.937015923e-09
3.204387319e-07
2.834110865e-06
2.373908594e-10
1.827447632e-10
2.230101714e-07
2.935469223e-08
1.283653118e-05
1.768922904e-09
1.023937016e-07
1.810088219e-10
1.898013563e-06
1.555535859e-07
2.081843632e-11
1.130189635e-07
5.044203176e-10
5.33181246e-12
8.383860763e-08
1.081590063e-05
4.255989963e-09
4.961325451e-09
4.722044461e-12
2.605827829e-07
6.302277862e-12
7.401607027e-11
5.826773899e-09
3.775837499e-06
7.970709192e-06
1.077514227e-08
4.597223421e-10
6.830064994e-09
2.132411988e-10
3.596324757e-11
6.420547725e-09
2.577236401e-09
3.303229315e-09
6.18774819e-09
4.85351932e-08
4.830745338e-10
3.720076558e-06
1.356666145e-08
2.380000545e-09
6.934200897e-11
3.297569231e-09
3.435257108e-10
5.29774949e-08
4.382596642e-09
2.15867853e-08
5.682155814e-12
2.87290739e-07
1.68447438e-10
3.702374348e-08
4.342623412e-10
6.279260244e-08
2.403676166e-08
6.186547757e-09
1.906369995e-07
4.632234292e-10
2.755342608e-09
1.375328487e-08
2.872876599e-08
6.101397141e-07
3.268512848e-10
2.065059787e-07
2.629402879e-09
4.962476519e-11
4.062857789e-10
6.600566104e-08
2.769850168e-07
8.89254745e-07
4.790857584e-07
8.513434203e-09
2.371056219e-11
1.969882054e-08
1.194180775e-06
2.536127752e-06
3.538123306e-11
7.535922813e-11
3.279643366e-09
5.915578256e-09
3.02572337e-12
2.425014125e-08
1.736441207e-08
9.148487522e-11
2.773967859e-09
6.97620218e-08
1.281934525e-09
3.47740087e-11
7.026364385e-12
3.159133311e-08
2.31268605e-06
1.021090566e-11
6.311611218e-10
3.496581068e-09
5.006407047e-09
1.086254132e-09
6.442837702e-07
3.40398573e-08
1.040533844e-07
2.015250069e-08
2.367399296e-08
1.455786943e-08
4.519557691e-07
6.28474478e-08
1.681943257e-08
3.350671592e-10
1.788948112e-11
8.659128848e-08
1.581729673e-10
3.16899147e-08
1.269452476e-09
3.307661094e-11
6.743047887e-11
4.30353612e-10
1.049230002e-07
7.726948943e-12
5.740327635e-07
7.011998453e-07
5.337360544e-09
5.348339963e-12
6.365307921e-10
1.704476599e-14
2.368415063e-08
1.761682817e-07
8.032268925e-11
5.988468898e-06
1.432690031e-12
7.971925347e-11
1.851944696e-07
1.679623381e-10
6.480885684e-06
9.759515118e-09
3.448768806e-09
7.000972121e-11
2.170432861e-09
5.016191785e-09
1.49812676e-09
5.87915813e-10
6.048689475e-10
1.149696962e-08
8.392285319e-11
2.869985172e-13
9.120207147e-09
1.023039271e-06
1.843349283e-11
1.00493373e-06
1.437966503e-09
3.104285805e-06
3.736848159e-10
9.594801088e-11
1.232024101e-11
3.535931521e-07
7.131216184e-08
2.337326742e-10
2.82410779e-07
3.138271299e-08
5.494774408e-12
3.07020008e-12
2.011245542e-10
1.46145366e-11
3.082048934e-09
1.460968738e-10
3.186600527e-07
8.613894657e-10
4.227362638e-07
1.132945088e-08
1.385877833e-08
2.535124307e-09
5.647614453e-08
1.712204167e-10
9.689130473e-07
2.23306088e-09
2.607948681e-11
1.678634611e-13
2.258659145e-06
1.695727868e-09
1.706090682e-11
4.758749529e-12
1.426536131e-05
1.449437791e-07
5.476310283e-08
3.617494916e-10
8.345680248e-12
7.09135827e-13
8.179812515e-09
3.323888637e-09
1.436305431e-08
5.885631636e-10
3.584884074e-08
9.15498181e-09
9.640861937e-11
1.908907816e-09
2.558536978e-10
5.232929043e-10
1.261308954e-09
2.535469757e-10
1.402765641e-06
9.405064753e-10
1.920855926e-07
4.40419547e-11
5.00514635e-10
1.406592654e-07
4.52773411e-10
4.146849421e-07
1.166300995e-07
7.051208905e-09
2.309332429e-09
9.269115967e-07
-1.602436999e-10
+1.080660039831632e-10
1.045618377e-06
4.94758944e-07
3.094362842e-13
1.351633908e-08
2.784256537e-09
1.934517844e-08
-7.880756508e-13
+5.314666752529101e-13
1.357165448e-09
4.977436192e-12
6.477152561e-11
4.234578629e-08
4.918295292e-13
3.683451803e-08
3.727223271e-09
1.040244287e-09
6.604553745e-10
3.806952288e-09
1.097089218e-08
1.705682648e-08
6.319891143e-07
1.574130537e-09
3.013782963e-08
4.772880086e-09
-2.994218614e-10
+2.019257173887625e-10
8.938611559e-10
5.601588127e-09
1.526952592e-05
3.497239887e-10
9.33354413e-08
1.359301132e-08
6.404370059e-11
1.535628496e-11
-3.32672467e-09
+2.2434943877662e-09
4.662210325e-09
1.022052539e-10
1.383917366e-07
1.619072147e-12
1.507995029e-08
2.056138551e-09
1.984766447e-05
5.313284134e-09
1.04749306e-08
1.156452533e-07
1.201443255e-06
1.584111891e-09
2.09858386e-11
3.445445203e-09
4.164676e-06
4.092776954e-12
9.944337818e-10
5.276009823e-09
3.628978704e-11
2.197020798e-10
1.052815682e-08
1.068627029e-08
3.831431452e-11
1.364029236e-07
3.403124343e-10
5.38162347e-08
6.729323819e-12
1.207769187e-07
4.017542328e-09
1.315683749e-08
7.126836658e-11
1.491126652e-07
7.209030864e-10
6.628062787e-08
1.21591155e-10
6.230923352e-09
0.000101234359
5.30942773e-10
5.741385346e-10
2.429893426e-08
7.325998025e-12
3.533987787e-09
7.926231814e-10
8.291189826e-11
3.746302962e-11
5.242302307e-11
8.718447375e-08
6.176050515e-11
6.09961364e-09
2.644343829e-09
2.306901235e-08
7.322344218e-10
1.884045094e-06
3.755632286e-08
1.659778169e-08
5.495025114e-08
3.546727594e-09
2.674246643e-07
1.028252112e-09
9.289249703e-06
7.959960572e-08
1.068081079e-11
5.674343265e-07
6.214463157e-09
1.601066932e-06
1.977738539e-06
3.462786645e-07
7.373697384e-09
5.726726673e-12
diff --git a/t/ME_data/ME_Wp.dat b/t/ME_data/ME_Wp.dat
index d9e12fc..53c7f25 100644
--- a/t/ME_data/ME_Wp.dat
+++ b/t/ME_data/ME_Wp.dat
@@ -1,1289 +1,1289 @@
0.0001355688725
3.608722302e-05
4.306341977e-07
9.75550325e-07
7.965530581e-08
1.158741173e-06
2.661843554e-06
2.305034582e-09
2.352135748e-07
1.82973514e-06
7.663338515e-10
1.356607225e-07
3.025163564e-07
9.674844388e-08
1.08577046e-08
1.140610128e-07
9.611408862e-06
9.926076721e-06
0.0002395617463
3.935060204e-07
1.135288401e-06
3.820198516e-06
9.707598452e-05
1.989871939e-08
4.670964333e-06
6.369460204e-08
5.769260099e-06
6.100771974e-07
4.177606056e-07
5.131150168e-05
2.623663346e-08
9.921267961e-07
2.81862487e-05
2.357728363e-06
1.52940914e-07
3.812119407e-09
4.089890436e-06
5.519579034e-08
3.108467376e-05
6.275083181e-07
0.0004711631212
3.321145292e-11
3.599907229e-07
0.0005693438579
5.344259311e-06
1.369519561e-07
4.610580273e-08
2.759136384e-08
0.0002011663129
8.79267475e-06
2.70389756e-06
1.871332754e-09
6.64520168e-07
0.0134591178
8.026976255e-06
1.13118716e-07
2.178470183e-07
3.510890511e-05
1.899844473e-07
1.045126454e-06
1.869095429e-06
2.768546174e-06
8.551014192e-08
7.773400333e-07
4.53144298e-06
1.601156732e-06
1.364165622e-07
1.896432416e-08
1.972891892e-07
7.870126436e-08
1.027083098e-09
1.68746816e-07
6.312089294e-10
3.505278148e-07
4.91589775e-06
1.471309784e-08
5.806965344e-07
1.107356095e-06
2.913951197e-05
2.124037647e-06
2.12990617e-07
1.356428765e-06
2.199225656e-08
9.711357262e-08
0.0001385574064
6.515223743e-07
2.535298939e-05
9.498820108e-06
7.275424133e-08
3.660521443e-08
1.057452388e-07
1.476424382e-10
1.535388646e-05
8.248395671e-07
3.681439971e-07
6.978801279e-07
2.186342703e-08
0.00317235351
0.0002184140294
1.851577146e-06
2.387181833e-07
0.0002408737015
2.563590676e-07
1.00063893e-09
4.85334472e-06
1.03755407e-07
7.119588285e-06
2.561890817e-06
6.552784786e-08
2.501651241e-08
3.220116707e-08
4.669231392e-07
1.713317361e-05
1.245254933e-09
7.470170779e-06
1.12179479e-10
2.566806243e-08
0.000690252632
1.619690789e-08
0.007842918829
3.718578613e-05
0.0007152644659
3.366872929e-05
3.590280598e-05
0.0001083353837
1.746053028e-07
2.714474998e-05
5.49513934e-05
6.535082158e-05
5.041690215e-07
3.546092261e-06
1.20179547e-08
2.615186068e-08
0.005544506156
2.669153415e-08
1.124413241e-07
2.097675904e-07
0.0004497021745
3.522902236e-06
3.609795144e-10
2.381543307e-06
0.0001723836125
8.005109609e-06
2.55897517e-05
1.661304559e-05
7.830076033e-06
8.689122826e-05
1.60422181e-05
1.56327932e-07
7.475227738e-08
3.375304186e-08
7.208470216e-08
-8.191390486e-10
+5.524153757135919e-10
5.4103721e-09
9.793444055e-05
1.21018852e-06
5.403237195e-06
9.396037827e-07
-1.366976742e-06
+9.218690916306924e-07
1.916819566e-05
2.177871221e-06
0.0003020897612
6.197213888e-08
2.556142456e-08
1.364194249e-05
1.038835664e-05
1.181892855e-07
5.005165507e-06
-5.35483996e-08
+3.611225631803286e-08
1.221974589e-08
1.942760585e-05
3.734309313e-06
0.0001093995311
2.068051727e-06
8.857501175e-06
4.173118322e-06
2.109563994e-08
1.44827703e-06
4.552053973e-07
1.371432834e-06
8.524886363e-07
1.659141663e-09
1.245863146e-07
0.0002242637793
3.699075972e-06
2.084644637e-10
2.212848552e-08
8.425415515e-08
2.259085562e-07
0.007259913649
4.467603509e-06
3.025351151e-06
1.224940717e-08
1.15612173e-07
1.910495778e-06
6.957098183e-08
3.814134721e-08
4.664811457e-07
1.051815649e-05
3.412719008e-13
6.131289047e-07
6.896616304e-06
2.90239901e-05
5.756729461e-09
1.528995487e-08
5.225540705e-09
3.429458529e-05
2.622252694e-09
0.0005682063516
7.453799242e-07
1.807444313e-05
9.740908615e-07
2.028479198e-06
9.039653406e-08
2.480442463e-06
0.0004804561921
1.082404311e-08
2.792579136e-06
0.0002121921202
9.812552393e-05
7.828388801e-07
0.001456303522
1.943547444e-08
1.718825287e-07
-1.018305194e-09
+6.867301068547712e-10
1.796235147e-07
3.651695186e-07
1.227417157e-07
-9.855784321e-06
+6.64659657953662e-06
8.406970187e-05
0.001885290476
0.0001570330885
2.314708659e-06
0.0002982727873
0.003229609179
2.42756184e-10
1.896039184e-06
-1.765970661e-08
+1.190944745856344e-08
4.219946428e-07
7.272059576e-08
1.591962088e-06
2.437096009e-07
2.465458718e-10
3.867407921e-06
7.355333709e-09
3.882793753e-06
0.003249948478
1.636247528e-05
0.0003223792165
2.270764742e-05
1.158210386e-07
4.803440794e-06
-2.87447101e-06
+1.938501144056279e-06
1.00011733e-06
7.966528549e-08
1.977164344e-05
3.697861147e-06
0.0001118273732
8.542125574e-09
3.549017984e-05
4.419576384e-08
0.000415723574
1.269135249e-06
1.336086492e-09
0.0001810768552
5.607526182e-08
7.073836906e-08
6.226377541e-05
3.842918719e-07
6.352284838e-05
1.767534549e-08
0.002554251878
1.646540099e-08
1.960645506e-08
0.0002646807702
1.199915361e-05
4.194897886e-08
0.000116967724
1.112499069e-08
4.229713001e-07
4.907009816e-08
1.380172172e-06
1.086563506e-05
4.639365797e-08
5.978111043e-07
7.468593061e-08
2.128107648e-08
3.566631743e-10
1.4537143e-09
3.088717592e-05
3.041172443e-05
2.604197881e-07
1.592445808e-06
1.544432896e-07
3.372379695e-06
0.0001468102323
5.574768414e-07
4.759127919e-10
0.00093067198
5.974247937e-06
0.0002033151473
2.248619959e-06
1.737481127e-05
3.508166199e-06
2.435498983e-09
6.369004222e-06
9.703729079e-09
1.740708919e-06
3.283306645e-07
1.701425138e-08
4.192916377e-08
4.498845691e-05
3.834128247e-05
1.152621968e-09
1.27759416e-05
0.000281123799
1.725068848e-07
1.077501014e-06
0.004189517285
1.659706255e-07
1.343358849e-07
8.327570633e-07
1.678461499e-08
5.667318085e-07
1.495443328e-09
6.283129102e-07
7.691168852e-09
1.613371844e-08
4.367733457e-08
1.328415042e-09
1.452913836e-06
2.025565342e-07
1.29314441e-06
3.674495466e-07
4.037738642e-06
2.086123437e-08
1.579986535e-06
7.231915707e-10
3.157015643e-07
0.0003660037209
0.0008949054777
5.638174345e-08
3.288585027e-06
5.432504634e-08
5.180759332e-10
5.121204016e-08
1.939451095e-08
1.979142335e-05
2.529002486e-08
0.0002538283082
3.80397913e-05
7.172326439e-06
2.759121311e-06
3.708767178e-06
4.71298987e-08
0.0001084418497
5.277822142e-08
-6.471171803e-08
+4.364063474995147e-08
1.209897825e-07
3.136890863e-08
1.044095011e-07
6.962956326e-08
6.454347396e-06
2.03105437e-09
4.064764112e-07
1.125205343e-06
3.916697214e-06
5.611785349e-07
1.848611741e-07
2.582403728e-07
8.839788397e-05
2.850326226e-06
2.108292593e-07
1.232805123e-08
9.539264139e-08
0.006539873531
2.250970922e-05
3.287682183e-08
3.209268362e-11
3.531432382e-10
1.155035491e-07
1.892464014e-07
1.339909548e-08
5.52939388e-08
1.151900346e-05
0.005889722048
8.216256879e-05
1.557120111e-07
9.089312296e-08
4.410189849e-07
4.12565868e-08
1.888812372e-07
5.726907174e-05
1.874558888e-09
9.008463668e-12
1.709044906e-05
1.67931414e-07
8.335566347e-06
1.14159826e-08
0.002680456252
1.068398085e-08
1.957173741e-06
2.244367604e-05
1.231910962e-05
4.629935916e-05
4.576226372e-06
1.980987784e-08
3.93911073e-07
2.650096111e-05
8.453749094e-06
3.871976822e-08
0.002862691182
4.135291611e-08
0.0002240310514
1.626024307e-07
1.446697749e-07
1.81572275e-06
3.82343645e-07
1.519921801e-06
7.180741242e-10
3.48543689e-05
7.358200622e-08
7.791803572e-06
7.340099307e-07
1.64882977e-05
5.701378479e-07
1.101811989e-07
2.180113388e-08
0.005048614688
1.4109536e-08
1.796629531e-05
3.537585755e-05
2.035011121e-05
1.003894855e-08
0.0007063054044
8.670678472e-05
6.82639369e-08
5.900706666e-06
1.447352541e-08
-2.170801903e-06
+1.463957005451528e-06
2.460311e-05
3.389946454e-08
3.777468751e-07
7.004890999e-07
3.633307241e-05
4.394910201e-07
7.401664237e-07
0.004667703381
4.822543865e-07
1.289604844e-08
0.0001296635825
1.387032601e-06
2.915976119e-05
1.716214639e-07
2.901236713e-08
7.505249455e-10
1.054930074e-09
5.765824783e-07
9.011588796e-11
1.700596863e-06
3.410463387e-08
7.817851611e-06
0.000171359498
4.092934015e-08
9.55513037e-08
2.412195475e-06
3.183516928e-09
0.0006930495435
1.228003613e-07
8.804634278e-05
2.13517622e-06
4.693017624e-06
5.984372321e-08
2.487744823e-06
1.776206816e-07
6.533684691e-07
1.015310185e-06
4.483048832e-06
1.390213091e-06
2.463272722e-09
0.0001373179166
4.692519624e-10
6.223649822e-09
3.026899135e-06
5.540599649e-08
5.422865018e-08
8.715360663e-08
3.005262078e-06
4.491518597e-06
1.296331911e-08
8.034189259e-07
8.445679444e-07
1.331376027e-08
1.808292903e-06
4.811816349e-06
3.356636075e-06
3.718204067e-07
1.63321067e-06
1.145748524e-06
1.082083666e-06
6.645503138e-07
7.690061877e-06
9.446655788e-07
0.0001628166399
0.0004161528763
2.733182188e-06
4.797973294e-06
0.0001243192011
0.0004489774904
4.018994469e-08
6.294131882e-05
0.0006300938801
3.276964999e-07
5.253447013e-07
-8.917191585e-09
+6.013623386455616e-09
6.133379962e-07
7.498082273e-08
9.70785668e-08
2.892020441e-07
0.0008858500922
5.834224897e-12
2.657258192e-05
1.937050603e-08
1.014207639e-08
2.078416478e-09
9.09266244e-09
1.171070677e-06
6.186008954e-07
7.316626453e-08
1.137928043e-05
3.926244745e-07
9.411347464e-08
2.358883261e-05
3.718597508e-06
8.746279414e-07
2.158976175e-07
7.682516501e-08
4.623476813e-08
1.822643012e-09
5.649636691e-07
8.861102166e-07
3.661799027e-06
0.005151284527
1.333261465e-08
0.01070078066
3.557312307e-07
7.956062001e-07
1.102026215e-05
9.099085434e-08
9.760469392e-07
4.026083515e-07
3.148718342e-05
3.344966932e-10
5.83259011e-08
0.000275860311
5.979810844e-06
2.336312158e-06
1.214278297e-06
1.257019376e-05
4.380094449e-07
8.01469852e-07
0.003169580958
9.221636814e-05
4.690760304e-07
7.648050266e-07
1.418690988e-05
2.220802016e-07
3.14201845e-07
1.911591064e-07
2.452984351e-08
0.0001445594492
3.388666373e-11
1.315444637e-09
3.270418535e-05
2.302140766e-07
4.002096106e-06
-1.469589992e-06
+9.910699640329833e-07
1.156198845e-08
4.49026944e-08
0.0002105189957
1.502882545e-10
0.0003938508954
2.359522705e-06
-8.550778195e-05
+5.766519563931604e-05
4.012744154e-07
1.060136498e-07
0.0002141768308
2.679472087e-06
1.069266147e-05
2.48432308e-09
0.0001019831675
1.193130845e-06
1.086711459e-06
1.985358343e-07
-1.739018449e-10
-7.388141233e-07
+1.172768568716882e-10
+4.982454226859217e-07
2.574857412e-05
0.001250253279
2.058059351e-09
1.333980604e-07
9.818799108e-09
5.15408822e-07
0.0001083932826
1.428558882e-06
1.74625784e-05
2.151559189e-05
0.01028504249
7.324000358e-07
2.685867501e-08
4.250514247e-05
0.0004356795916
3.289457433e-07
0.00323686003
3.007561717e-08
2.793887221e-06
1.897911806e-05
1.61972882e-07
8.953729168e-05
4.802035152e-06
7.858920744e-06
7.610825707e-07
2.314286712e-05
4.93977485e-09
4.92614674e-09
8.469370311e-07
3.80925189e-05
8.402086849e-09
4.426792197e-07
2.426545735e-06
3.814110561e-08
1.746240227e-06
1.52446212e-07
7.921180881e-05
2.10135866e-07
4.839236788e-07
5.352097408e-07
4.622274946e-05
0.001900140645
-2.85310478e-06
+1.924092071480528e-06
0.001430541179
4.980473899e-08
0.001646505316
1.600615719e-05
2.692315114e-06
4.189143283e-06
9.988552464e-08
3.325410404e-05
1.063381388e-06
4.230272945e-06
3.9885381e-06
0.0001050685482
5.986074981e-08
2.025564094e-06
4.345191943e-07
8.910595474e-05
2.543279105e-06
0.0002129944726
8.738486781e-06
4.913962453e-09
4.914524134e-10
3.330632607e-05
1.95151118e-06
1.384784661e-05
-3.18562271e-07
+2.148337292837028e-07
1.685318586e-07
1.114845673e-07
6.020836887e-05
1.280889433e-07
4.215297019e-07
3.062431208e-07
6.986710594e-07
1.615334662e-07
1.317663183e-07
0.0008030217322
2.015927418e-05
5.19186806e-06
5.500044264e-07
3.714750104e-08
0.0002157285841
9.440044669e-07
6.51295213e-08
1.03903577e-08
0.0001795918463
1.962893454e-05
7.871495075e-10
2.12576483e-08
1.743490745e-05
1.541416437e-06
3.391980216e-06
8.359096999e-06
1.415613908e-09
7.887561636e-07
1.745081946e-07
1.932727685e-07
9.828062711e-07
1.291461662e-07
9.278171187e-05
4.264104098e-06
3.754457748e-07
3.452634763e-05
8.101309385e-06
4.056256865e-09
2.25304955e-07
1.172523704e-05
4.796449304e-07
-5.612956712e-09
+3.785295788394707e-09
2.390481542e-05
5.308836572e-08
4.670595432e-05
6.758171771e-05
-1.509575766e-09
+1.01803578427425e-09
1.661568152e-05
5.893560303e-08
3.900939999e-05
3.504262189e-08
8.675964674e-05
0.0001335992511
1.302314855e-06
1.146394419e-06
6.394209116e-06
2.519496842e-05
0.0001367050958
4.154554201e-09
2.099823438e-06
2.146448613e-05
1.233863159e-09
0.001728586801
4.09215386e-05
0.007151464488
4.760189124e-09
0.0007525692076
1.763213602e-08
2.706320225e-07
4.168939514e-06
6.813774053e-07
6.199680386e-06
1.450358601e-09
0.001178944008
3.552524877e-07
0.001006338796
3.68488544e-07
2.167674892e-07
2.537842791e-06
7.658583099e-07
-3.842627831e-07
+2.591411922436858e-07
4.604369581e-08
4.168373998e-07
7.238778845e-09
1.170127806e-06
4.393582652e-07
0.0007919330066
1.02877764e-06
6.918341723e-08
9.436339417e-08
1.591183988e-06
2.089478745e-06
1.613369096e-07
3.637102293e-07
9.604977934e-07
9.357975878e-07
2.88751538e-06
6.126247948e-06
1.015660592e-05
5.755046456e-06
0.002323684387
8.012994582e-07
5.26519905e-07
0.0002032569161
1.451031711e-07
5.477022502e-07
7.468911459e-07
1.172272701e-05
5.000191366e-06
1.572889827e-06
6.442064306e-08
2.908881721e-05
4.569175874e-05
0.0003671971117
1.009483851e-07
1.696662855e-09
1.138346461e-08
7.716922378e-07
3.529252959e-06
6.363107685e-07
1.932801483e-07
5.985897033e-07
3.856271125e-06
8.483829284e-07
7.519170568e-09
6.046535116e-06
8.332016051e-05
0.0004929003218
1.066348909e-06
6.354668796e-05
2.325179903e-05
3.330122089e-07
7.907070499e-09
3.514991564e-06
6.69519932e-06
7.791955281e-08
3.007406815e-09
1.687262337e-06
2.966015995e-06
2.551884972e-08
0.02757085539
1.902124521e-06
1.858863499e-06
1.174589642e-07
6.265121517e-07
0.0001498056382
2.510077007e-07
0.003704558135
7.86457155e-05
1.776365133e-08
0.0002011705313
4.239853347e-06
1.056573955e-06
0.0008521527377
4.768494505e-07
5.244547978e-06
0.0001425058496
4.507288271e-07
3.092538831e-07
0.007681533942
1.571813043e-07
3.341399818e-08
9.577252789e-07
1.135356146e-09
1.20086557e-07
4.167792379e-06
7.827764185e-06
1.444125375e-06
2.567040739e-05
3.703109657e-07
3.826055602e-08
2.223174582e-07
1.771598439e-05
1.154475457e-07
4.676107532e-07
1.045592241e-08
5.174223151e-07
8.734482888e-07
1.146406702e-08
8.784052908e-07
7.500525642e-10
1.921656528e-05
5.899142263e-08
7.30476082e-05
1.070817446e-08
2.432332991e-05
-3.657113723e-09
+2.466303925798377e-09
1.587801136e-05
5.253580132e-08
1.186388137e-06
2.513054181e-08
2.878235795e-05
1.268594957e-09
-5.985219895e-07
+4.036344625890604e-07
1.480245097e-05
0.007764649525
1.266626553e-05
2.673316923e-07
2.052285557e-07
3.026411522e-09
3.958784188e-09
5.153150957e-10
2.304903761e-10
4.00173635e-09
5.422672396e-12
1.0797419e-08
7.876103608e-07
2.04929605e-13
2.452584831e-10
6.723994081e-09
1.881390676e-09
8.769218159e-09
1.325727129e-07
2.663023215e-10
8.598804962e-09
1.838865267e-09
1.411793749e-09
3.180735014e-09
1.657163946e-09
6.291662125e-07
1.407385591e-08
2.446407045e-10
3.694683988e-06
7.168303932e-08
1.515131945e-07
5.976833117e-07
9.943757029e-11
1.770333656e-10
6.826832537e-10
6.831079527e-10
1.003726202e-11
2.246655664e-10
1.56313715e-08
1.258803718e-08
1.510698504e-09
1.860286446e-08
7.436041453e-11
8.891989836e-11
8.343726115e-07
7.82113065e-10
7.372199042e-09
1.73794743e-08
4.796628283e-10
1.006843364e-11
2.297522292e-07
2.964339936e-09
5.119472399e-06
1.217098723e-09
1.247855962e-10
2.379923838e-10
6.847023572e-08
1.126084853e-08
9.36637247e-09
8.815353934e-13
1.248942677e-09
1.288006978e-06
1.862006179e-05
5.986024169e-09
3.576379278e-09
1.396207867e-07
5.597776909e-09
1.227869194e-07
3.184774841e-10
1.281546379e-09
3.785370242e-08
2.653293919e-11
1.797326593e-08
4.870914289e-09
1.481397924e-08
6.148653826e-13
2.60768604e-08
7.503399778e-09
7.319978311e-10
3.917579972e-06
2.456000696e-09
1.6358728e-09
1.504286785e-07
4.626076109e-12
2.544705854e-10
2.326334049e-09
1.209466619e-07
9.766234969e-09
2.844168649e-07
1.201220077e-06
7.29513531e-08
2.673472495e-08
3.217664769e-11
3.051103391e-08
8.284753217e-08
2.058274741e-10
2.829620954e-09
6.31774405e-05
5.647654886e-08
3.939849241e-09
6.945357028e-09
2.688851327e-09
1.246699666e-07
1.659903866e-07
2.236198294e-05
9.480526263e-09
9.500853178e-07
7.271132888e-08
4.043174967e-07
7.378138213e-10
1.151223659e-07
3.256246906e-08
1.170385775e-12
5.503748995e-08
1.04382963e-11
1.203555924e-08
1.075568292e-10
8.108085296e-07
4.847854101e-11
7.343632975e-06
1.187566413e-11
6.085942402e-09
7.018111299e-12
2.463652527e-10
1.851835159e-10
1.513115041e-09
8.018679868e-09
3.929445577e-06
2.084075835e-07
5.951216671e-10
3.507041185e-06
5.267975954e-12
3.372609083e-06
1.694025855e-11
1.867026183e-10
2.116504705e-10
3.831483929e-10
6.969465117e-05
1.539468292e-08
1.101455262e-11
1.388260351e-10
5.247893014e-08
4.682122453e-08
6.387349339e-14
2.19360099e-08
1.553226092e-09
2.996160973e-08
2.723720456e-11
1.94231476e-09
1.292957006e-11
3.828519911e-12
1.610161407e-08
1.687340225e-06
1.925507425e-09
4.366763058e-13
6.454545293e-07
6.024146742e-07
2.580279916e-10
1.225824367e-05
1.999348472e-11
4.860969609e-11
2.362388708e-08
1.373681959e-08
3.790417277e-12
2.360779892e-11
3.1293099e-07
1.364398624e-09
1.715330241e-05
3.172978773e-08
1.823609884e-08
3.114529471e-09
3.261473002e-10
2.126205049e-05
5.67083409e-08
6.583915678e-14
7.465321603e-10
1.354621048e-10
2.460861072e-08
1.087512798e-09
2.736687033e-07
5.113035943e-10
6.306214627e-07
4.283779224e-11
2.044825034e-08
6.895612973e-10
6.338164558e-12
5.398174211e-09
1.015670208e-07
5.845735119e-10
4.836785182e-08
5.100824595e-10
1.572989395e-08
1.787382013e-10
1.632652152e-08
7.01728011e-08
7.320289671e-12
3.029928398e-11
1.409128876e-12
1.238566734e-06
8.502740668e-08
3.438552357e-08
4.863721359e-07
5.154442037e-07
2.08972413e-09
3.431301186e-10
5.565393637e-13
1.544285017e-09
4.30867667e-09
2.652046726e-09
1.604339305e-11
1.019689627e-05
1.123301038e-09
2.815571339e-12
1.072295152e-10
6.773569516e-09
2.399435949e-09
1.754778047e-07
6.118111465e-11
8.592783136e-06
6.887731648e-11
1.799971692e-09
2.674451282e-10
1.267736949e-08
1.616784109e-09
2.41843279e-10
2.161775938e-07
3.594767539e-06
1.314807907e-05
3.035193063e-06
5.146951027e-12
8.12518692e-07
2.904383885e-07
8.278942839e-08
1.491259815e-12
3.232983294e-10
3.014861189e-07
3.869348398e-08
4.060145232e-09
2.722588887e-09
1.252225008e-06
8.534088341e-10
1.170269802e-13
1.013048217e-07
1.173680983e-10
1.230880253e-06
3.603057796e-07
2.340496864e-10
1.622753822e-07
1.855111476e-10
1.873444017e-13
1.477986621e-10
5.974804124e-09
5.502900798e-08
6.238803852e-09
2.281648706e-08
2.708668262e-09
8.942323868e-07
8.600240612e-08
5.313999062e-08
1.150328412e-06
4.275682161e-10
2.535797082e-12
1.472728159e-07
2.05454394e-10
2.452447725e-06
2.149816554e-09
1.302067004e-08
1.354519139e-05
1.493773984e-10
4.155192241e-09
1.975441615e-09
4.417745266e-09
3.017038977e-08
6.10674516e-07
7.43296019e-09
3.802076048e-08
1.356542979e-08
7.945938899e-11
3.638927521e-12
1.061272859e-06
3.425627163e-09
4.724471432e-09
2.264988832e-06
5.206993128e-12
1.903839604e-11
5.585175902e-09
9.146019191e-13
8.796344869e-08
2.018357672e-09
7.316136952e-07
6.365052389e-11
5.084402628e-08
4.955636143e-10
7.579083396e-10
1.041589492e-10
1.067922679e-10
2.629804234e-10
7.253503328e-09
3.644406678e-09
3.524116213e-06
1.911509485e-10
7.70936902e-09
9.497529697e-06
1.776186622e-07
7.672088841e-12
1.053749536e-07
4.528511033e-07
2.689071764e-11
2.172737406e-09
1.388580447e-10
8.063951339e-10
5.664084936e-08
1.744805297e-09
2.052101897e-11
2.110860659e-11
4.471965163e-07
1.547872897e-09
1.123990826e-08
1.220466442e-08
7.869401405e-07
1.643506226e-10
1.008594225e-08
3.353990229e-07
6.498833539e-10
1.784315183e-10
3.805120267e-07
8.427310035e-07
1.498076756e-09
1.039474025e-10
1.111563184e-08
1.747513597e-09
2.355165317e-09
1.678233595e-08
3.217683158e-12
5.649816397e-08
7.451025295e-08
3.278979106e-11
3.654252804e-11
1.935635959e-11
4.454493836e-09
1.468843362e-11
6.571277705e-09
6.056988847e-10
1.451019442e-07
1.129031793e-10
1.042132383e-11
3.257388686e-11
9.125965424e-09
1.157787061e-08
1.347062065e-10
1.77255365e-07
2.066007246e-07
1.183381206e-10
4.34066107e-11
1.724089241e-07
1.017527409e-09
4.44819419e-08
9.298034227e-10
1.457982094e-08
1.344707431e-06
3.942998333e-10
6.431805784e-10
2.598344406e-09
9.723411144e-11
5.888094133e-07
1.250960845e-08
1.918965372e-08
4.048461902e-10
2.358561388e-08
9.312555038e-07
4.777092354e-09
1.101320525e-05
4.464641617e-10
7.607867556e-08
5.154678729e-10
2.264518558e-08
8.764203667e-09
8.615359173e-07
1.400305016e-06
2.971220423e-07
1.703068601e-09
9.085894816e-08
7.676840903e-09
3.495641886e-07
5.3486215e-10
8.550499106e-10
5.359621176e-10
3.270585775e-08
2.199341422e-11
3.481833841e-07
1.663118365e-12
1.079199186e-09
2.042388686e-11
1.522460286e-08
1.814687048e-09
3.79827064e-08
3.251582144e-06
1.339564566e-08
1.067207794e-10
2.907408741e-09
2.121056487e-09
-2.860605344e-10
+1.929150342315743e-10
5.485146663e-10
2.944552971e-10
3.400695511e-09
2.722315508e-07
1.103475015e-07
4.918180357e-10
1.403530754e-08
1.186111191e-10
2.383851903e-09
6.366433869e-10
2.884801368e-08
3.233653883e-09
0.0001662977498
2.949009865e-10
2.818141905e-08
4.590633891e-10
2.878567552e-11
1.429984095e-11
2.151490092e-11
4.991815907e-07
7.608561997e-11
4.322238014e-08
8.515965863e-10
4.522028585e-09
4.658966814e-08
4.872474226e-10
1.225955173e-08
9.732456294e-09
-8.233972895e-08
+5.552870709834466e-08