diff --git a/include/HEJ/MatrixElement.hh b/include/HEJ/MatrixElement.hh
index 29a6306..7f9e372 100644
--- a/include/HEJ/MatrixElement.hh
+++ b/include/HEJ/MatrixElement.hh
@@ -1,187 +1,187 @@
/** \file
* \brief Contains the MatrixElement Class
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <functional>
#include <vector>
#include "fastjet/PseudoJet.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Parameters.hh"
#include "HEJ/config.hh"
namespace CLHEP {
class HepLorentzVector;
}
namespace HEJ{
class Event;
class Particle;
//! Class to calculate the squares of matrix elements
class MatrixElement{
public:
/** \brief MatrixElement Constructor
* @param alpha_s Function taking the renormalisation scale
* and returning the strong coupling constant
* @param conf General matrix element settings
*/
MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
);
/**
* \brief squares of regulated HEJ matrix elements
* @param event The event for which to calculate matrix elements
* @returns The squares of HEJ matrix elements including virtual corrections
*
* This function returns one value for the central parameter choice
* and one additional value for each entry in \ref Event.variations().
* See eq. (22) in \cite Andersen:2011hs for the definition of the squared
* matrix element.
*
* \internal Relation to standard HEJ Met2: MatrixElement = Met2*shat^2/(pdfta*pdftb)
*/
Weights operator()(Event const & event) const;
//! Squares of HEJ tree-level matrix elements
/**
* @param event The event for which to calculate matrix elements
* @returns The squares of HEJ matrix elements without virtual corrections
*
* cf. eq. (22) in \cite Andersen:2011hs
*/
Weights tree(Event const & event) const;
/**
* \brief Virtual corrections to matrix element squares
* @param event The event for which to calculate matrix elements
* @returns The virtual corrections to the squares of the matrix elements
*
* The all order virtual corrections to LL in the MRK limit is
* given by replacing 1/t in the scattering amplitude according to the
* lipatov ansatz.
*
* cf. second-to-last line of eq. (22) in \cite Andersen:2011hs
* note that indices are off by one, i.e. out[0].p corresponds to p_1
*/
Weights virtual_corrections(Event const & event) const;
/**
* \brief Scale-dependent part of tree-level matrix element squares
* @param event The event for which to calculate matrix elements
* @returns The scale-dependent part of the squares of the
* tree-level matrix elements
*
* The tree-level matrix elements factorises into a renormalisation-scale
* dependent part, given by the strong coupling to some power, and a
* scale-independent remainder. This function only returns the former parts
* for the central scale choice and all \ref Event.variations().
*
* @see tree, tree_kin
*/
- Weights tree_param(
- Event const & event
- ) const;
+ Weights tree_param(Event const & event) const;
/**
* \brief Kinematic part of tree-level matrix element squares
* @param event The event for which to calculate matrix elements
* @returns The kinematic part of the squares of the
* tree-level matrix elements
*
* The tree-level matrix elements factorises into a renormalisation-scale
* dependent part, given by the strong coupling to some power, and a
* scale-independent remainder. This function only returns the latter part.
* Since it does not depend on the parameter variations, only a single value
* is returned.
*
* @see tree, tree_param
*/
double tree_kin(Event const & event) const;
private:
double tree_param(
Event const & event,
double mur
) const;
double virtual_corrections_W(
Event const & event,
double mur,
Particle const & WBoson
) const;
double virtual_corrections(
Event const & event,
double mur
) const;
//! \internal cf. last line of eq. (22) in \cite Andersen:2011hs
double omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const;
double tree_kin_jets(
Event const & ev
) const;
double tree_kin_W(
Event const & ev
) const;
double tree_kin_Higgs(
Event const & ev
) const;
double tree_kin_Higgs_first(
Event const & ev
) const;
double tree_kin_Higgs_last(
Event const & ev
) const;
/**
* \internal
* \brief Higgs inbetween extremal partons.
*
* Note that in the case of unordered emission, the Higgs is *always*
* treated as if in between the extremal (FKL) partons, even if its
* rapidity is outside the extremal parton rapidities
*/
double tree_kin_Higgs_between(
Event const & ev
) const;
double tree_param_partons(
double alpha_s, double mur,
std::vector<Particle> const & partons
) const;
std::vector<int> in_extremal_jet_indices(
std::vector<fastjet::PseudoJet> const & partons
) const;
std::vector<Particle> tag_extremal_jet_partons(
Event const & ev
) const;
double MH2_forwardH(
- CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
+ CLHEP::HepLorentzVector const & p1out,
+ CLHEP::HepLorentzVector const & p1in,
pid::ParticleID type2,
- CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
- CLHEP::HepLorentzVector pH,
+ CLHEP::HepLorentzVector const & p2out,
+ CLHEP::HepLorentzVector const & p2in,
+ CLHEP::HepLorentzVector const & pH,
double t1, double t2
) const;
std::function<double (double)> alpha_s_;
MatrixElementConfig param_;
};
}
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 1396d29..830ab3e 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1755 +1,1726 @@
/**
* \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 {
+ Weights MatrixElement::operator()(Event const & event) const {
return tree(event)*virtual_corrections(event);
}
- Weights MatrixElement::tree(
- Event const & event
- ) const {
+ Weights MatrixElement::tree(Event const & event) const {
return tree_param(event)*tree_kin(event);
}
- Weights MatrixElement::tree_param(
- Event const & event
- ) const {
+ Weights MatrixElement::tree_param(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
- Weights MatrixElement::virtual_corrections(
- Event const & event
- ) const {
+ Weights MatrixElement::virtual_corrections(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = virtual_corrections(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = virtual_corrections(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
double mur,
Particle const & WBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
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)
- {
+ double C2Lipatov(
+ CLHEP::HepLorentzVector const & qav,
+ CLHEP::HepLorentzVector const & qbv,
+ CLHEP::HepLorentzVector const & p1,
+ CLHEP::HepLorentzVector const & p2
+ ){
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(
- CLHEP::HepLorentzVector qav,
- CLHEP::HepLorentzVector qbv,
- CLHEP::HepLorentzVector p1,
- CLHEP::HepLorentzVector p2,
+ CLHEP::HepLorentzVector const & qav,
+ CLHEP::HepLorentzVector const & qbv,
+ CLHEP::HepLorentzVector const & p1,
+ CLHEP::HepLorentzVector const & p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
- else {
- double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
- return Cls-4./(kperp*kperp);
- }
+
+ 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
- {
+ double C2Lipatov( // B
+ CLHEP::HepLorentzVector const & qav,
+ CLHEP::HepLorentzVector const & qbv,
+ CLHEP::HepLorentzVector const & pim,
+ CLHEP::HepLorentzVector const & pip,
+ CLHEP::HepLorentzVector const & pom,
+ CLHEP::HepLorentzVector const & pop
+ ){
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(
- CLHEP::HepLorentzVector qav,
- CLHEP::HepLorentzVector qbv,
- CLHEP::HepLorentzVector pa,
- CLHEP::HepLorentzVector pb,
- CLHEP::HepLorentzVector p1,
- CLHEP::HepLorentzVector p2,
+ CLHEP::HepLorentzVector const & qav,
+ CLHEP::HepLorentzVector const & qbv,
+ CLHEP::HepLorentzVector const & pa,
+ CLHEP::HepLorentzVector const & pb,
+ CLHEP::HepLorentzVector const & p1,
+ CLHEP::HepLorentzVector const & p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
- else {
- double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
- double temp=Cls-4./(kperp*kperp);
- return temp;
- }
+
+ 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 ME_gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_qg(pn,pb,p1,pa);
else
return ME_qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return ME_qg(p1,pa,pn,pb);
else
return ME_qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_qQ(pn,pb,p1,pa);
else
return ME_qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return ME_qQbar(p1,pa,pn,pb);
else
return ME_qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_W_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// We know it cannot be gg incoming.
assert(!(aptype==21 && bptype==21));
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_W_qg(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qbarg(pn,plbar,pl,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return ME_W_qg(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qbarg(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 ME_W_qQ(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qQbar(pn,plbar,pl,pb,p1,pa);
}
else {
if (aptype>0)
return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa);
}
}
else{ // emission off a, (first argument paout)
if (aptype > 0) {
if (bptype > 0)
return ME_W_qQ(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qQbar(p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0)
return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qbarQbar(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 ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl);
}
else { // they are both quark
- if (wc==true) {// emission off b, i.e. b is first current
+ if (wc) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl);
}
else{
if (aptype>0)
return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl);
else //qqbar
return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
else //qbarqbar
return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
}
}
}
+ throw std::logic_error("unknown particle types");
}
/** 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
- ){
-
+ 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 ME_Wuno_qg(pn, pb, p1, pa, pg, plbar, pl);
else
return ME_Wuno_qbarg(pn, pb, p1, pa, pg, plbar, pl);
}
else { // they are both quark
- if (wc==true) {// emission off b, i.e. b is first current
+ if (wc) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
- return ME_Wuno_qQ(pn,pb,p1,pa,pg,plbar,pl);
+ return ME_Wuno_qQ(pn, pb, p1, pa, pg, plbar, pl);
else
- return ME_Wuno_qQbar(pn,pb,p1,pa,pg,plbar,pl);
+ return ME_Wuno_qQbar(pn, pb, p1, pa, pg, plbar, pl);
}
else{
if (aptype>0)
- return ME_Wuno_qbarQ(pn,pb,p1,pa,pg,plbar,pl);
+ return ME_Wuno_qbarQ(pn, pb, p1, pa, pg, plbar, pl);
else
- return ME_Wuno_qbarQbar(pn,pb,p1,pa,pg,plbar,pl);
+ return ME_Wuno_qbarQbar(pn, pb, p1, pa, pg, plbar, pl);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
- return ME_W_unof_qQ(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_W_unof_qQ(p1, pa, pn, pb, pg, plbar, pl);
// return ME_W_unof_qQ(pn,pb,p1,pa,pg,plbar,pl);
else //qqbar
- return ME_W_unof_qQbar(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_W_unof_qQbar(p1, pa, pn, pb, pg, plbar, pl);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
- return ME_W_unof_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_W_unof_qbarQ(p1, pa, pn, pb, pg, plbar, pl);
else //qbarqbar
- return ME_W_unof_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_W_unof_qbarQbar(p1, pa, pn, pb, pg, plbar, pl);
}
}
}
+ throw std::logic_error("unknown particle types");
}
/** \brief Matrix element squared for backward qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param 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
- ){
+ 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 ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn,pb)*CFbackward;}
- else {
- return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn,pb)*CFbackward;}
+ if (swapQuarkAntiquark)
+ return ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb)*CFbackward;
+ else
+ return ME_WExqqx_qbarqg(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 ME_WExqqx_qqbarQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
- else{
- return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
+ if (!wc){ // W Emitted from backwards qqx
+ if (swapQuarkAntiquark)
+ return ME_WExqqx_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb)*CFbackward;
+ else
+ return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;
}
else { // W Must be emitted from forwards leg.
if(bptype > 0){
- if (swapQuarkAntiquark){
- return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, false)*CFbackward;}
- else{
- return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, false)*CFbackward;}
+ if (swapQuarkAntiquark)
+ return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, false)*CFbackward;
+ else
+ return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, false)*CFbackward;
} else {
- if (swapQuarkAntiquark){
- return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, true)*CFbackward;}
- else{
- return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, true)*CFbackward;}
+ if (swapQuarkAntiquark)
+ return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, true)*CFbackward;
+ else
+ return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, true)*CFbackward;
}
}
}
- else{
- throw std::logic_error("Incompatible incoming particle types with qqxb");
- }
+ 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
- ){
+ 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 ME_WExqqx_qqbarg(pb, pqbar, plbar, pl, pq, p1,pa)*CFforward;}
- else {
- return ME_WExqqx_qbarqg(pb, pq, plbar, pl, pqbar, p1,pa)*CFforward;}
+ if (swapQuarkAntiquark)
+ return ME_WExqqx_qqbarg(pb, pqbar, plbar, pl, pq, p1, pa)*CFforward;
+ else
+ return ME_WExqqx_qbarqg(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 ME_WExqqx_qbarqQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
- else {
- return ME_WExqqx_qqbarQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
+ if (wc){ // W Emitted from forwards qqx
+ if (swapQuarkAntiquark)
+ return ME_WExqqx_qqbarQ(pb, pqbar, plbar, pl, pq, p1, pa)*CFforward;
+ else
+ return ME_WExqqx_qbarqQ(pb, pq, plbar, pl, pqbar, p1, pa)*CFforward;
}
// W Must be emitted from backwards leg.
if (aptype > 0){
- if (swapQuarkAntiquark){
- return ME_W_Exqqx_QQq(pa,pb, p1, pqbar, pq, plbar, pl, false)*CFforward;}
- else{
- return ME_W_Exqqx_QQq(pa,pb, p1, pq, pqbar, plbar, pl, false)*CFforward;}
- } else
- {
- if (swapQuarkAntiquark){
- return ME_W_Exqqx_QQq(pa,pb, p1, pqbar, pq, plbar, pl, true)*CFforward;}
- else{
- return ME_W_Exqqx_QQq(pa,pb, p1, pq, pqbar, plbar, pl, true)*CFforward;}
- }
- }
- else{
- throw std::logic_error("Incompatible incoming particle types with qqxf");
+ if (swapQuarkAntiquark)
+ return ME_W_Exqqx_QQq(pa, pb, p1, pqbar, pq, plbar, pl, false)*CFforward;
+ else
+ return ME_W_Exqqx_QQq(pa, pb, p1, pq, pqbar, plbar, pl, false)*CFforward;
+ } else {
+ if (swapQuarkAntiquark)
+ return ME_W_Exqqx_QQq(pa, pb, p1, pqbar, pq, plbar, pl, true)*CFforward;
+ else
+ return ME_W_Exqqx_QQq(pa, pb, p1, pq, pqbar, plbar, pl, true)*CFforward;
+ }
}
+ 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
- ){
+ int aptype, int bptype,
+ int nabove, int nbelow,
+ CLHEP::HepLorentzVector const & pa,
+ CLHEP::HepLorentzVector const & pb,
+ CLHEP::HepLorentzVector const & pq,
+ CLHEP::HepLorentzVector const & pqbar,
+ std::vector<HLV> const & partons,
+ CLHEP::HepLorentzVector const & plbar,
+ CLHEP::HepLorentzVector const & pl,
+ bool const wqq, bool const wc
+ ){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
bool swapQuarkAntiquark=false;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
}
double wt=1.;
if (aptype==21) wt*=K_g(partons.front(),pa)/HEJ::C_F;
if (bptype==21) wt*=K_g(partons.back(),pb)/HEJ::C_F;
if (aptype <=0 && bptype <=0){ // Both External AntiQuark
if (wqq==1){//emission from central qqbar
- return wt*ME_WCenqqx_qq(pa, pb, pl,plbar, partons,true,true,
+ return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons,true,true,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
- return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, true,true,
+ return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, true,true,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(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*ME_WCenqqx_qq(pa, pb, pl,plbar, partons, false, true,
+ return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons, false, true,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
- return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons,false,true,
+ return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons,false,true,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
- return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, false, true,
+ return wt*ME_W_Cenqqx_qq(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*ME_WCenqqx_qq(pa, pb, pl,plbar, partons, true, false,
+ return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons, true, false,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
- return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, true, false,
+ return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, true, false,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
- return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, true, false,
+ return wt*ME_W_Cenqqx_qq(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*ME_WCenqqx_qq(pa, pb, pl, plbar, partons, false, false,
swapQuarkAntiquark, nabove);}
else if (wc==1){//emission from b leg
- return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, false, false,
+ return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, false, false,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
- return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, false, false,
+ return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, false, false,
swapQuarkAntiquark, nabove, nbelow, false);
}
-
}
+ throw std::logic_error("unknown particle types");
}
/** \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 ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_H_qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
else
return ME_H_qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return ME_H_qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
else
return ME_H_qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return ME_H_qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return ME_H_qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return ME_H_qbarQbar(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 ME_H_unof_qg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unof_qbarg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_H_unof_qQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unof_qQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (aptype>0)
return ME_H_unof_qbarQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unof_qbarQbar(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 ME_H_unob_gQ(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unob_gQbar(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return ME_H_unob_qQ(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unob_qbarQ(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (bptype>0)
return ME_H_unob_qQbar(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unob_qbarQbar(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_resummable(ev.type())) return 0.;
auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing()))
return tree_kin_jets(ev);
switch(AWZH_boson->type){
case pid::Higgs:
return tree_kin_Higgs(ev);
case pid::Wp:
case pid::Wm:
return tree_kin_W(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
namespace{
constexpr int extremal_jet_idx = 1;
constexpr int no_extremal_jet_idx = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == extremal_jet_idx;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
double lambda
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
Event const & ev
) const{
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
// 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 {
+ 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
- ) {
+ 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;
+ const auto begin_ladder = cbegin(partons) + 1;
+ const auto end_ladder = cend(partons) - 1;
- bool wc = true;
+ bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
- if (aptype!=partons[0].type) { //leg a emits w
- wc = false;
- q0 -=pl + plbar;
- }
+ if(!wc)
+ 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;
+ const auto begin_ladder = cbegin(partons) + 2;
+ const auto end_ladder = cend(partons) - 1;
- bool wc = true;
+ bool wc = aptype==partons[1].type; //leg b emits w
auto q0 = pa - p1 -pg;
- if (aptype!=partons[1].type) { //leg a emits w
- wc = false;
- q0 -=pl + plbar;
- }
+ if(!wc)
+ 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;
+ const auto begin_ladder = cbegin(partons) + 1;
+ const auto end_ladder = cend(partons) - 2;
- bool wc = true;
+ bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
- if (aptype!=partons[0].type) { //leg a emits w
- wc = false;
- q0 -=pl + plbar;
- }
+ if(!wc)
+ 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;
+ const auto begin_ladder = cbegin(partons) + 2;
+ const auto end_ladder = cend(partons) - 1;
- bool wc = true;
+ bool wc = bptype!=partons.back().type; //leg b emits w
auto q0 = pa - pq - pqbar;
- if (partons[1].type!=partons[0].type) { //leg a emits w
- wc = false;
- q0 -=pl + plbar;
- }
+ if(!wc)
+ 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;
+ const auto begin_ladder = cbegin(partons) + 1;
+ const auto end_ladder = cend(partons) - 2;
- bool wc = true;
+ bool wc = aptype==partons.front().type; //leg b emits w
auto q0 = pa - p1;
- if (aptype!=partons[0].type) { //leg a emits w
- wc = false;
- q0 -=pl + plbar;
- }
+ if(!wc)
+ 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;
+ const auto begin_ladder = cbegin(partons) + 1;
+ const auto end_ladder_1 = (backmidquark);
+ const auto begin_ladder_2 = (backmidquark+2);
+ const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
int nabove = std::distance(begin_ladder, backmidquark);
int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc
);
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 {
+ 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,
+ CLHEP::HepLorentzVector const & p1out,
+ CLHEP::HepLorentzVector const & p1in,
ParticleID type2,
- CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
- CLHEP::HepLorentzVector pH,
+ CLHEP::HepLorentzVector const & p2out,
+ CLHEP::HepLorentzVector const & p2in,
+ CLHEP::HepLorentzVector const & pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
// 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*ME_Houtside_gq(
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 {
+ 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 {
+ 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 {
+ 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;
}
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{
+ double MatrixElement::tree_param(Event const & ev, double mur) const {
assert(is_resummable(ev.type()));
const auto begin_partons = ev.begin_partons();
const auto end_partons = ev.end_partons();
const auto num_partons = std::distance(begin_partons, end_partons);
const double alpha_s = alpha_s_(mur);
const double gs2 = 4.*M_PI*alpha_s;
double res = std::pow(gs2, num_partons);
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(num_partons >= 2);
const auto first_emission = std::next(begin_partons);
const auto last_emission = std::prev(end_partons);
for(auto parton = first_emission; parton != last_emission; ++parton){
res *= 1. + alpha_s/(2.*M_PI)*beta0*log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s)*res;
}
} // namespace HEJ
diff --git a/src/Wjets.cc b/src/Wjets.cc
index 89511e6..7f53581 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,1295 +1,1295 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/currents.hh"
#include "HEJ/utility.hh"
#include "HEJ/Tensor.hh"
#include "HEJ/Constants.hh"
#include <array>
#include <iostream>
namespace { // Helper Functions
// FKL W Helper Functions
double WProp (const HLV & plbar, const HLV & pl){
COM propW = COM(0.,-1.)/( (pl+plbar).m2() - HEJ::MW*HEJ::MW + COM(0.,1.)*HEJ::MW*HEJ::GammaW);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
CCurrent jW (HLV pout, bool helout, HLV plbar, bool hellbar, HLV pl, bool hell,
HLV pin, bool helin
){
COM cur[4];
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
CCurrent sum(0.,0.,0.,0.);
// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is
// anti-lepton(5)
// Need to swap e and nu for events with W- --> e- nubar!
if (helin==helout && hellbar==hell) {
HLV qa=pout+plbar+pl;
HLV qb=pin-plbar-pl;
double ta(qa.m2()),tb(qb.m2());
CCurrent temp2,temp3,temp5;
CCurrent t65 = joo(pl,hell,plbar,hellbar);
CCurrent vout(pout.e(),pout.x(),pout.y(),pout.z());
CCurrent vin(pin.e(),pin.x(),pin.y(),pin.z());
COM brac615=t65.dot(vout);
COM brac645=t65.dot(vin);
// prod1565 and prod6465 are zero for Ws (not Zs)!!
temp2 = joo(pout,helout,pl,helout);
COM prod1665=temp2.dot(t65);
temp3 = joi(plbar,helin,pin,helin);
COM prod5465=temp3.dot(t65);
temp2=joo(pout,helout,plbar,helout);
temp3=joi(pl,hell,pin,helin);
temp5=joi(pout,helout,pin,helin);
CCurrent term1,term2,term3;
term1=(2.*brac615/ta+2.*brac645/tb)*temp5;
term2=(prod1665/ta)*temp3;
term3=(-prod5465/tb)*temp2;
sum=term1+term2+term3;
}
return sum;
}
CCurrent jWbar (HLV pout, bool helout, HLV plbar, bool hellbar, HLV pl, bool hell,
HLV pin, bool helin
){
COM cur[4];
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
CCurrent sum(0.,0.,0.,0.);
// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is
// anti-lepton(5)
// Need to swap e and nu for events with W- --> e- nubar!
if (helin==helout && hellbar==hell) {
HLV qa=pout+plbar+pl;
HLV qb=pin-plbar-pl;
double ta(qa.m2()),tb(qb.m2());
CCurrent temp2,temp3,temp5;
CCurrent t65 = joo(pl,hell,plbar,hellbar);
CCurrent vout(pout.e(),pout.x(),pout.y(),pout.z());
CCurrent vin(pin.e(),pin.x(),pin.y(),pin.z());
COM brac615=t65.dot(vout);
COM brac645=t65.dot(vin);
// prod1565 and prod6465 are zero for Ws (not Zs)!!
temp2 = joo(plbar,helout,pout,helout); // temp2 is <5|alpha|1>
COM prod5165=temp2.dot(t65);
temp3 = jio(pin,helin,pl,helin); // temp3 is <4|alpha|6>
COM prod4665=temp3.dot(t65);
temp2=joo(pl,helout,pout,helout); // temp2 is now <6|mu|1>
temp3=jio(pin,helin,plbar,helin); // temp3 is now <4|mu|5>
temp5=jio(pin,helin,pout,helout); // temp5 is <4|mu|1>
CCurrent term1,term2,term3;
term1 =(-2.*brac615/ta-2.*brac645/tb)*temp5;
term2 =(-prod5165/ta)*temp3;
term3 =(prod4665/tb)*temp2;
sum = term1 + term2 + term3;
}
return sum;
}
// Extremal quark current with W emission.
// Using Tensor class rather than CCurrent
Tensor <1> jW(HLV pin, HLV pout, HLV plbar, HLV pl, bool aqline){
// Build the external quark line W Emmision
Tensor<1> ABCurr = TCurrent(pl, false, plbar, false);
Tensor<1> Tp4W = Construct1Tensor((pout+pl+plbar));//p4+pw
Tensor<1> TpbW = Construct1Tensor((pin-pl-plbar));//pb-pw
Tensor<3> J4bBlank;
if (aqline){
J4bBlank = T3Current(pin,false,pout,false);
}
else{
J4bBlank = T3Current(pout,false,pin,false);
}
double t4AB = (pout+pl+plbar).m2();
double tbAB = (pin-pl-plbar).m2();
Tensor<2> J4b1 = (J4bBlank.contract(Tp4W,2))/t4AB;
Tensor<2> J4b2 = (J4bBlank.contract(TpbW,2))/tbAB;
Tensor<2> T4bmMom(0.);
if (aqline){
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
T4bmMom(mu, nu) = (J4b1(nu,mu) + J4b2(mu,nu))*COM(0,-1);
}
}
}
else{
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
T4bmMom(nu,mu) = (J4b1(nu,mu) + J4b2(mu,nu))*COM(0,1);
}
}
}
Tensor<1> T4bm = T4bmMom.contract(ABCurr,1);
return T4bm;
}
// Relevant W+Jets Unordered Contribution Helper Functions
double jM2Wuno(HLV pg, HLV p1,HLV plbar, HLV pl, HLV pa, bool h1,
HLV p2, HLV pb, bool h2, bool pol
){
//@TODO Simplify the below (less Tensor class?)
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;
}
HLV pW = pl+plbar;
HLV q1g=pa-pW-p1-pg;
HLV q1 = pa-p1-pW;
HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double taW1 = (pa-pW-p1).m2();
const double tb2 = (pb-p2).m2();
const double tb2g = (pb-p2-pg).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+pg).m2();
const double tag = (pa-pg).m2();
const double taWg = (pa-pW-pg).m2();
//use p1 as ref vec in pol tensor
Tensor<1> epsg = eps(pg,p2,pol);
Tensor<1> epsW = TCurrent(pl,false,plbar,false);
Tensor<1> j2b = TCurrent(p2,h2,pb,h2);
Tensor<1> Tq1q2 = Construct1Tensor((q1+q2)/taW1 + (pb/pb.dot(pg)
+p2/p2.dot(pg)) * tb2/(2*tb2g));
Tensor<1> Tq1g = Construct1Tensor((-pg-q1))/taW1;
Tensor<1> Tq2g = Construct1Tensor((pg-q2));
Tensor<1> TqaW = Construct1Tensor((pa-pW));//pa-pw
Tensor<1> Tqag = Construct1Tensor((pa-pg));
Tensor<1> TqaWg = Construct1Tensor((pa-pg-pW));
Tensor<1> Tp1g = Construct1Tensor((p1+pg));
Tensor<1> Tp1W = Construct1Tensor((p1+pW));//p1+pw
Tensor<1> Tp1gW = Construct1Tensor((p1+pg+pW));//p1+pw+pg
Tensor<2> g=Metric();
Tensor<3> J31a = T3Current(p1, h1, pa, h1);
Tensor<2> J2_qaW =J31a.contract(TqaW/taW, 2);
Tensor<2> J2_p1W =J31a.contract(Tp1W/s1W, 2);
Tensor<3> L1a = outer(Tq1q2, J2_qaW);
Tensor<3> L1b = outer(Tq1q2, J2_p1W);
Tensor<3> L2a = outer(Tq1g,J2_qaW);
Tensor<3> L2b = outer(Tq1g, J2_p1W);
Tensor<3> L3 = outer(g, J2_qaW.contract(Tq2g,1)+J2_p1W.contract(Tq2g,2))/taW1;
Tensor<3> L(0.);
Tensor<5> J51a = T5Current(p1, h1, pa, h1);
Tensor<4> J_qaW = J51a.contract(TqaW,4);
Tensor<4> J_qag = J51a.contract(Tqag,4);
Tensor<4> J_p1gW = J51a.contract(Tp1gW,4);
Tensor<3> U1a = J_qaW.contract(Tp1g,2);
Tensor<3> U1b = J_p1gW.contract(Tp1g,2);
Tensor<3> U1c = J_p1gW.contract(Tp1W,2);
Tensor<3> U1(0.);
Tensor<3> U2a = J_qaW.contract(TqaWg,2);
Tensor<3> U2b = J_qag.contract(TqaWg,2);
Tensor<3> U2c = J_qag.contract(Tp1W,2);
Tensor<3> U2(0.);
for(int nu=0; nu<4;nu++){
for(int mu=0;mu<4;mu++){
for(int rho=0;rho<4;rho++){
L(nu, mu, rho) = L1a(nu,mu,rho) + L1b(nu,rho,mu)
+ L2a(mu,nu,rho) + L2b(mu,rho,nu) + L3(mu,nu,rho);
U1(nu, mu, rho) = U1a(nu, mu, rho) / (s1g*taW)
+ U1b(nu,rho,mu) / (s1g*s1gW) + U1c(rho,nu,mu) / (s1W*s1gW);
U2(nu,mu,rho) = U2a(mu,nu,rho) / (taWg*taW)
+ U2b(mu,rho,nu) / (taWg*tag) + U2c(rho,mu,nu) / (s1W*tag);
}
}
}
COM X = ((((U1-L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
COM Y = ((((U2+L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
double amp = HEJ::C_A*HEJ::C_F*HEJ::C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
double t1 = q1g.m2();
double t2 = q2.m2();
double WPropfact = WProp(plbar, pl);
//Divide by WProp
amp*=WPropfact;
//Divide by t-channels
amp/=(t1*t2);
//Average over initial states
amp/=(4.*HEJ::C_A*HEJ::C_A);
return amp;
}
// Relevant Wqqx Helper Functions.
//g->qxqlxl (Calculates gluon to qqx Current. See JV_\mu in WSubleading Notes)
Tensor <1> gtqqxW(HLV pq,HLV pqbar,HLV pl,HLV plbar){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
Tensor<1> Tpq = Construct1Tensor(pq);
Tensor<1> Tpqbar = Construct1Tensor(pqbar);
Tensor<1> TAB = Construct1Tensor(pl+plbar);
// Define llx current.
Tensor<1> ABCur = TCurrent(pl, false, plbar, false);
//blank 3 Gamma Current
Tensor<3> JV23 = T3Current(pq,false,pqbar,false);
// Components of g->qqW before W Contraction
Tensor<2> JV1 = JV23.contract((Tpq + TAB),2)/(s2AB);
Tensor<2> JV2 = JV23.contract((Tpqbar + TAB),2)/(s3AB);
// g->qqW Current. Note Minus between terms due to momentum flow.
// Also note: (-I)^2 from W vert. (I) from Quark prop.
Tensor<1> JVCur = (JV1.contract(ABCur,1) - JV2.contract(ABCur,2))*COM(0.,-1.);
return JVCur;
}
// Helper Functions Calculate the Crossed Contribution
Tensor <2> MCrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar, HLV pl,
HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCross?
// Useful propagator factors
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV q1, q3;
q1=pa;
for(int i=0; i<nabove+1;i++){
q1=q1-partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double tcro1=(q3+pq).m2();
double tcro2=(q1-pqbar).m2();
Tensor<1> Tpq = Construct1Tensor(pq);
Tensor<1> Tpqbar = Construct1Tensor(pqbar);
Tensor<1> TAB = Construct1Tensor(pl+plbar);
Tensor<1> Tq1 = Construct1Tensor(q1);
Tensor<1> Tq3 = Construct1Tensor(q3);
// Define llx current.
Tensor<1> ABCur = TCurrent(pl, false, plbar,false);
//Blank 5 gamma Current
Tensor<5> J523 = T5Current(pq,false,pqbar,false);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_q3q = J523.contract((Tq3+Tpq),2);
Tensor<4> J_2AB = J523.contract((Tpq+TAB),2);
// Components of Crossed Vertex Contribution
Tensor<3> Xcro1 = J_q3q.contract((Tpqbar + TAB),3);
Tensor<3> Xcro2 = J_q3q.contract((Tq1-Tpqbar),3);
Tensor<3> Xcro3 = J_2AB.contract((Tq1-Tpqbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xcro1Cont = Xcro1.contract(ABCur,3)/(tcro1*s3AB);
Tensor<2> Xcro2Cont = Xcro2.contract(ABCur,2)/(tcro1*tcro2);
Tensor<2> Xcro3Cont = Xcro3.contract(ABCur,1)/(s2AB*tcro2);
//Initialise the Crossed Vertex Object
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro(mu,nu) = -(-Xcro1Cont(nu,mu)+Xcro2Cont(nu,mu)+Xcro3Cont(nu,mu));
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar,
HLV pl, HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncross?
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV q1, q3;
q1=pa;
for(int i=0; i<nabove+1;i++){
q1=q1-partons.at(i);
}
q3 = q1 - pl - plbar - pq - pqbar;
double tunc1 = (q1-pq).m2();
double tunc2 = (q3+pqbar).m2();
Tensor<1> Tpq = Construct1Tensor(pq);
Tensor<1> Tpqbar = Construct1Tensor(pqbar);
Tensor<1> TAB = Construct1Tensor(pl+plbar);
Tensor<1> Tq1 = Construct1Tensor(q1);
Tensor<1> Tq3 = Construct1Tensor(q3);
// Define llx current.
Tensor<1> ABCur = TCurrent(pl, false, plbar, false);
//Blank 5 gamma Current
Tensor<5> J523 = T5Current(pq,false,pqbar,false);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_2AB = J523.contract((Tpq+TAB),2);
Tensor<4> J_q1q = J523.contract((Tq1-Tpq),2);
// 2 Contractions taken care of.
Tensor<3> Xunc1 = J_2AB.contract((Tq3+Tpqbar),3);
Tensor<3> Xunc2 = J_q1q.contract((Tq3+Tpqbar),3);
Tensor<3> Xunc3 = J_q1q.contract((Tpqbar+TAB),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xunc1Cont = Xunc1.contract(ABCur,1)/(s2AB*tunc2);
Tensor<2> Xunc2Cont = Xunc2.contract(ABCur,2)/(tunc1*tunc2);
Tensor<2> Xunc3Cont = Xunc3.contract(ABCur,3)/(tunc1*s3AB);
//Initialise the Uncrossed Vertex Object
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc(mu,nu) = -(- Xunc1Cont(mu,nu)+Xunc2Cont(mu,nu) +Xunc3Cont(mu,nu));
}
}
return Xunc;
}
// Helper Functions Calculate the g->qqxW (Eikonal) Contributions
Tensor <2> MSymW(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
HLV pl,HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MSym?
double sa2=(pa+pq).m2();
double s12=(p1+pq).m2();
double sa3=(pa+pqbar).m2();
double s13=(p1+pqbar).m2();
double saA=(pa+pl).m2();
double s1A=(p1+pl).m2();
double saB=(pa+plbar).m2();
double s1B=(p1+plbar).m2();
double sb2=(pb+pq).m2();
double s42=(p4+pq).m2();
double sb3=(pb+pqbar).m2();
double s43=(p4+pqbar).m2();
double sbA=(pb+pl).m2();
double s4A=(p4+pl).m2();
double sbB=(pb+plbar).m2();
double s4B=(p4+plbar).m2();
double s23AB=(pl+plbar+pq+pqbar).m2();
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3=q1-pq-pqbar-plbar-pl;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Define Tensors to be used
Tensor<1> Tp1 = Construct1Tensor(p1);
Tensor<1> Tp4 = Construct1Tensor(p4);
Tensor<1> Tpa = Construct1Tensor(pa);
Tensor<1> Tpb = Construct1Tensor(pb);
Tensor<1> Tpq = Construct1Tensor(pq);
Tensor<1> Tpqbar = Construct1Tensor(pqbar);
Tensor<1> TAB = Construct1Tensor(pl+plbar);
Tensor<1> Tq1 = Construct1Tensor(q1);
Tensor<1> Tq3 = Construct1Tensor(q3);
Tensor<2> g=Metric();
// g->qqW Current (Factors of sqrt2 dealt with in this function.)
Tensor<1> JV = gtqqxW(pq,pqbar,pl,plbar);
// 1a gluon emisson Contribution
Tensor<3> X1a = outer(g, Tp1*(t1/(s12+s13+s1A+s1B))
+ Tpa*(t1/(sa2+sa3+saA+saB)) );
Tensor<2> X1aCont = X1a.contract(JV,3);
//4b gluon emission Contribution
Tensor<3> X4b = outer(g, Tp4*(t3/(s42+s43+s4A+s4B))
+ Tpb*(t3/(sb2+sb3+sbA+sbB)) );
Tensor<2> X4bCont = X4b.contract(JV,3);
//Set up each term of 3G diagram.
Tensor<3> X3g1 = outer(Tq1+Tpq+Tpqbar+TAB, g);
Tensor<3> X3g2 = outer(Tq3-Tpq-Tpqbar-TAB, g);
Tensor<3> X3g3 = outer(Tq1+Tq3, g);
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(JV,3);
Tensor<2> X3g2Cont = X3g2.contract(JV,2);
Tensor<2> X3g3Cont = X3g3.contract(JV,1);
// XSym is an amalgamation of x1a, X4b and X3g.
// Makes sense from a colour factor point of view.
Tensor<2>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym(mu,nu) = (X3g1Cont(nu,mu) + X3g2Cont(mu,nu) - X3g3Cont(nu,mu))
+ (X1aCont(mu,nu) - X4bCont(mu,nu));
}
}
return Xsym/s23AB;
}
Tensor <2> MCross(HLV pa, HLV pq, HLV pqbar, std::vector<HLV> partons,
bool hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2=(q1-pqbar).m2();
Tensor<1> Tq1 = Construct1Tensor(q1-pqbar);
//Blank 3 gamma Current
Tensor<3> J323 = T3Current(pq,hq,pqbar,hq);
// 2 gamma current (with 1 contraction already).
Tensor<2> XCroCont = J323.contract((Tq1),2)/(t2);
//Initialise the Crossed Vertex
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro(mu,nu) = XCroCont(nu,mu);
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncross(HLV pa, HLV pq,HLV pqbar, std::vector<HLV> partons,
bool hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2 = (q1-pq).m2();
Tensor<1> Tq1 = Construct1Tensor(q1-pq);
//Blank 3 gamma Current
Tensor<3> J323 = T3Current(pq,hq,pqbar,hq);
// 2 gamma currents (with 1 contraction already).
Tensor<2> XUncCont = J323.contract((Tq1),2)/t2;
//Initialise the Uncrossed Vertex
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc(mu,nu) = -XUncCont(mu,nu);
}
}
return Xunc;
}
// Helper Functions Calculate the Eikonal Contributions
Tensor <2> MSym(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
std::vector<HLV> partons, bool hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MsymW?
HLV q1, q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1-pq-pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
double s23 = (pq+pqbar).m2();
double sa2 = (pa+pq).m2();
double sa3 = (pa+pqbar).m2();
double s12 = (p1+pq).m2();
double s13 = (p1+pqbar).m2();
double sb2 = (pb+pq).m2();
double sb3 = (pb+pqbar).m2();
double s42 = (p4+pq).m2();
double s43 = (p4+pqbar).m2();
//Define Tensors to be used
Tensor<1> Tp1 = Construct1Tensor(p1);
Tensor<1> Tp4 = Construct1Tensor(p4);
Tensor<1> Tpa = Construct1Tensor(pa);
Tensor<1> Tpb = Construct1Tensor(pb);
Tensor<1> Tpq = Construct1Tensor(pq);
Tensor<1> Tpqbar = Construct1Tensor(pqbar);
Tensor<1> Tq1 = Construct1Tensor(q1);
Tensor<1> Tq3 = Construct1Tensor(q3);
Tensor<2> g=Metric();
Tensor<1> qqxCur = TCurrent(pq, hq, pqbar, hq);
// // 1a gluon emisson Contribution
Tensor<3> X1a = outer(g, Tp1*(t1/(s12+s13))+Tpa*(t1/(sa2+sa3)));
Tensor<2> X1aCont = X1a.contract(qqxCur,3);
// //4b gluon emission Contribution
Tensor<3> X4b = outer(g, Tp4*(t3/(s42+s43)) + Tpb*(t3/(sb2+sb3)));
Tensor<2> X4bCont = X4b.contract(qqxCur,3);
// New Formulation Corresponding to New Analytics
Tensor<3> X3g1 = outer(Tq1+Tpq+Tpqbar, g);
Tensor<3> X3g2 = outer(Tq3-Tpq-Tpqbar, g);
Tensor<3> X3g3 = outer(Tq1+Tq3, g);
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(qqxCur,3);
Tensor<2> X3g2Cont = X3g2.contract(qqxCur,2);
Tensor<2> X3g3Cont = X3g3.contract(qqxCur,1);
Tensor<2>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym(mu, nu) = COM(0,1) * ( (X3g1Cont(nu,mu) + X3g2Cont(mu,nu)
- X3g3Cont(nu,mu)) + (X1aCont(mu,nu) - X4bCont(mu,nu)) );
}
}
return Xsym/s23;
}
} // Anonymous Namespace helper functions
//! W+Jets FKL Contributions
/**
* @brief W+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_\mu.
* Handles all possible incoming states.
*/
double jW_j( HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
bool aqlineb, bool aqlinef
){
CCurrent mj1m,mj2p,mj2m;
HLV q1=p1in-p1out-plbar-pl;
HLV q2=-(p2in-p2out);
if(aqlineb) mj1m=jWbar(p1out,false,plbar,false,pl,false,p1in,false);
else mj1m=jW(p1out,false,plbar,false,pl,false,p1in,false);
if(aqlinef){
mj2p=jio(p2in,true,p2out,true);
mj2m=jio(p2in,false,p2out,false);
} else{
mj2p=joi(p2out,true,p2in,true);
mj2m=joi(p2out,false,p2in,false);
}
COM Mmp=mj1m.dot(mj2p);
COM Mmm=mj1m.dot(mj2m);
// sum of spinor strings ||^2
double a2Mmp=abs2(Mmp);
double a2Mmm=abs2(Mmm);
double WPropfact = WProp(plbar, pl);
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
double ME_W_qQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false);
}
double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true);
}
double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false);
}
double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true);
}
double ME_W_qg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false)*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false)*K_g(p2out, p2in)/HEJ::C_F;
}
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param pg Unordered Gluon momenta
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_{uno}_\mu. Ie, unordered with W emission opposite side.
* Handles all possible incoming states.
*/
double jW_juno(HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out,
HLV p2in, HLV pg, bool aqlineb, bool aqlinef){
CCurrent mj1m,mj2p,mj2m, jgbm,jgbp,j2gm,j2gp;
HLV q1=p1in-p1out-plbar-pl;
HLV q2=-(p2in-p2out-pg);
HLV q3=-(p2in-p2out);
if(aqlineb) mj1m=jWbar(p1out,false,plbar,false,pl,false,p1in,false);
else mj1m=jW(p1out,false,plbar,false,pl,false,p1in,false);
//@TODO Is aqlinef necessary? Gives same results.
if(aqlinef){
mj2p=jio(p2in,true,p2out,true);
mj2m=jio(p2in,false,p2out,false);
j2gp=joo(pg,true,p2out,true);
j2gm=joo(pg,false,p2out,false);
jgbp=jio(p2in,true,pg,true);
jgbm=jio(p2in,false,pg,false);
} else{
mj2p=joi(p2out,true,p2in,true);
mj2m=joi(p2out,false,p2in,false);
j2gp=joo(p2out,true,pg,true);
j2gm=joo(p2out,false,pg,false);
jgbp=joi(pg,true,p2in,true);
jgbm=joi(pg,false,p2in,false);
}
// Dot products of these which occur again and again
COM MWmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM MWmm=mj1m.dot(mj2m);
CCurrent qsum(q2+q3);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
CCurrent p2o(p2out);
CCurrent p2i(p2in);
Lmm=( (-1.)*qsum*(MWmm) + (-2.*mj1m.dot(pg))*mj2m + 2.*mj2m.dot(pg)*mj1m
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmm/2. ) )/q3.m2();
Lmp=( (-1.)*qsum*(MWmp) + (-2.*mj1m.dot(pg))*mj2p + 2.*mj2p.dot(pg)*mj1m
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmp/2. ) )/q3.m2();
U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
double amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
ampsq*=WProp(plbar, pl);
return ampsq/((16)*(q2.m2()*q1.m2()));
}
double ME_W_unob_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false);
}
double ME_W_unob_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true);
}
double ME_W_unob_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false);
}
double ME_W_unob_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true);
}
double ME_W_unof_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, false, false);
}
double ME_W_unof_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
- return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, true, false);
+ return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, false, true);
}
double ME_W_unof_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
- return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, false, true);
+ return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, true, false);
}
double ME_W_unof_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, true, true);
}
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (Quark - W and Uno emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (Quark - W and Uno emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
*
* Calculates j_W_{uno} ^\mu j_\mu. Ie, unordered with W emission same side.
* Handles all possible incoming states. Note this handles both forward and back-
* -ward Wuno emission. For forward, ensure p1out is the uno and W emission parton.
* @TODO: Include separate wrapper functions for forward and backward to clean up
* ME_W_unof_current in `MatrixElement.cc`.
*/
double jWuno_j(HLV pg, HLV p1out, HLV plbar, HLV pl, HLV p1in,
HLV p2out, HLV p2in, bool aqlineb
){
//Calculate different Helicity choices
double ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,aqlineb,p2out,p2in,true,true);
double ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,aqlineb,p2out,p2in,true,false);
double ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,aqlineb,p2out,p2in,false,true);
double ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,aqlineb,p2out,p2in,false,false);
return ME2mpp + ME2mpm + ME2mmp + ME2mmm;
}
double ME_Wuno_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false);
}
double ME_Wuno_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false);
}
double ME_Wuno_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true);
}
double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true);
}
double ME_Wuno_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false)*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_Wuno_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true)*K_g(p2out, p2in)/HEJ::C_F;
}
/**
* @brief W+Jets Extremal qqx Contributions, function to handle all incoming types.
* @param pgin Incoming gluon which will split into qqx.
* @param pqout Quark of extremal qqx outgoing (W-Emission).
* @param plbar Outgoing anti-lepton momenta
* @param pl Outgoing lepton momenta
* @param pqbarout Anti-quark of extremal qqx pair. (W-Emission)
* @param pout Outgoing Particle 2 (end of FKL chain)
* @param p2in Incoming Particle 2
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W_{qqx} ^\mu j_\mu. Ie, Ex-QQX with W emission same side.
* Handles all possible incoming states. Calculated via crossing symmetry from jWuno_j.
*/
-double jWqqx_j(HLV pgin, HLV pqout,HLV plbar,HLV pl,
+double jWqqx_j(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in, bool aqlinef){
//Calculate Different Helicity Configurations.
double ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,aqlinef,p2out,p2in,true,true);
double ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,aqlinef,p2out,p2in,true,false);
double ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,aqlinef,p2out,p2in,false,true);
double ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,aqlinef,p2out,p2in,false,false);
//Helicity sum
double ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
-double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout,HLV plbar,HLV pl,
+double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false);
}
-double ME_WExqqx_qqbarQ(HLV pgin, HLV pqbarout,HLV plbar,HLV pl,
+double ME_WExqqx_qqbarQ(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true);
}
-double ME_WExqqx_qbarqg(HLV pgin, HLV pqout,HLV plbar,HLV pl,
+double ME_WExqqx_qbarqg(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false)*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_WExqqx_qqbarg(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true)*K_g(p2out,p2in)/HEJ::C_F;
}
namespace {
//Function to calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm1(HLV pb, HLV p2, HLV p3, bool hel2, bool helg, HLV refmom){
//@TODO Simplify the calculation below. (Less Tensor class use?)
double t1 = (p3-pb)*(p3-pb);
Tensor<1> Tp3 = Construct1Tensor((p3));//p3
Tensor<1> Tpb = Construct1Tensor((pb));//pb
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank = T3Current(p2,hel2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(Tp3-Tpb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,2)/t1;
return gqqCur*(-1);
}
//Function to calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm2(HLV pb, HLV p2, HLV p3, bool hel2, bool helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double t1 = (p2-pb)*(p2-pb);
Tensor<1> Tp2 = Construct1Tensor((p2));//p2
Tensor<1> Tpb = Construct1Tensor((pb));//pb
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb,refmom, helg);
Tensor<3> qqCurBlank = T3Current(p2,hel2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(Tp2-Tpb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,1)/t1;
return gqqCur;
}
//Function to calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm3(HLV pb, HLV p2, HLV p3, bool hel2, bool helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s23 = (p2+p3)*(p2+p3);
Tensor<1> Tp2 = Construct1Tensor((p2));//p2
Tensor<1> Tp3 = Construct1Tensor((p3));//p3
Tensor<1> Tpb = Construct1Tensor((pb));//pb
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<2> g=Metric();
Tensor<3> qqCurBlank1 = outer(Tp2+Tp3, g)/s23;
Tensor<3> qqCurBlank2 = outer(Tpb, g)/s23;
Tensor<1> Cur23 = TCurrent(p2,hel2, p3,hel2);
Tensor<2> qqCur1 = qqCurBlank1.contract(Cur23,3);
Tensor<2> qqCur2 = qqCurBlank2.contract(Cur23,3);
Tensor<2> qqCur3 = qqCurBlank2.contract(Cur23,1);
Tensor<1> gqqCur = (qqCur1.contract(epsg,1)
- qqCur2.contract(epsg,2)
+ qqCur3.contract(epsg,1))*2*COM(0,1);
return gqqCur;
}
}
// no wqq emission
double ME_W_Exqqx_QQq(HLV pa, HLV pb, HLV p1, HLV p2,
HLV p3,HLV plbar, HLV pl, bool aqlinepa
){
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;}
// 2 independent helicity choices (complex conjugation related).
Tensor<1> TMmmm1 = qggm1(pb,p2,p3,false,false, pa);
Tensor<1> TMmmm2 = qggm2(pb,p2,p3,false,false, pa);
Tensor<1> TMmmm3 = qggm3(pb,p2,p3,false,false, pa);
Tensor<1> TMpmm1 = qggm1(pb,p2,p3,false,true, pa);
Tensor<1> TMpmm2 = qggm2(pb,p2,p3,false,true, pa);
Tensor<1> TMpmm3 = qggm3(pb,p2,p3,false,true, pa);
// Build the external quark line W Emmision
Tensor<1> cur1a = jW(pa,p1,plbar,pl, aqlinepa);
//Contract with the qqxCurrent.
COM Mmmm1 = TMmmm1.contract(cur1a,1);
COM Mmmm2 = TMmmm2.contract(cur1a,1);
COM Mmmm3 = TMmmm3.contract(cur1a,1);
COM Mpmm1 = TMpmm1.contract(cur1a,1);
COM Mpmm2 = TMpmm2.contract(cur1a,1);
COM Mpmm3 = TMpmm3.contract(cur1a,1);
//Colour factors:
COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
cm1m1=8./3.;
cm2m2=8./3.;
cm3m3=6.;
cm1m2 =-1./3.;
cm1m3 = -3.*COM(0.,1.);
cm2m3 = 3.*COM(0.,1.);
//Sqaure and sum for each helicity config:
double Mmmm = real( cm1m1*pow(abs(Mmmm1),2) + cm2m2*pow(abs(Mmmm2),2)
+ cm3m3*pow(abs(Mmmm3),2) + 2.*real(cm1m2*Mmmm1*conj(Mmmm2))
+ 2.*real(cm1m3*Mmmm1*conj(Mmmm3))
+ 2.*real(cm2m3*Mmmm2*conj(Mmmm3)) );
double Mpmm = real( cm1m1*pow(abs(Mpmm1),2) + cm2m2*pow(abs(Mpmm2),2)
+ cm3m3*pow(abs(Mpmm3),2) + 2.*real(cm1m2*Mpmm1*conj(Mpmm2))
+ 2.*real(cm1m3*Mpmm1*conj(Mpmm3))
+ 2.*real(cm2m3*Mpmm2*conj(Mpmm3)) );
// Divide by WProp
double WPropfact = WProp(plbar, pl);
return (2*WPropfact*(Mmmm+Mpmm)/24./4.)/(pa-p1-pl-plbar).m2()/(p2+p3-pb).m2();
}
// W+Jets qqxCentral
double ME_WCenqqx_qq(HLV pa, HLV pb,HLV pl, HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove
){
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;}
HLV pq, pqbar, p1, p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am, T4bm, T1ap, T4bp;
if(!(aqlinepa)){
T1ap = TCurrent(p1, true, pa, true);
T1am = TCurrent(p1, false, pa, false);}
else if(aqlinepa){
T1ap = TCurrent(pa, true, p1, true);
T1am = TCurrent(pa, false, p1, false);}
if(!(aqlinepb)){
T4bp = TCurrent(p4, true, pb, true);
T4bm = TCurrent(p4, false, pb, false);}
else if(aqlinepb){
T4bp = TCurrent(pb, true, p4, true);
T4bm = TCurrent(pb, false, p4, false);}
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc = MUncrossW(pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xcro = MCrossW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xsym = MSymW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is
// also the choice in qqbar helicity.
// (- - hel choice)
COM M_mmUnc = (((Xunc).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc).contract(T1am,1)).contract(T4bp,1));
COM M_mpCro = (((Xcro).contract(T1am,1)).contract(T4bp,1));
COM M_mpSym = (((Xsym).contract(T1am,1)).contract(T4bp,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc).contract(T1ap,1)).contract(T4bp,1));
COM M_ppCro = (((Xcro).contract(T1ap,1)).contract(T4bp,1));
COM M_ppSym = (((Xsym).contract(T1ap,1)).contract(T4bp,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
double WPropfact = WProp(plbar, pl);
amp*=WPropfact;
return amp;
}
// no wqq emission
double ME_W_Cenqqx_qq(HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
int nbelow, bool forwards
){
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;
}
if (!forwards){ //If Emission from Leg a instead, flip process.
HLV dummymom = pa;
bool dummybool= aqlinepa;
int dummyint = nabove;
pa = pb;
pb = dummymom;
std::reverse(partons.begin(),partons.end());
qqxmarker = !(qqxmarker);
aqlinepa = aqlinepb;
aqlinepb = dummybool;
nabove = nbelow;
nbelow = dummyint;
}
HLV pq, pqbar, p1,p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am(0.), T1ap(0.);
if(!(aqlinepa)){
T1ap = TCurrent(p1, true, pa, true);
T1am = TCurrent(p1, false, pa, false);}
else if(aqlinepa){
T1ap = TCurrent(pa, true, p1, true);
T1am = TCurrent(pa, false, p1, false);}
Tensor <1> T4bm = jW(pb, p4, plbar, pl, aqlinepb);
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc_m = MUncross(pa, pq, pqbar,partons, false, nabove);
Tensor<2> Xcro_m = MCross( pa, pq, pqbar,partons, false, nabove);
Tensor<2> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, false, nabove);
Tensor<2> Xunc_p = MUncross(pa, pq, pqbar,partons, true, nabove);
Tensor<2> Xcro_p = MCross( pa, pq, pqbar,partons, true, nabove);
Tensor<2> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, true, nabove);
// (- - hel choice)
COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
double WPropfact = WProp(plbar, pl);
amp*=WPropfact;
return amp;
}
diff --git a/t/ME_data/ME_Wm.dat b/t/ME_data/ME_Wm.dat
index 95149c7..cb0328e 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.510830208e-07
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
-1.381537057737695e-08
+3.025959937e-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
+8.291732575e-07
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
-3.538161996944155e-08
+1.546221661e-07
0.0004236100189
3.617420289e-06
4.516879429e-05
3.955718614e-05
0.0001365328248
5.126170478e-05
-2.005036424e-07
+1.077742001e-06
1.018758481e-08
-6.630491325177325e-09
+2.310396111e-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
+2.377796617e-09
0.0001490868071
9.612832588e-11
1.654742249e-05
1.57073488e-07
1.609407803e-06
-2.782539582e-09
+6.946367688e-08
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
+4.136266286e-07
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
+5.57330681e-10
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
-1.905147471814478e-09
+1.905147472e-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.312972541815244e-07
+1.312972542e-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
+3.9300806e-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
-3.247374364864003e-05
+1.393919618e-07
+7.338405876e-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
-4.513834365350505e-11
+6.701394914e-12
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
-8.773550961839615e-07
+2.009104635e-08
+8.773550962e-07
9.824802351e-10
-1.631528338714485e-07
+1.631528339e-07
2.288354101e-05
2.035248598e-07
0.0001019103655
3.912289372e-06
5.329003459e-08
3.815264848e-10
1.388142926e-05
-2.066904312077937e-08
+2.066904312e-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
+1.345426643e-06
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
-2.619064730671114e-08
+2.619064731e-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.003162717e-06
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
-1.98544141222215e-07
+1.985441412e-07
6.929088895e-08
-2.565897108e-07
+2.512260458e-06
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
+2.314556691e-08
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
+9.443687582e-10
+1.363026913e-06
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
-1.509519246185828e-07
+1.509519246e-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
+7.778718617e-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
+1.924149804e-05
7.698409214e-06
-1.111533673613149e-06
+1.221205448e-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
-6.869368158598052e-08
+6.869368159e-08
4.672444205e-09
7.02639721e-05
7.305366536e-07
2.162269825e-05
0.0004328624246
1.653264992e-07
1.553411281e-05
5.165593047e-10
-2.153731042127394e-10
-1.044312456634803e-08
+2.153731042e-10
+8.592408471e-09
8.974979512e-09
7.520620911e-07
-6.416418106e-08
+2.080333069e-07
1.653991504e-07
9.9043617e-09
6.542874961e-09
0.0005318553158
1.144555776e-05
-8.760870822e-10
+6.936210932e-05
7.62010085e-08
-3.426636308e-07
+1.834443898e-08
8.099206646e-06
-1.154675564e-09
+1.946663573e-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
+2.519238034e-06
3.662564708e-07
1.646271735e-07
5.373793619e-07
4.353476735e-08
8.727968448e-08
1.053026857e-05
-8.908419895e-07
+1.083603462e-06
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
-4.652418777294805e-10
-1.639422266502427e-06
+4.652418777e-10
+1.639422267e-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
+4.69742523e-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
+1.224885671e-06
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
+5.066538235e-11
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
+8.494930984e-08
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
+2.708628413e-06
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
-3.356011806202945e-07
-4.178663868e-10
+3.356011806e-07
+5.603569089e-09
9.856721257e-08
2.793861685e-07
5.195375958e-06
2.632243695e-06
6.350302361e-07
-3.302377879428893e-10
+3.302377879e-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
+2.297705942e-12
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
-9.719373728155117e-11
-1.732767627245824e-11
+9.719373728e-11
+1.732767627e-11
4.41946888e-11
5.657095809e-11
5.288646625e-07
5.324238867e-07
2.221589565e-07
2.48971658e-09
-1.285025395e-12
+2.994471157e-09
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
+1.251189973e-08
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
-4.995240386235438e-08
+4.995240386e-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
-1.542558411660617e-10
+1.542558412e-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
+7.129850446e-08
1.481220144e-09
9.13398584e-07
9.416502413e-11
1.481517933e-08
2.317808347e-06
-2.429491433e-08
+3.901257454e-08
4.627625989e-12
2.494572686e-08
9.094749026e-10
-1.806091619e-08
+1.731947437e-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
+1.037626188e-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
-3.444930056306798e-09
+3.444930056e-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.0756188221237e-09
-8.186541654e-13
+1.075618822e-09
+1.712525448e-09
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
-7.736937051663231e-10
+7.736937052e-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
+3.68209028e-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
-1.712560508754375e-13
+5.890823115e-10
+1.712560509e-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
+8.959130596e-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
+4.419441096e-10
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
+1.388892596e-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.364224996e-09
6.365307921e-10
-1.704476599e-14
+7.724490255e-10
2.368415063e-08
1.761682817e-07
-8.032268925e-11
+3.074395438e-08
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
+4.244764945e-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
+4.11857413e-09
9.269115967e-07
-1.080660039831632e-10
+1.08066004e-10
1.045618377e-06
4.94758944e-07
-3.094362842e-13
+6.160229109e-12
1.351633908e-08
2.784256537e-09
1.934517844e-08
-5.314666752529101e-13
+5.314666753e-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.019257173887625e-10
+3.679484278e-09
8.938611559e-10
5.601588127e-09
1.526952592e-05
3.497239887e-10
9.33354413e-08
1.359301132e-08
6.404370059e-11
1.535628496e-11
-2.2434943877662e-09
+2.243494388e-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 53c7f25..79ae7bf 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
+1.064845939e-10
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
+5.96828212e-08
1.972891892e-07
7.870126436e-08
1.027083098e-09
1.68746816e-07
6.312089294e-10
-3.505278148e-07
+3.391835464e-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
+6.037059083e-08
1.535388646e-05
8.248395671e-07
3.681439971e-07
6.978801279e-07
-2.186342703e-08
+0.0001134623229
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
-5.524153757135919e-10
+5.524153757e-10
5.4103721e-09
9.793444055e-05
1.21018852e-06
5.403237195e-06
9.396037827e-07
-9.218690916306924e-07
+9.218690916e-07
1.916819566e-05
2.177871221e-06
0.0003020897612
6.197213888e-08
-2.556142456e-08
+7.288410467e-08
1.364194249e-05
1.038835664e-05
1.181892855e-07
5.005165507e-06
-3.611225631803286e-08
+3.611225632e-08
1.221974589e-08
1.942760585e-05
3.734309313e-06
0.0001093995311
2.068051727e-06
8.857501175e-06
4.173118322e-06
-2.109563994e-08
+6.957949482e-09
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
+8.361499159e-09
6.131289047e-07
6.896616304e-06
2.90239901e-05
-5.756729461e-09
+7.595782532e-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
-6.867301068547712e-10
+8.760574342e-08
+6.867301069e-10
1.796235147e-07
3.651695186e-07
1.227417157e-07
-6.64659657953662e-06
+6.64659658e-06
8.406970187e-05
0.001885290476
0.0001570330885
2.314708659e-06
0.0002982727873
0.003229609179
2.42756184e-10
1.896039184e-06
-1.190944745856344e-08
+1.190944746e-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
-1.938501144056279e-06
+1.938501144e-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
+2.978511607e-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
+3.477529199e-06
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
+2.346067359e-05
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
-4.364063474995147e-08
+4.364063475e-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
+1.89348364e-07
+2.619177604e-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
-1.463957005451528e-06
+1.072050501e-05
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
+7.295054917e-08
5.765824783e-07
-9.011588796e-11
+3.265782722e-08
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
+1.608251404e-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
-6.013623386455616e-09
+6.013623386e-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
+4.562416832e-07
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
+3.652137087e-08
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
+2.06149244e-07
1.315444637e-09
3.270418535e-05
2.302140766e-07
4.002096106e-06
-9.910699640329833e-07
+1.609778389e-06
1.156198845e-08
4.49026944e-08
0.0002105189957
1.502882545e-10
0.0003938508954
2.359522705e-06
-5.766519563931604e-05
+5.766519564e-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.172768568716882e-10
-4.982454226859217e-07
+1.172768569e-10
+4.982454227e-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
-1.924092071480528e-06
+1.924092071e-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.3533254e-08
3.330632607e-05
1.95151118e-06
1.384784661e-05
-2.148337292837028e-07
+2.148337293e-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
+7.15420681e-08
0.0002157285841
9.440044669e-07
6.51295213e-08
1.03903577e-08
0.0001795918463
1.962893454e-05
-7.871495075e-10
+1.739516144e-07
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
-3.785295788394707e-09
+3.785295788e-09
2.390481542e-05
5.308836572e-08
4.670595432e-05
6.758171771e-05
-1.01803578427425e-09
+3.147635808e-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
-2.591411922436858e-07
+2.591411922e-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.683391447e-05
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
-2.466303925798377e-09
+2.466303926e-09
1.587801136e-05
5.253580132e-08
1.186388137e-06
2.513054181e-08
2.878235795e-05
1.268594957e-09
-4.036344625890604e-07
+4.036344626e-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.79685454e-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
+3.501999566e-09
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
+8.065384687e-10
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
+3.504032934e-08
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
+1.641630928e-10
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
+1.087174385e-08
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
+9.666274841e-10
7.465321603e-10
1.354621048e-10
2.460861072e-08
1.087512798e-09
2.736687033e-07
5.113035943e-10
6.306214627e-07
-4.283779224e-11
+3.166204088e-09
2.044825034e-08
6.895612973e-10
-6.338164558e-12
+6.672832434e-11
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.781922807e-08
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
+7.883829674e-09
3.017038977e-08
6.10674516e-07
7.43296019e-09
3.802076048e-08
1.356542979e-08
7.945938899e-11
-3.638927521e-12
+6.146929649e-09
1.061272859e-06
3.425627163e-09
4.724471432e-09
2.264988832e-06
5.206993128e-12
1.903839604e-11
5.585175902e-09
-9.146019191e-13
+5.143392872e-12
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
+2.202228399e-10
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
+1.266172515e-09
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
+6.965999199e-10
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
-1.929150342315743e-10
+1.929150342e-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.146479379e-09
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
-5.552870709834466e-08
+5.55287071e-08
diff --git a/t/test_ME_generic.cc b/t/test_ME_generic.cc
index 290fc1c..e2ad6d9 100644
--- a/t/test_ME_generic.cc
+++ b/t/test_ME_generic.cc
@@ -1,132 +1,199 @@
/**
* \brief Generic tester for the ME for a given set of PSP
*
* \note reference weights and PSP (as LHE file) have to be given as
* _individual_ files
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <algorithm>
#include <fstream>
#include <math.h>
#include <random>
#include "HEJ/Event.hh"
#include "HEJ/EventReader.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/stream.hh"
#include "HEJ/YAMLreader.hh"
constexpr double alpha_s = 0.118;
constexpr double ep = 1e-5;
+constexpr double ep_mirror = 1e-3;
void shuffle_particles(HEJ::Event::EventData & ev) {
static std::mt19937_64 ran{0};
+ // incoming
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
- std::shuffle(begin(ev.outgoing), end(ev.outgoing), ran);
+ // outgoing (through index)
+ auto old_outgoing = std::move(ev.outgoing);
+ std::vector<size_t> idx(old_outgoing.size());
+ std::iota(idx.begin(), idx.end(), 0);
+ std::shuffle(begin(idx), end(idx), ran);
+ ev.outgoing.clear();
+ ev.outgoing.reserve(old_outgoing.size());
+ for(size_t i: idx) {
+ ev.outgoing.emplace_back(std::move(old_outgoing[i]));
+ }
+ // find decays again
+ if(!ev.decays.empty()){
+ auto old_decays = std::move(ev.decays);
+ ev.decays.clear();
+ for(size_t i=0; i<idx.size(); ++i) {
+ auto decay = old_decays.find(idx[i]);
+ if(decay != old_decays.end())
+ ev.decays.emplace(i, std::move(decay->second));
+ }
+ }
}
enum MEComponent {tree, virt};
MEComponent guess_component(std::string const & data_file) {
if(data_file.find("virt") != data_file.npos) return MEComponent::virt;
return MEComponent::tree;
}
+HEJ::Event::EventData mirror_event(HEJ::Event::EventData ev){
+ for(auto & part: ev.incoming){
+ auto & p{ part.p };
+ p.reset(p.px(),p.py(),-p.pz(),p.E());
+ }
+ for(auto & part: ev.outgoing){
+ auto & p{ part.p };
+ p.reset(p.px(),p.py(),-p.pz(),p.E());
+ }
+ for(auto & decay: ev.decays){
+ for(auto & part: decay.second){
+ auto & p{ part.p };
+ p.reset(p.px(),p.py(),-p.pz(),p.E());
+ }
+ }
+ return ev;
+}
+
int main(int argn, char** argv){
if(argn != 4 && argn != 5){
std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml ME_weights input_file.lhe\n\n";
return EXIT_FAILURE;
}
bool OUTPUT_MODE = false;
if(argn == 5 && std::string("OUTPUT")==std::string(argv[4]))
OUTPUT_MODE = true;
const HEJ::Config config = HEJ::load_config(argv[1]);
std::fstream wgt_file;
if ( OUTPUT_MODE ) {
std::cout << "_______________________USING OUTPUT MODE!_______________________" << std::endl;
wgt_file.open(argv[2], std::fstream::out);
wgt_file.precision(10);
} else {
wgt_file.open(argv[2], std::fstream::in);
}
auto reader{ HEJ::make_reader(argv[3])};
const auto component = guess_component(argv[2]);
HEJ::MatrixElement ME{
[](double){ return alpha_s; },
HEJ::to_MatrixElementConfig(config)
};
double max_ratio = 0.;
size_t idx_max_ratio = 0;
HEJ::Event ev_max_ratio(HEJ::Event::EventData{}.cluster(
config.resummation_jets.def,0
)
);
double av_ratio = 0;
size_t i = 0;
while(reader->read_event()){
++i;
HEJ::Event::EventData data{reader->hepeup()};
shuffle_particles(data);
+ HEJ::Event::EventData data_mirror{mirror_event(data)};
+ shuffle_particles(data_mirror);
+
HEJ::Event event{
data.cluster(
config.resummation_jets.def,
config.resummation_jets.min_pt
)
};
+ HEJ::Event event_mirror{
+ data_mirror.cluster(
+ config.resummation_jets.def,
+ config.resummation_jets.min_pt
+ )
+ };
const double our_ME = (component == MEComponent::tree)?
ME.tree(event).central:
ME.virtual_corrections(event).central
;
if(!isfinite(our_ME)){
std::cerr << "Found non-finite ME ("<< our_ME <<")\n" << event << std::endl;
return EXIT_FAILURE;
}
+ const double ME_mirror = (component == MEComponent::tree)?
+ ME.tree(event_mirror).central:
+ ME.virtual_corrections(event_mirror).central
+ ;
+ if(!isfinite(ME_mirror)){
+ std::cerr << "Found non-finite ME ("<< ME_mirror <<")\n" << event_mirror << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ if(std::abs(our_ME/ME_mirror-1.)>ep_mirror){
+ size_t precision(std::cout.precision());
+ std::cerr.precision(16);
+ std::cerr<< "z-Mirrored ME gives different result " << i << "\n"
+ <<our_ME << " vs " << ME_mirror << " => difference: "
+ << std::abs(our_ME/ME_mirror-1.) << "\n" << event
+ << "\nmirrored:\n" << event_mirror << std::endl;
+ std::cerr.precision(precision);
+ return EXIT_FAILURE;
+ }
if ( OUTPUT_MODE ) {
wgt_file << our_ME << std::endl;
} else {
std::string line;
if(!std::getline(wgt_file,line)) break;
const double ref_ME = std::stod(line);
const double diff = std::abs(our_ME/ref_ME-1.);
av_ratio+=diff;
if( diff > max_ratio ) {
max_ratio = diff;
idx_max_ratio = i;
ev_max_ratio = event;
}
if( diff > ep ){
size_t precision(std::cout.precision());
std::cerr.precision(16);
std::cerr<< "Large difference in PSP " << i << "\nis: "<<our_ME
- << " should: " << ref_ME << " => difference: " << diff
+ << " should: " << ref_ME << " => difference: " << diff << "\n"
<< event << std::endl;
std::cerr.precision(precision);
return EXIT_FAILURE;
}
}
}
wgt_file.close();
if ( i<100 )
throw std::invalid_argument{"Not enough PSP tested"};
if ( !OUTPUT_MODE ) {
size_t precision(std::cout.precision());
std::cout.precision(16);
std::cout << "Avg ratio after " << i << " PSP: " << av_ratio/i << std::endl;
std::cout << "maximal ratio at " << idx_max_ratio << ": " << max_ratio << std::endl;
std::cout.precision(precision);
}
return EXIT_SUCCESS;
}
diff --git a/t/test_classify.cc b/t/test_classify.cc
index 2351e3c..91a3245 100644
--- a/t/test_classify.cc
+++ b/t/test_classify.cc
@@ -1,835 +1,855 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <iostream>
#include <random>
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#define ASSERT(x) if(!(x)) { \
throw std::logic_error("Assertion '" #x "' failed."); \
}
namespace {
const fastjet::JetDefinition jet_def{fastjet::JetAlgorithm::antikt_algorithm, 0.4};
const double min_jet_pt{30.};
const std::array<std::string, 6> all_quarks{"-4","-1","1","2","3","4"};
const std::array<std::string, 7> all_partons{"g","-2","-1","1","2","3","4"};
const std::array<std::string, 3> all_bosons{"h", "Wp", "Wm"};
static std::mt19937_64 ran{0};
void shuffle_particles(HEJ::Event::EventData & ev) {
+ // incoming
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
- std::shuffle(begin(ev.outgoing), end(ev.outgoing), ran);
+ // outgoing (through index)
+ auto old_outgoing = std::move(ev.outgoing);
+ std::vector<size_t> idx(old_outgoing.size());
+ std::iota(idx.begin(), idx.end(), 0);
+ std::shuffle(begin(idx), end(idx), ran);
+ ev.outgoing.clear();
+ ev.outgoing.reserve(old_outgoing.size());
+ for(size_t i: idx) {
+ ev.outgoing.emplace_back(std::move(old_outgoing[i]));
+ }
+ // find decays again
+ if(!ev.decays.empty()){
+ auto old_decays = std::move(ev.decays);
+ ev.decays.clear();
+ for(size_t i=0; i<idx.size(); ++i) {
+ auto decay = old_decays.find(idx[i]);
+ if(decay != old_decays.end())
+ ev.decays.emplace(i, std::move(decay->second));
+ }
+ }
}
// if pos_boson = -1 (or not implemented) -> no boson
// njet==7 is special: has less jets, i.e. multiple parton in one jet,
// pos_boson < 0 to select process (see list for details)
HEJ::Event::EventData get_process(int const njet, int const pos_boson){
HEJ::Event::EventData ev;
if(njet == 2){
switch(pos_boson){
case 0:
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 198, 33, -170, 291}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-154, 68, 44, 174}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -44, -101, 88, 141}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -322, 322}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 284, 284}, {}};
return ev;
case 1:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -6, 82, -159, 179}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 195, -106, 74, 265}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-189, 24, 108, 219}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -320, 320}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 343, 343}, {}};
return ev;
case 2:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -80, -80, -140, 180}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -60, -32, 0, 68}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 140, 112, 177, 281}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -246, 246}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 283, 283}, {}};
return ev;
default:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -72, 24, 18, 78}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 72, -24, 74, 106}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -46, 46}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 138, 138}, {}};
return ev;
}
}
if(njet == 3){
switch(pos_boson){
case 0:
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 152, -117, -88, 245}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-146, 62, -11, 159}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 126, -72, 96, 174}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-132, 127, 144, 233}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -335, 335}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 476, 476}, {}};
return ev;
case 1:
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-191, 188, -128, 297}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 199, 72, -76, 257}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 184, -172, -8, 252}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-192, -88, 54, 218}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -591, 591}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 433, 433}, {}};
return ev;
case 2:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -42, 18, -49, 67}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, -54, -28, 62}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 99, 32, -16, 163}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -45, 4, 72, 85}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -199, 199}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 178, 178}, {}};
return ev;
case 3:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -65, -32, -76, 105}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -22, 31, -34, 51}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, -67, -36, 77}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 99, 68, -4, 173}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -278, 278}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 128, 128}, {}};
return ev;
default:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -90, -135, 30, 165}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-108, 198, 76, 238}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 198, -63, 126, 243}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -207, 207}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 439, 439}, {}};
return ev;
}
}
if(njet == 4){
switch(pos_boson){
case 0:
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 199, 72, -76, 257}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-200, -155, -64, 261}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 198, 194, 57, 283}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 1, 32, 8, 33}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-198, -143, 186, 307}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -515, 515}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 626, 626}, {}};
return ev;
case 1:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 198, 61, -162, 263}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 199, 72, -76, 257}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-200, 135, 144, 281}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-198, -186, 171, 321}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 1, -82, 122, 147}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -535, 535}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 734, 734}, {}};
return ev;
case 2:
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-180, -27, -164, 245}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-108, 78, -36, 138}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 196, -189, 68, 307}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-107, 136, 76, 189}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 199, 2, 178, 267}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -512, 512}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 634, 634}, {}};
return ev;
case 3:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, -30, -84, 90}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -72, 22, -96, 122}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 0, -51, 85}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 64, 72, -81, 177}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -64, 84, 116}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -409, 409}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 181, 181}, {}};
return ev;
case 4:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -72, -49, -72, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, 0, -36, 60}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 54, -36, 66}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, -77, -56, 117}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 64, 72, -81, 177}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -407, 407}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 126, 126}, {}};
return ev;
default:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 248, -56, -122, 282}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 249, 30, -10, 251}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-249, -18, 26, 251}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-248, 44, 199, 321}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -506, 506}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 599, 599}, {}};
return ev;
}
}
if(njet == 6){
switch(pos_boson){
case 0:
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 349, 330, -94, 505}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-315, -300, 0, 435}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 347, 306, 18, 463}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-249, -342, 162, 453}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 345, 312, 284, 545}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-324, -126, 292, 454}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-153, -180, 304, 385}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -1137, 1137}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 2103, 2103}, {}};
return ev;
case 1:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 242, 241, -182, 387}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 243, 238, -190, 409}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-218, -215, -74, 315}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-224, -224, 112, 336}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 241, 182, 154, 339}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -53, -234, 126, 271}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-231, 12, 156, 279}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -1117, 1117}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 1219, 1219}, {}};
return ev;
case 2:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 151, 102, -42, 187}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -86, -46, -17, 99}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 152, 153, 0, 249}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -60, -135, 64, 161}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 150, 123, 110, 223}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-154, -49, 98, 189}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-153, -148, 144, 257}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -504, 504}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 861, 861}, {}};
return ev;
case 3:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 198, 197, -66, 287}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-198, -189, -54, 279}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-200, -64, 2, 210}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 199, 158, 6, 283}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-199, -184, 172, 321}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 196, 168, 177, 313}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 4, -86, 92, 126}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -745, 745}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 1074, 1074}, {}};
return ev;
case 4:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 151, 102, -42, 187}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -86, -133, -14, 159}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-154, -104, -8, 186}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -60, 11, 0, 61}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 152, 153, 0, 249}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 150, 125, 90, 215}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-153, -154, 126, 251}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -578, 578}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 730, 730}, {}};
return ev;
case 5:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -90, -94, 131}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 82, -74, 111}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -80, -64, 105}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -25, -36, 65}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -16, 101}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 68, 92, -18, 170}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -78, 54, 95}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -513, 513}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 265, 265}, {}};
return ev;
case 6:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 198, 197, -66, 287}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 4, -84, -18, 86}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-198, -60, -36, 210}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 196, -78, -36, 214}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-200, 45, 0, 205}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-199, -178, 2, 267}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 199, 158, 6, 283}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -850, 850}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 702, 702}, {}};
return ev;
default:
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-350, -112, -280, 462}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 347, 266, -322, 543}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-349, -314, -38, 471}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 349, 348, 12, 493}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-342, -54, 23, 347}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 345, -134, 138, 395}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -1589, 1589}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 1122, 1122}, {}};
return ev;
}
}
if(njet == 7){
switch(pos_boson){
case -1: // jet idx: -1 0 1 2 3 4 5
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, -18, -42, 47}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, 26, -18, 35}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, -24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -54, -6, 59}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -44, 8, 45}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -96, 44, 116}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, 88, 133}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -249, 249}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 299, 299}, {}};
return ev;
case -2: // jet idx: 0 1 2 3 4 -1 -1
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -80, -64, 105}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -84, -12, 85}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, 88, 133}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -24, 62, 82}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, -18, 42, 47}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, 20, 60, 65}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -215, 215}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 415, 415}, {}};
return ev;
case -3: // jet idx: 0 0 1 2 3 4 5
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -86, -70, 111}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -52, -36, 65}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -44, 109}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -60, 5, 77}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 92, 8, 93}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -80, 64, 105}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -361, 361}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 312, 312}, {}};
return ev;
case -4: // jet idx: 0 1 2 3 4 5 5
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -40, -56, 69}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -88, 133}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -84, -72, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 92, 8, 93}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -58, 30, 67}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -96, 96, 144}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -395, 395}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 337, 337}, {}};
return ev;
case -5: // jet idx: 0 1 -1 -1 2 3 4
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -64, -80, 105}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -16, 101}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, 20, 0, 25}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, -10, 2, 15}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -72, 54, 102}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -60, 48, 77}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -253, 253}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 285, 285}, {}};
return ev;
case -6: // jet idx: 0 0 0 1 2 2 3
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -60, -48, 77}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -52, -36, 65}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -96, -58, 122}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, -24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -16, 101}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 92, 8, 93}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -70, 14, 75}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -403, 403}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 243, 243}, {}};
return ev;
case -7: // jet idx: 0 1 1 2 2 3 4
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -46, -46, 69}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -90, -70, 115}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -78, -54, 95}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 88, -28, 93}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -16, 101}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -60, 5, 77}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -424, 424}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 239, 239}, {}};
return ev;
case -8: // jet idx: 0 1 2 2 2 3 4
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -84, -84, 121}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 92, -8, 93}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -36, 0, 39}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -62, 10, 63}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -96, 19, 109}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, 88, 133}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -311, 311}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 360, 360}, {}};
return ev;
}
}
throw HEJ::unknown_option{"unkown process"};
}
bool couple_quark(std::string const & boson, std::string & quark){
if(abs(HEJ::to_ParticleID(boson)) == HEJ::ParticleID::Wp){
auto qflav{ HEJ::to_ParticleID(quark) };
if(!HEJ::is_anyquark(qflav)) return false;
const int W_charge = HEJ::to_ParticleID(boson)>0?1:-1;
if(W_charge*qflav < 0 && !(abs(qflav)%2)) return false; // not anti-down
if(W_charge*qflav > 0 && (abs(qflav)%2)) return false; // not up
quark=std::to_string(qflav-W_charge);
}
return true;
}
// create event corresponding from given Configuration
// overwrite_boson to force a specific boson position, indepentent from input
// (useful for njet == 7)
HEJ::Event parse_configuration(
std::array<std::string,2> const & in, std::vector<std::string> const & out,
int const overwrite_boson = 0
){
auto boson = std::find_if(out.cbegin(), out.cend(),
[](std::string id){ return !HEJ::is_parton(HEJ::to_ParticleID(id)); });
int const pos_boson = (overwrite_boson!=0)?overwrite_boson:
((boson == out.cend())?-1:std::distance(out.cbegin(), boson));
size_t njets = out.size();
if(boson != out.cend()) --njets;
HEJ::Event::EventData ev{get_process(njets, pos_boson)};
ASSERT((pos_boson<0) || (ev.outgoing[pos_boson].type == HEJ::ParticleID::higgs));
for(size_t i=0; i<out.size(); ++i){
ev.outgoing[i].type = HEJ::to_ParticleID(out[i]);
}
for(size_t i=0; i<in.size(); ++i){
ev.incoming[i].type = HEJ::to_ParticleID(in[i]);
}
shuffle_particles(ev);
return ev.cluster(jet_def, min_jet_pt);
}
bool match_expectation( HEJ::event_type::EventType expected,
std::array<std::string,2> const & in, std::vector<std::string> const & out,
int const overwrite_boson = 0
){
HEJ::Event ev{parse_configuration(in,out,overwrite_boson)};
if(ev.type() != expected){
std::cerr << "Expected type " << HEJ::event_type::name(expected)
<< " but found " << HEJ::event_type::name(ev.type()) << "\n" << ev;
auto jet_idx{ ev.particle_jet_indices(ev.jets()) };
std::cout << "Particle Jet indices: ";
for(int const i: jet_idx)
std::cout << i << " ";
std::cout << std::endl;
return false;
}
return true;
}
bool check_fkl(){
using namespace HEJ;
for(std::string const & first: all_partons) // all quark flavours
for(std::string const & last: all_partons){
for(int njet=2; njet<=6; ++njet){ // all multiplicities
if(njet==5) continue;
std::array<std::string,2> base_in{first,last};
std::vector<std::string> base_out(njet, "g");
base_out.front() = first;
base_out.back() = last;
if(!match_expectation(event_type::FKL, base_in, base_out))
return false;
for(auto const & boson: all_bosons) // any boson
for(int pos=0; pos<=njet; ++pos){ // at any position
auto in{base_in};
auto out{base_out};
// change quark flavours for W
int couple_idx{ std::uniform_int_distribution<int>{0,1}(ran) };
if(!couple_quark(boson, couple_idx?out.back():out.front()))
continue;
out.insert(out.begin()+pos, boson);
if(!match_expectation(event_type::FKL, in, out))
return false;
}
}
for(int i=-1; i>-9;--i){ // jet setup doesn't matter for FKL
std::array<std::string,2> base_in{first,last};
std::vector<std::string> base_out(7, "g");
base_out.front() = first;
base_out.back() = last;
if(!match_expectation(event_type::FKL, base_in, base_out, i))
return false;
}
}
return true;
}
bool check_uno(){
using namespace HEJ;
auto const b{ event_type::unob };
auto const f{ event_type::unof };
for(std::string const & uno: all_quarks) // all quark flavours
for(std::string const & fkl: all_partons){
for(int njet=3; njet<=6; ++njet){ // all multiplicities >2
if(njet==5) continue;
for(int i=0; i<2; ++i){ // forward & backwards
std::array<std::string,2> base_in;
std::vector<std::string> base_out(njet, "g");
const int uno_pos = i?1:(njet-2);
const int fkl_pos = i?(njet-1):0;
base_in[i?0:1] = uno;
base_in[i?1:0] = fkl;
base_out[uno_pos] = uno;
base_out[fkl_pos] = fkl;
auto expectation{ i?b:f };
if( !match_expectation(expectation, base_in, base_out) )
return false;
for(auto const & boson: all_bosons){ // any boson
// at any position (higgs only inside FKL chain)
int start = 0;
int end = njet;
if(to_ParticleID(boson) == pid::higgs){
start = i?(uno_pos+1):fkl_pos;
end = i?(fkl_pos+1):uno_pos;
}
for(int pos=start; pos<=end; ++pos){
auto in{base_in};
auto out{base_out};
// change quark flavours for W
bool couple_idx{ std::uniform_int_distribution<int>{0,1}(ran) };
if(!couple_quark(boson, couple_idx?out[fkl_pos]:out[uno_pos]))
continue;
out.insert(out.begin()+pos, boson);
if(!match_expectation(expectation, in, out))
return false;
}
}
}
}
// test allowed jet configurations
if(!(match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -1)
&& match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -2)
&& match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -3)
&& match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -4)
&& match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -5)
&& match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -5)
&& match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -6)
&& match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -7)
&& match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -7)
&& match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -8)
&& match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -8)
))
return false;
}
return true;
}
bool check_extremal_qqx(){
using namespace HEJ;
auto const b{ event_type::qqxexb };
auto const f{ event_type::qqxexf };
for(std::string const & qqx: all_quarks) // all quark flavours
for(std::string const & fkl: all_partons){
std::string const qqx2{ std::to_string(HEJ::to_ParticleID(qqx)*-1) };
for(int njet=3; njet<=6; ++njet){ // all multiplicities >2
if(njet==5) continue;
for(int i=0; i<2; ++i){ // forward & backwards
std::array<std::string,2> base_in;
std::vector<std::string> base_out(njet, "g");
const int qqx_pos = i?0:(njet-2);
const int fkl_pos = i?(njet-1):0;
base_in[i?0:1] = "g";
base_in[i?1:0] = fkl;
base_out[fkl_pos] = fkl;
base_out[qqx_pos] = qqx;
base_out[qqx_pos+1] = qqx2;
auto expectation{ i?b:f };
if( !match_expectation(expectation, base_in, base_out) )
return false;
for(auto const & boson: all_bosons){ // any boson
// at any position (higgs only inside FKL chain)
int start = 0;
int end = njet;
if(to_ParticleID(boson) == pid::higgs){
start = i?(qqx_pos+2):fkl_pos;
end = i?(fkl_pos+1):qqx_pos;
}
for(int pos=start; pos<=end; ++pos){
auto in{base_in};
auto out{base_out};
// change quark flavours for W
bool couple_idx{ std::uniform_int_distribution<int>{0,1}(ran) };
if(couple_idx || !couple_quark(boson, out[fkl_pos]) ){
// (randomly) try couple to FKL, else fall-back to qqx
if(!couple_quark(boson, out[qqx_pos]))
couple_quark(boson, out[qqx_pos+1]);
}
out.insert(out.begin()+pos, boson);
if(!match_expectation(expectation, in, out))
return false;
}
}
}
}
// test allowed jet configurations
if(!(match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -1)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -2)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -3)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -4)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -5)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -5)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -6)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -7)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -7)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -8)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -8)
))
return false;
}
return true;
}
bool check_central_qqx(){
using namespace HEJ;
auto const t{ event_type::qqxmid };
for(std::string const & qqx: all_quarks) // all quark flavours
for(std::string const & fkl1: all_partons)
for(std::string const & fkl2: all_partons){
std::string const qqx2{ std::to_string(HEJ::to_ParticleID(qqx)*-1) };
for(int njet=4; njet<=6; ++njet){ // all multiplicities >3
if(njet==5) continue;
for(int qqx_pos=1; qqx_pos<njet-2; ++qqx_pos){ // any qqx position
std::array<std::string,2> base_in;
std::vector<std::string> base_out(njet, "g");
base_in[0] = fkl1;
base_in[1] = fkl2;
base_out.front() = fkl1;
base_out.back() = fkl2;
base_out[qqx_pos] = qqx;
base_out[qqx_pos+1] = qqx2;
if( !match_expectation(t, base_in, base_out) )
return false;
for(auto const & boson: all_bosons) // any boson
for(int pos=0; pos<=njet; ++pos){ // at any position
if( to_ParticleID(boson) == pid::higgs
&& (pos==qqx_pos || pos==qqx_pos+1) )
continue;
auto in{base_in};
auto out{base_out};
// change quark flavours for W
int couple_idx{ std::uniform_int_distribution<int>{0,2}(ran) };
// (randomly) try couple to FKL, else fall-back to qqx
if( couple_idx == 0 && couple_quark(boson, out.front()) ){}
else if( couple_idx == 1 && couple_quark(boson, out.back()) ){}
else {
if(!couple_quark(boson, out[qqx_pos]))
couple_quark(boson, out[qqx_pos+1]);
}
out.insert(out.begin()+pos, boson);
if(!match_expectation(t, in, out))
return false;
}
}
}
// test allowed jet configurations (non exhaustive collection)
if(!(match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -1)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -2)
&& match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -3)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -3)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -4)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -4)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -5)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -6)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -6)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -7)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -7)
&& match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -8)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -8)
))
return false;
}
return true;
}
// this checks a (non excessive) list of non-resummable states
bool check_non_resummable(){
auto type{ HEJ::event_type::non_resummable};
return
// 2j - crossing lines
match_expectation(type, {"g","2"}, {"2","g"})
&& match_expectation(type, {"-1","g"}, {"g","-1"})
&& match_expectation(type, {"1","-1"}, {"-1","1"})
&& match_expectation(type, {"g","2"}, {"2","g","h"})
&& match_expectation(type, {"1","2"}, {"2","h","1"})
&& match_expectation(type, {"1","-1"}, {"h","-1","1"})
&& match_expectation(type, {"g","2"}, {"Wp","1","g"})
&& match_expectation(type, {"1","-1"}, {"-2","Wp","1"})
&& match_expectation(type, {"4","g"}, {"g","3","Wp"})
&& match_expectation(type, {"1","-2"}, {"-1","Wm","1"})
&& match_expectation(type, {"g","3"}, {"4","g","Wm"})
&& match_expectation(type, {"1","3"}, {"Wm","4","1"})
// 2j - qqx
&& match_expectation(type, {"g","g"}, {"1","-1"})
&& match_expectation(type, {"g","g"}, {"-2","2","h"})
&& match_expectation(type, {"g","g"}, {"-4","Wp","3"})
&& match_expectation(type, {"g","g"}, {"Wm","-1","2"})
// 3j - crossing lines
&& match_expectation(type, {"g","4"}, {"4","g","g"})
&& match_expectation(type, {"-1","g"}, {"g","g","-1"})
&& match_expectation(type, {"1","3"}, {"3","g","1"})
&& match_expectation(type, {"-2","2"}, {"2","g","-2","h"})
&& match_expectation(type, {"-3","g"}, {"g","g","Wp","-4"})
&& match_expectation(type, {"1","-2"}, {"Wm","-1","g","1"})
&& match_expectation(type, {"-1","g"}, {"1","-1","-1"})
// higgs inside uno
&& match_expectation(type, {"-1","g"}, {"g","h","-1","g"})
&& match_expectation(type, {"-1","1"}, {"g","h","-1","1"})
&& match_expectation(type, {"g","2"}, {"g","2","h","g"})
&& match_expectation(type, {"-1","1"}, {"-1","1","h","g"})
// higgs outside uno
&& match_expectation(type, {"-1","g"}, {"h","g","-1","g"})
&& match_expectation(type, {"-1","1"}, {"-1","1","g","h"})
// higgs inside qqx
&& match_expectation(type, {"g","g"}, {"-1","h","1","g","g"})
&& match_expectation(type, {"g","g"}, {"g","-1","h","1","g"})
&& match_expectation(type, {"g","g"}, {"g","g","2","h","-2"})
// higgs outside qqx
&& match_expectation(type, {"g","g"}, {"h","-1","1","g","g"})
&& match_expectation(type, {"g","g"}, {"g","g","2","-2","h"})
// 4j - two uno
&& match_expectation(type, {"-2","2"}, {"g","-2","2","g"})
&& match_expectation(type, {"1","3"}, {"g","1","h","3","g"})
&& match_expectation(type, {"1","2"}, {"g","1","3","Wp","g"})
&& match_expectation(type, {"1","-2"}, {"g","Wm","1","-1","g"})
// 4j - two gluon outside
&& match_expectation(type, {"g","4"}, {"g","4","g","g"})
&& match_expectation(type, {"1","3"}, {"1","3","h","g","g"})
&& match_expectation(type, {"1","2"}, {"1","3","g","Wp","g"})
&& match_expectation(type, {"1","-2"}, {"1","Wm","-1","g","g"})
&& match_expectation(type, {"-1","g"}, {"g","g","-1","g"})
&& match_expectation(type, {"1","3"}, {"g","g","1","3","h"})
&& match_expectation(type, {"1","2"}, {"g","g","1","Wp","3"})
&& match_expectation(type, {"1","-2"}, {"Wm","g","g","1","-1"})
// 4j - ggx+uno
&& match_expectation(type, {"g","4"}, {"1","-1","4","g"})
&& match_expectation(type, {"2","g"}, {"g","2","-3","3"})
&& match_expectation(type, {"g","4"}, {"1","-1","h","4","g"})
&& match_expectation(type, {"2","g"}, {"g","2","-3","3","h"})
&& match_expectation(type, {"g","4"}, {"Wp","1","-1","3","g"})
&& match_expectation(type, {"2","g"}, {"g","2","-4","Wp","3"})
&& match_expectation(type, {"g","4"}, {"2","Wm","-1","4","g"})
&& match_expectation(type, {"2","g"}, {"g","2","Wp","-3","4"})
// 3j - crossing+uno
&& match_expectation(type, {"1","4"}, {"g","4","1"})
&& match_expectation(type, {"1","4"}, {"4","1","g"})
&& match_expectation(type, {"1","4"}, {"g","h","4","1"})
&& match_expectation(type, {"-1","-3"},{"Wm","g","-4","-1"})
&& match_expectation(type, {"1","4"}, {"3","1","Wp","g"})
&& match_expectation(type, {"1","4"}, {"3","1","g","h"})
// 3j - crossing+qqx
&& match_expectation(type, {"1","g"}, {"-1","1","g","1"})
&& match_expectation(type, {"1","g"}, {"-1","1","1","g"})
&& match_expectation(type, {"g","1"}, {"1","g","1","-1"})
&& match_expectation(type, {"g","1"}, {"g","1","1","-1"})
&& match_expectation(type, {"1","g"}, {"2","-2","g","1"})
&& match_expectation(type, {"1","g"}, {"2","-2","1","g"})
&& match_expectation(type, {"g","1"}, {"1","g","-2","2"})
&& match_expectation(type, {"g","1"}, {"g","1","-2","2"})
&& match_expectation(type, {"1","g"}, {"-1","1","h","g","1"})
&& match_expectation(type, {"1","g"}, {"-1","h","1","1","g"})
&& match_expectation(type, {"g","1"}, {"1","g","1","h","-1"})
&& match_expectation(type, {"g","1"}, {"h","g","1","1","-1"})
&& match_expectation(type, {"1","g"}, {"2","-2","1","g","h"})
&& match_expectation(type, {"g","1"}, {"g","h","1","-2","2"})
&& match_expectation(type, {"1","g"}, {"Wp","3","-4","g","1"})
&& match_expectation(type, {"3","g"}, {"-2","Wm","1","3","g"})
&& match_expectation(type, {"g","1"}, {"1","g","Wm","-3","4"})
&& match_expectation(type, {"g","-3"}, {"g","-3","-1","Wp","2"})
// 4j- gluon in qqx
&& match_expectation(type, {"g","1"}, {"1","g","-1","1"})
&& match_expectation(type, {"1","g"}, {"1","-1","g","1"})
&& match_expectation(type, {"g","1"}, {"1","g","Wm","-2","1"})
&& match_expectation(type, {"2","g"}, {"2","-2","g","Wp","1"})
&& match_expectation(type, {"g","g"}, {"Wp","3","g","-4","g"})
&& match_expectation(type, {"1","g"}, {"1","h","-1","g","1"})
// 6j - two qqx
&& match_expectation(type, {"g","g"}, {"1","-1","g","g","1","-1"})
&& match_expectation(type, {"g","g"}, {"1","-1","g","1","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","-1","g","1","-1"})
&& match_expectation(type, {"g","g"}, {"g","1","-1","1","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","1","-1","-1","g"})
&& match_expectation(type, {"g","g"}, {"h","1","-1","g","g","1","-1"})
&& match_expectation(type, {"g","g"}, {"1","Wp","-2","g","1","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","Wp","-1","g","1","-2"})
&& match_expectation(type, {"g","g"}, {"g","1","-1","Wm","2","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","2","-1","Wm","-1","g"})
// random stuff (can be non-physical)
&& match_expectation(type, {"g","g"}, {"1","-2","2","-1"}) // != 2 qqx
&& match_expectation(type, {"g","g"}, {"1","-2","2","g"}) // could be qqx
&& match_expectation(type, {"e+","e-"},{"1","-1"}) // bad initial state
&& match_expectation(type, {"1","e-"}, {"g","1","Wm"}) // bad initial state
&& match_expectation(type, {"h","g"}, {"g","g"}) // bad initial state
&& match_expectation(type, {"-1","g"}, {"-1","1","1"}) // bad qqx
&& match_expectation(type, {"-1","g"}, {"1","1","-1"}) // crossing in bad qqx
&& match_expectation(type, {"-1","g"}, {"-2","1","1","Wp"}) // bad qqx
&& match_expectation(type, {"1","2"}, {"1","-1","g","g","g","2"}) // bad qqx
&& match_expectation(type, {"1","2"}, {"1","-1","-2","g","g","2"}) // gluon in bad qqx
&& match_expectation(type, {"g","g"}, {"-1","2","g","g"}) // wrong back qqx
&& match_expectation(type, {"g","g"}, {"g","g","2","1"}) // wrong forward qqx
&& match_expectation(type, {"g","g"}, {"g","-2","1","g"}) // wrong central qqx
;
}
// Events failing the jet requirements, e.g. qqx inside one jet
//! @FIXME they are currently allowed
bool check_illegal_jets(){
auto type{ HEJ::event_type::non_resummable};
return true
// uno backward not in jet
&& match_expectation(type, {"1","g"}, {"g","1","g","g","g","g","g"}, -1)
// & also legal uno on other side
&& match_expectation(type, {"1","1"}, {"g","1","g","g","g","1","g"}, -1)
// qqx backward not in jet
&& match_expectation(type, {"g","2"}, {"-1","1","g","g","g","g","2"}, -1)
// uno forward not in jet
&& match_expectation(type, {"3","3"}, {"3","g","g","g","g","3","g"}, -2)
// qqx backward not in jet
&& match_expectation(type, {"g","g"}, {"g","g","g","g","g","-2","2"}, -2)
// uno backward in same jet
&& match_expectation(type, {"1","g"}, {"g","1","g","g","g","g","g"}, -3)
// & also legal uno on other side
&& match_expectation(type, {"1","1"}, {"g","1","g","g","g","1","g"}, -3)
// qqx backward in same jet
&& match_expectation(type, {"g","2"}, {"-4","4","g","g","g","g","2"}, -3)
// uno forward in same jet
&& match_expectation(type, {"3","2"}, {"3","g","g","g","g","2","g"}, -4)
// qqx backward in same jet
&& match_expectation(type, {"g","g"}, {"g","g","g","g","g","-2","2"}, -4)
// central qqx not in jet
&& match_expectation(type, {"1","2"}, {"1","g","-1","1","g","g","2"}, -5)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","-1","1","g","g","g","2"}, -6)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","g","g","g","2","-2","2"}, -6)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","-1","1","g","g","g","2"}, -7)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","g","g","3","-3","g","2"}, -7)
// central qqx in same jet
&& match_expectation(type, {"g","3"}, {"g","g","-2","2","g","g","3"}, -8)
// central qqx in same jet
&& match_expectation(type, {"g","-2"}, {"g","g","g","2","-2","g","-2"}, -8)
;
}
// Two boson states, that are currently not implemented
bool check_bad_FS(){
auto type{ HEJ::event_type::bad_final_state};
return
match_expectation(type, {"g","g"}, {"g","h","h","g"})
&& match_expectation(type, {"g","g"}, {"h","g","h","g"})
&& match_expectation(type, {"g","-1"}, {"g","h","Wp","-2"})
&& match_expectation(type, {"-3","-1"},{"-4","g","Wp","Wp","-2"})
&& match_expectation(type, {"-4","-1"},{"-3","Wp","g","Wm","-2"})
&& match_expectation(type, {"-4","-1"},{"g","-3","Wp","Wm","-2"})
;
}
}
int main() {
// tests for "no false negatives"
// i.e. all HEJ-configurations get classified correctly
if(!check_fkl()) return EXIT_FAILURE;
if(!check_uno()) return EXIT_FAILURE;
if(!check_extremal_qqx()) return EXIT_FAILURE;
if(!check_central_qqx()) return EXIT_FAILURE;
// test for "no false positive"
// i.e. non-resummable gives non-resummable
if(!check_non_resummable()) return EXIT_FAILURE;
if(!check_illegal_jets()) return EXIT_FAILURE;
if(!check_bad_FS()) return EXIT_FAILURE;
//! @TODO missing:
//! checks for partons sharing a jet
return EXIT_SUCCESS;
}
diff --git a/t/test_classify_ref.cc b/t/test_classify_ref.cc
index e9b0067..a53a0a3 100644
--- a/t/test_classify_ref.cc
+++ b/t/test_classify_ref.cc
@@ -1,75 +1,95 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <fstream>
#include <random>
#include <algorithm>
#include "LHEF/LHEF.h"
#include "HEJ/YAMLreader.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReader.hh"
namespace{
// this is deliberately chosen bigger than in the generation,
// to cluster multiple partons in one jet
constexpr double min_jet_pt = 40.;
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
void shuffle_particles(HEJ::Event::EventData & ev) {
static std::mt19937_64 ran{0};
+ // incoming
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
- std::shuffle(begin(ev.outgoing), end(ev.outgoing), ran);
+ // outgoing (through index)
+ auto old_outgoing = std::move(ev.outgoing);
+ std::vector<size_t> idx(old_outgoing.size());
+ std::iota(idx.begin(), idx.end(), 0);
+ std::shuffle(begin(idx), end(idx), ran);
+ ev.outgoing.clear();
+ ev.outgoing.reserve(old_outgoing.size());
+ for(size_t i: idx) {
+ ev.outgoing.emplace_back(std::move(old_outgoing[i]));
+ }
+ // find decays again
+ if(!ev.decays.empty()){
+ auto old_decays = std::move(ev.decays);
+ ev.decays.clear();
+ for(size_t i=0; i<idx.size(); ++i) {
+ auto decay = old_decays.find(idx[i]);
+ if(decay != old_decays.end())
+ ev.decays.emplace(i, std::move(decay->second));
+ }
+ }
}
}
int main(int argn, char** argv) {
if(argn != 3 && argn != 4){
std::cerr << "Usage: " << argv[0]
<< " reference_classification input_file.lhe\n";
return EXIT_FAILURE;
}
bool OUTPUT_MODE = false;
if(argn == 4 && std::string("OUTPUT")==std::string(argv[3]))
OUTPUT_MODE = true;
std::fstream ref_file;
if ( OUTPUT_MODE ) {
std::cout << "_______________________USING OUTPUT MODE!_______________________" << std::endl;
ref_file.open(argv[1], std::fstream::out);
} else {
ref_file.open(argv[1], std::fstream::in);
}
auto reader{ HEJ::make_reader(argv[2]) };
while(reader->read_event()){
HEJ::Event::EventData data{ reader->hepeup() };
shuffle_particles(data);
const HEJ::Event ev{
data.cluster(
jet_def, min_jet_pt
)
};
if ( OUTPUT_MODE ) {
ref_file << ev.type() << std::endl;
} else {
std::string line;
if(!std::getline(ref_file,line)) break;
const auto expected{static_cast<HEJ::event_type::EventType>(std::stoi(line))};
if(ev.type() != expected){
using HEJ::event_type::name;
std::cerr << "wrong classification of event:\n" << ev
<< "classified as " << name(ev.type())
<< ", is " << name(expected) << '\n';
return EXIT_FAILURE;
}
}
}
}
diff --git a/t/test_colours.cc b/t/test_colours.cc
index d47b225..e85ae4c 100644
--- a/t/test_colours.cc
+++ b/t/test_colours.cc
@@ -1,281 +1,301 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <random>
#include <stdexcept>
#include <utility>
#include "HEJ/Event.hh"
#include "HEJ/RNG.hh"
#define ASSERT(x) if(!(x)) { \
throw std::logic_error("Assertion '" #x "' failed."); \
}
/// biased RNG to connect always to colour
class dum_rnd: public HEJ::DefaultRNG {
public:
dum_rnd() = default;
double flat() override {
return 0.;
};
};
void shuffle_particles(HEJ::Event::EventData & ev) {
static std::mt19937_64 ran{0};
+ // incoming
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
- std::shuffle(begin(ev.outgoing), end(ev.outgoing), ran);
+ // outgoing (through index)
+ auto old_outgoing = std::move(ev.outgoing);
+ std::vector<size_t> idx(old_outgoing.size());
+ std::iota(idx.begin(), idx.end(), 0);
+ std::shuffle(begin(idx), end(idx), ran);
+ ev.outgoing.clear();
+ ev.outgoing.reserve(old_outgoing.size());
+ for(size_t i: idx) {
+ ev.outgoing.emplace_back(std::move(old_outgoing[i]));
+ }
+ // find decays again
+ if(!ev.decays.empty()){
+ auto old_decays = std::move(ev.decays);
+ ev.decays.clear();
+ for(size_t i=0; i<idx.size(); ++i) {
+ auto decay = old_decays.find(idx[i]);
+ if(decay != old_decays.end())
+ ev.decays.emplace(i, std::move(decay->second));
+ }
+ }
}
void dump_event(HEJ::Event const & ev){
for(auto const & in: ev.incoming()){
std::cerr << "in type=" << in.type
<< ", colour={" << (*in.colour).first
<< ", " << (*in.colour).second << "}\n";
}
for(auto const & out: ev.outgoing()){
std::cerr << "out type=" << out.type << ", colour={";
if(out.colour)
std::cerr << (*out.colour).first << ", " << (*out.colour).second;
else
std::cerr << "non, non";
std::cerr << "}\n";
}
}
/// true if colour is allowed for particle
bool correct_colour(HEJ::Particle const & part){
if(HEJ::is_AWZH_boson(part) && !part.colour) return true;
if(!part.colour) return false;
int const colour = part.colour->first;
int const anti_colour = part.colour->second;
if(part.type == HEJ::ParticleID::gluon)
return colour != anti_colour && colour > 0 && anti_colour > 0;
if(HEJ::is_quark(part))
return anti_colour == 0 && colour > 0;
return colour == 0 && anti_colour > 0;
}
bool correct_colour(HEJ::Event const & ev){
for(auto const & part: ev.incoming()){
if(!correct_colour(part))
return false;
}
for(auto const & part: ev.outgoing()){
if(!correct_colour(part))
return false;
}
return true;
}
bool match_expected(
HEJ::Event const & ev,
std::vector<HEJ::Colour> const & expected
){
ASSERT(ev.outgoing().size()+2==expected.size());
for(size_t i=0; i<ev.incoming().size(); ++i){
ASSERT(ev.incoming()[i].colour);
if( *ev.incoming()[i].colour != expected[i])
return false;
}
for(size_t i=2; i<ev.outgoing().size()+2; ++i){
if( ev.outgoing()[i-2].colour ){
if( *ev.outgoing()[i-2].colour != expected[i] )
return false;
} else if( expected[i].first != 0 || expected[i].second != 0)
return false;
}
return true;
}
void check_event(
HEJ::Event::EventData unc_ev, std::vector<HEJ::Colour> const & expected_colours
){
shuffle_particles(unc_ev); // make sure incoming order doesn't matter
HEJ::Event ev{unc_ev.cluster(
fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30.)
};
ASSERT(HEJ::event_type::is_resummable(ev.type()));
dum_rnd rng;
ASSERT(ev.generate_colours(rng));
if(!correct_colour(ev)){
std::cerr << "Found illegal colours for event\n";
dump_event(ev);
throw std::invalid_argument("Illegal colour set");
}
if(!match_expected(ev, expected_colours)){
std::cerr << "Colours didn't match expectation. Found\n";
dump_event(ev);
std::cerr << "but expected\n";
for(auto const & col: expected_colours){
std::cerr << "colour={" << col.first << ", " << col.second << "}\n";
}
throw std::logic_error("Colours did not match expectation");
}
}
int main() {
HEJ::Event::EventData ev;
std::vector<HEJ::Colour> expected_colours(7);
/// pure gluon
ev.incoming[0] = { HEJ::ParticleID::gluon, { 0, 0,-427, 427}, {}};
ev.incoming[1] = { HEJ::ParticleID::gluon, { 0, 0, 851, 851}, {}};
ev.outgoing.push_back({ HEJ::ParticleID::gluon, { 196, 124, -82, 246}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-167,-184, 16, 249}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::higgs, { 197, 180, 168, 339}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-190, -57, 126, 235}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::gluon, { -36, -63, 196, 209}, {}});
expected_colours[0] = {502, 501};
expected_colours[1] = {509, 502};
expected_colours[2] = {503, 501};
expected_colours[3] = {505, 503};
expected_colours[4] = { 0, 0};
expected_colours[5] = {507, 505};
expected_colours[6] = {509, 507};
check_event(ev, expected_colours);
/// last g to Qx (=> gQx -> g ... Qx)
ev.incoming[1].type = HEJ::ParticleID::d_bar;
ev.outgoing[4].type = HEJ::ParticleID::d_bar;
// => only end changes
expected_colours[1].first = 0;
expected_colours[6].first = 0;
check_event(ev, expected_colours);
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno forward (=> gQx -> g ... Qx g)
std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
// => uno quarks eats colour and gluon connects to anti-colour
new_expected[5] = {0, expected_colours[3].first};
new_expected[6] = {expected_colours[0].first, expected_colours[0].first+2};
new_expected[1].second += 2; // one more anti-colour in line
check_event(new_ev, new_expected);
}
/// swap Qx <-> Q (=> gQ -> g ... Q)
ev.incoming[1].type = HEJ::ParticleID::d;
ev.outgoing[4].type = HEJ::ParticleID::d;
// => swap: colour<->anti && inital<->final
std::swap(expected_colours[1], expected_colours[6]);
std::swap(expected_colours[1].first, expected_colours[1].second);
std::swap(expected_colours[6].first, expected_colours[6].second);
check_event(ev, expected_colours);
/// first g to qx (=> qxQ -> qx ... Q)
ev.incoming[0].type = HEJ::ParticleID::u_bar;
ev.outgoing[0].type = HEJ::ParticleID::u_bar;
expected_colours[0] = { 0, 501};
// => shift anti-colour index one up
expected_colours[1].first -= 2;
expected_colours[5] = expected_colours[3];
expected_colours[3] = expected_colours[2];
expected_colours[2] = { 0, 502};
check_event(ev, expected_colours);
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno backward (=> qxQ -> g qx ... Q)
std::swap(new_ev.outgoing[0].type, new_ev.outgoing[1].type);
// => uno gluon connects to quark colour
new_expected[3] = expected_colours[2];
new_expected[2] = {expected_colours[0].second+2, expected_colours[0].second};
check_event(new_ev, new_expected);
/// swap qx <-> q (=> qQ -> g q ... Q)
new_ev.incoming[0].type = HEJ::ParticleID::u;
new_ev.outgoing[1].type = HEJ::ParticleID::u;
// => swap: colour<->anti && inital<->final
std::swap(new_expected[0], new_expected[3]);
std::swap(new_expected[0].first, new_expected[0].second);
std::swap(new_expected[3].first, new_expected[3].second);
// => & connect first gluon with remaining anti-colour
new_expected[2] = {new_expected[0].first, new_expected[0].first+2};
// shift colour line one down
new_expected[1].first-=2;
new_expected[5].first-=2;
new_expected[5].second-=2;
// shift anti-colour line one up
new_expected[6].first+=2;
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno forward (=> qxQ -> qx ... Q g)
std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
// => uno gluon connects to remaining colour
new_expected[5] = expected_colours[6];
new_expected[6] = {expected_colours[3].first+2, expected_colours[3].first};
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// qqx backward (=> gQ -> qx q ... Q)
// => swap: incoming q <-> outgoing gluon
std::swap(new_ev.incoming[0].type, new_ev.outgoing[1].type);
new_ev.outgoing[1].type=static_cast<HEJ::ParticleID>(
-1*new_ev.outgoing[1].type);
// incoming q -> outgoing q (colour<->anti)
std::swap(new_expected[0], new_expected[3]);
std::swap(new_expected[3].first, new_expected[3].second);
new_expected[3].first+=2;
new_expected[0].first-=1; // skip one index
// couple first in to first out
new_expected[2].second=new_expected[0].second;
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// qqx forward (=> qx g -> qx ... Qx Q)
// => swap: incoming Q <-> outgoing gluon
std::swap(new_ev.incoming[1].type, new_ev.outgoing[3].type);
new_ev.outgoing[3].type=static_cast<HEJ::ParticleID>(
-1*new_ev.outgoing[3].type);
// incoming q -> outgoing q (colour<->anti)
std::swap(new_expected[1], new_expected[5]);
std::swap(new_expected[5].first, new_expected[5].second);
new_expected[5].second-=2;
new_expected[1].second-=1; // skip one index
// couple last in to last out
new_expected[6].first=new_expected[1].first;
check_event(new_ev, new_expected);
// move Higgs to position 1 (=> qx g -> qx h g Qx Q)
std::swap(new_ev.outgoing[1].type, new_ev.outgoing[2].type);
std::swap(new_expected[3], new_expected[4]); // trivial
check_event(new_ev, new_expected);
// central qqx (=> qx g -> qx h Q Qx g)
// => swap: Q <-> g
std::swap(new_ev.outgoing[2].type, new_ev.outgoing[4].type);
std::swap(new_expected[4], new_expected[6]);
// gluon was connected on left side, i.e. doesn't matter for QQx
// => couple Q to out qx
new_expected[4].first = new_expected[2].second;
// Qx next in line
new_expected[5].second = new_expected[4].first+2;
// incoming g shifted by one position in line
new_expected[1].first-=2;
new_expected[1].second+=2;
check_event(new_ev, new_expected);
}
return EXIT_SUCCESS;
}
diff --git a/t/test_scale_arithmetics.cc b/t/test_scale_arithmetics.cc
index ff27844..3a66bab 100644
--- a/t/test_scale_arithmetics.cc
+++ b/t/test_scale_arithmetics.cc
@@ -1,100 +1,120 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <fstream>
#include <random>
#include <algorithm>
#include "LHEF/LHEF.h"
#include "HEJ/EventReweighter.hh"
#include "HEJ/make_RNG.hh"
#include "HEJ/Event.hh"
#include "HEJ/YAMLreader.hh"
#include "HEJ/stream.hh"
constexpr double alpha_s = 0.118;
constexpr double ep = 1e-13;
void dump(HEJ::Event const & ev){
{
LHEF::Writer writer{std::cout};
std::cout << std::setprecision(6);
writer.hepeup = to_HEPEUP(std::move(ev), nullptr);
writer.writeEvent();
}
std::cout << "Rapidity ordering:\n";
for(const auto & part: ev.outgoing()){
std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl;
}
}
void shuffle_particles(HEJ::Event::EventData & ev) {
static std::mt19937_64 ran{0};
+ // incoming
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
- std::shuffle(begin(ev.outgoing), end(ev.outgoing), ran);
+ // outgoing (through index)
+ auto old_outgoing = std::move(ev.outgoing);
+ std::vector<size_t> idx(old_outgoing.size());
+ std::iota(idx.begin(), idx.end(), 0);
+ std::shuffle(begin(idx), end(idx), ran);
+ ev.outgoing.clear();
+ ev.outgoing.reserve(old_outgoing.size());
+ for(size_t i: idx) {
+ ev.outgoing.emplace_back(std::move(old_outgoing[i]));
+ }
+ // find decays again
+ if(!ev.decays.empty()){
+ auto old_decays = std::move(ev.decays);
+ ev.decays.clear();
+ for(size_t i=0; i<idx.size(); ++i) {
+ auto decay = old_decays.find(idx[i]);
+ if(decay != old_decays.end())
+ ev.decays.emplace(i, std::move(decay->second));
+ }
+ }
}
int main(int argn, char** argv){
if(argn != 3){
std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml input_file.lhe\n\n";
return EXIT_FAILURE;
}
HEJ::Config config = HEJ::load_config(argv[1]);
config.scales = HEJ::to_ScaleConfig(
YAML::Load("scales: [H_T, 1 * H_T, 2/2 * H_T, 2*H_T/2, H_T/2*2, H_T/2/2*4, H_T*H_T/H_T]")
);
HEJ::istream in{argv[2]};
LHEF::Reader reader{in};
auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);
HEJ::ScaleGenerator scale_gen{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
};
HEJ::EventReweighter resum{
reader.heprup,
std::move(scale_gen),
to_EventReweighterConfig(config),
*ran
};
size_t i = 0;
while(reader.readEvent()){
++i;
HEJ::Event::EventData data{reader.hepeup};
shuffle_particles(data);
HEJ::Event event{
data.cluster(
config.resummation_jets.def,
config.resummation_jets.min_pt
)
};
auto resummed = resum.reweight(event, config.trials);
for(auto && ev: resummed) {
for(auto &&var: ev.variations()) {
if(std::abs(var.muf - ev.central().muf) > ep) {
std::cerr
<< std::setprecision(15)
<< "unequal scales: " << var.muf
<< " != " << ev.central().muf << '\n'
<< "in resummed event:\n";
dump(ev);
std::cerr << "\noriginal event:\n";
dump(event);
return EXIT_FAILURE;
}
}
}
}
}