Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8725449
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
61 KB
Subscribers
None
View Options
diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
index e7111fe..0e31d49 100644
--- a/FixedOrderGen/src/PhaseSpacePoint.cc
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -1,471 +1,476 @@
#include "PhaseSpacePoint.hh"
#include <random>
#include "RHEJ/kinematics.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/exceptions.hh"
#include "Process.hh"
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
using namespace RHEJ;
namespace HEJFOG{
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
RHEJ::UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
RHEJ::UnclusteredEvent result;
result.incoming = psp.incoming();
std::sort(
begin(result.incoming), end(result.incoming),
[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
);
assert(result.incoming.size() == 2);
result.outgoing = psp.outgoing();
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
RHEJ::rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.decays = psp.decays();
result.central.mur = NaN;
result.central.muf = NaN;
result.central.weight = psp.weight();
return result;
}
namespace{
bool can_swap_to_uno(
RHEJ::Particle const & p1, RHEJ::Particle const & p2
){
return is_parton(p1)
&& p1.type != RHEJ::pid::gluon
&& p2.type == RHEJ::pid::gluon;
}
}
void PhaseSpacePoint::maybe_turn_to_uno(
double chance,
RHEJ::RNG & ran
){
assert(outgoing_.size() >= 2);
const size_t nout = outgoing_.size();
const bool can_be_uno_backward = can_swap_to_uno(
outgoing_[0], outgoing_[1]
);
const bool can_be_uno_forward = can_swap_to_uno(
outgoing_[nout-1], outgoing_[nout-2]
);
if(!can_be_uno_backward && !can_be_uno_forward) return;
if(ran.flat() < chance){
weight_ /= chance;
if(can_be_uno_backward && can_be_uno_forward){
if(ran.flat() < 0.5){
std::swap(outgoing_[0].type, outgoing_[1].type);
}
else{
std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
}
weight_ *= 2.;
}
else if(can_be_uno_backward){
std::swap(outgoing_[0].type, outgoing_[1].type);
}
else{
assert(can_be_uno_forward);
std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
}
}
else weight_ /= 1 - chance;
}
PhaseSpacePoint::PhaseSpacePoint(
Process const & proc,
JetParameters const & jet_param,
RHEJ::PDF & pdf, double E_beam,
double uno_chance,
HiggsProperties const & h,
RHEJ::RNG & ran
)
{
assert(proc.njets >= 2);
const int nout = proc.njets + (proc.boson?1:0);
status_ = good;
weight_ = 1;
weight_ /= std::tgamma(nout);
outgoing_.reserve(nout);
for(auto&& p_out: gen_LO_partons(proc.njets, jet_param, E_beam, ran)){
outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out)});
}
if(status_ != good) return;
if(proc.boson && *proc.boson == pid::Higgs){
// The Higgs
auto Hparticle=gen_boson(pid::higgs, h.mass, h.width, ran);
auto pos=std::upper_bound(
begin(outgoing_),end(outgoing_),Hparticle,rapidity_less{}
);
outgoing_.insert(pos, Hparticle);
if(! h.decays.empty()){
const int boson_idx = std::distance(begin(outgoing_), pos);
decays_.emplace(
boson_idx,
decay_boson(outgoing_[boson_idx], h.decays, ran)
);
}
}
+ // normalisation of momentum-conserving delta function
+ weight_ *= pow(2*M_PI, 4);
reconstruct_incoming(proc.incoming, pdf, E_beam, jet_param.min_pt, ran);
if(status_ != good) return;
// set outgoing states
most_backward_FKL(outgoing_).type = incoming_[0].type;
most_forward_FKL(outgoing_).type = incoming_[1].type;
if(proc.njets > 2) maybe_turn_to_uno(uno_chance, ran);
}
double PhaseSpacePoint::gen_hard_pt(
int np , double ptmin, double ptmax, double y,
RHEJ::RNG & ran
) {
// heuristic parameters for pt sampling
const double ptpar = ptmin + np/5.;
const double arg_small_y = atan((ptmax - ptmin)/ptpar);
const double y_cut = 3.;
const double r1 = ran.flat();
if(y < y_cut){
const double pt = ptmin + ptpar*tan(r1*arg_small_y);
const double temp = cos(r1*arg_small_y);
weight_ *= pt*ptpar*arg_small_y/(temp*temp);
return pt;
}
const double ptpar2 = ptpar/(1 + 5*(y-y_cut));
const double temp = 1. - std::exp((ptmin-ptmax)/ptpar2);
const double pt = ptmin - ptpar2*std::log(1-r1*temp);
weight_ *= pt*ptpar2*temp/(1-r1*temp);
return pt;
}
double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, RHEJ::RNG & ran) {
constexpr double ptpar = 4.;
const double r = ran.flat();
const double pt = max_pt + ptpar/np*std::log(r);
weight_ *= pt*ptpar/(np*r);
return pt;
}
double PhaseSpacePoint::gen_parton_pt(
int count, JetParameters const & jet_param, double max_pt, double y,
RHEJ::RNG & ran
) {
constexpr double p_small_pt = 0.02;
if(! jet_param.peak_pt) {
return gen_hard_pt(count, jet_param.min_pt, max_pt, y, ran);
}
const double r = ran.flat();
if(r > p_small_pt) {
weight_ /= 1. - p_small_pt;
return gen_hard_pt(count, *jet_param.peak_pt, max_pt, y, ran);
}
weight_ /= p_small_pt;
const double pt = gen_soft_pt(count, *jet_param.peak_pt, ran);
if(pt < jet_param.min_pt) {
weight_=0.0;
status_ = not_enough_jets;
return jet_param.min_pt;
}
return pt;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_LO_partons(
int np ,
JetParameters const & jet_param,
double max_pt,
RHEJ::RNG & ran
){
if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"};
- weight_ /= pow(16.*pow(M_PI,3),np-2);
+ weight_ /= pow(16.*pow(M_PI,3),np);
std::vector<fastjet::PseudoJet> partons;
partons.reserve(np);
for(int i = 0; i < np; ++i){
const double y = -jet_param.max_y + 2*jet_param.max_y*ran.flat();
weight_ *= 2*jet_param.max_y;
const double phi = 2*M_PI*ran.flat();
weight_ *= 2.0*M_PI;
const double pt = gen_parton_pt(np, jet_param, max_pt, y, ran);
if(weight_ == 0.0) return {};
partons.emplace_back(fastjet::PtYPhiM(pt, y, phi));
assert(jet_param.min_pt <= partons[i].pt());
assert(partons[i].pt() <= max_pt+1e-5);
}
// Need to check that at LO, the number of jets = number of partons;
fastjet::ClusterSequence cs(partons, jet_param.def);
auto cluster_jets=cs.inclusive_jets(jet_param.min_pt);
if (cluster_jets.size()!=unsigned(np)){
weight_=0.0;
status_ = not_enough_jets;
return {};
}
std::sort(begin(partons), end(partons), rapidity_less{});
return partons;
}
Particle PhaseSpacePoint::gen_boson(
RHEJ::ParticleID bosonid, double mass, double width,
RHEJ::RNG & ran
){
- // The Higgs:
+ // Usual phase space measure
+ weight_ /= 16.*pow(M_PI, 3);
+
// Generate a y Gaussian distributed around 0
+ // TODO: magic number only for Higgs
const double y = random_normal(1.6, ran);
const double r1 = ran.flat();
const double sH = mass*(
mass + width*tan(M_PI/2.*r1 + (r1-1.)*atan(mass/width))
);
auto p = gen_last_momentum(sH, y);
return Particle{bosonid, std::move(p)};
}
fastjet::PseudoJet PhaseSpacePoint::gen_last_momentum(
const double mass_square, const double y
) const {
std::array<double,2> pt{0.,0.};
for (auto const & parton:outgoing_) {
pt[0]-=parton.px();
pt[1]-=parton.py();
}
const double mperp = sqrt(pt[0]*pt[0]+pt[1]*pt[1]+mass_square);
const double pz=mperp*sinh(y);
const double E=mperp*cosh(y);
return {pt[0], pt[1], pz, E};
}
Particle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> const & partons
) const{
if(!RHEJ::is_parton(partons[0])) return partons[1];
return partons[0];
}
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
const size_t last_idx = partons.size() - 1;
if(!RHEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
return partons[last_idx];
}
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
if(!RHEJ::is_parton(partons[0])) return partons[1];
return partons[0];
}
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
const size_t last_idx = partons.size() - 1;
if(!RHEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
return partons[last_idx];
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<RHEJ::ParticleID, 2> const & ids,
RHEJ::PDF & pdf, double E_beam,
double uf,
RHEJ::RNG & ran
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
// calculate xa, xb
const double sqrts=2*E_beam;
const double xa=(incoming_[0].p.e()-incoming_[0].p.pz())/sqrts;
const double xb=(incoming_[1].p.e()+incoming_[1].p.pz())/sqrts;
// abort if phase space point is outside of collider energy reach
if (xa>1. || xb>1.){
weight_=0;
status_ = too_much_energy;
return;
}
// pick pdfs
/** TODO:
* ufa, ufb don't correspond to our final scale choice.
* The reversed HEJ scale generators currently expect a full event as input,
* so fixing this is not completely trivial
*/
if(ids[0] == pid::proton || ids[0] == pid::p_bar){
const double ufa=uf;
incoming_[0].type = generate_incoming_id(xa, ufa, pdf, ran);
}
else {
incoming_[0].type = ids[0];
}
if(ids[1] == pid::proton || ids[1] == pid::p_bar){
const double ufb=uf;
incoming_[1].type = generate_incoming_id(xb, ufb, pdf, ran);
}
else {
incoming_[1].type = ids[1];
}
assert(momentum_conserved(1e-7));
}
RHEJ::ParticleID PhaseSpacePoint::generate_incoming_id(
double x, double uf, RHEJ::PDF & pdf,
RHEJ::RNG & ran
){
const double pdfg=fabs(pdf.pdfpt(0,x,uf,pid::gluon));
const double pdfu=fabs(pdf.pdfpt(0,x,uf,pid::up));
const double pdfd=fabs(pdf.pdfpt(0,x,uf,pid::down));
const double pdfux=fabs(pdf.pdfpt(0,x,uf,pid::u_bar));
const double pdfdx=fabs(pdf.pdfpt(0,x,uf,pid::d_bar));
const double pdfc=fabs(pdf.pdfpt(0,x,uf,pid::charm));
const double pdfs=fabs(pdf.pdfpt(0,x,uf,pid::strange));
const double pdfsx=fabs(pdf.pdfpt(0,x,uf,pid::s_bar));
const double pdfb=fabs(pdf.pdfpt(0,x,uf,pid::b));
const double pdftot=pdfg+4.0/9.0*(pdfu + pdfd + pdfux + pdfdx +pdfs +pdfsx + 2.0*(pdfc +pdfb ));
const double r1=pdftot*ran.flat();
double sum;
if (r1<(sum=pdfg)) {
weight_*=pdftot/pdfg;
return pid::gluon;
}
if (r1<(sum+=(4./9.)*pdfu)) {
weight_*=pdftot/((4./9.)*pdfu);
return pid::up;
}
else if (r1<(sum+=(4./9.)*pdfd)) {
weight_*=pdftot/((4./9.)*pdfd);
return pid::down;
}
else if (r1<(sum+=(4./9.)*pdfux)) {
weight_*=pdftot/((4./9.)*pdfux);
return pid::u_bar;
}
else if (r1<(sum+=(4./9.)*pdfdx)) {
weight_*=pdftot/((4./9.)*pdfdx);
return pid::d_bar;
}
else if (r1<(sum+=(4./9.)*pdfc)) {
weight_*=pdftot/((4./9.)*pdfc);
return pid::c;
}
else if (r1<(sum+=(4./9.)*pdfc)){
weight_*=pdftot/((4./9.)*pdfc);
return pid::c_bar;
}
else if (r1<(sum+=(4./9.)*pdfs)) {
weight_*=pdftot/((4./9.)*pdfs);
return pid::s;
}
else if (r1<(sum+=(4./9.)*pdfsx)) {
weight_*=pdftot/((4./9.)*pdfsx);
return pid::s_bar;
}
else if (r1<(sum+=(4./9.)*pdfb)) {
weight_*=pdftot/((4./9.)*pdfb);
return pid::b;
}
else if (r1<=(sum+=(4./9.)*pdfb)) {
weight_*=pdftot/((4./9.)*pdfb);
return pid::b_bar;
}
std::cout << "Error in choosing incoming parton: "<<x<<" "<<uf<<" "<<sum<<" "<<pdftot<<" "<<r1;
std::cout << " "<<pdfg+4./9.*(pdfu+pdfux+pdfd+pdfdx+pdfs+pdfsx+2.*(pdfc+pdfb))<<std::endl;
throw std::logic_error{"Failed to choose parton flavour"};
}
double PhaseSpacePoint::random_normal(
double stddev,
RHEJ::RNG & ran
){
const double r1 = ran.flat();
const double r2 = ran.flat();
const double lninvr1 = -log(r1);
const double result = stddev*sqrt(2.*lninvr1)*cos(2.*M_PI*r2);
weight_ *= exp(result*result/(2*stddev*stddev))*sqrt(2.*M_PI)*stddev;
return result;
}
bool PhaseSpacePoint::momentum_conserved(double ep) const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
for(auto const & out: outgoing()) diff -= out.p;
return nearby_ep(diff, fastjet::PseudoJet{}, ep);
}
Decay PhaseSpacePoint::select_decay_channel(
std::vector<Decay> const & decays,
RHEJ::RNG & ran
){
double br_total = 0.;
for(auto const & decay: decays) br_total += decay.branching_ratio;
// adjust weight
// this is given by (channel branching ratio)/(chance to pick channel)
// where (chance to pick channel) =
// (channel branching ratio)/(total branching ratio)
weight_ *= br_total;
const double r1 = br_total*ran.flat();
double br_sum = 0.;
for(auto const & decay: decays){
br_sum += decay.branching_ratio;
if(r1 < br_sum) return decay;
}
throw std::logic_error{"unreachable"};
}
std::vector<Particle> PhaseSpacePoint::decay_boson(
RHEJ::Particle const & parent,
std::vector<Decay> const & decays,
RHEJ::RNG & ran
){
const auto channel = select_decay_channel(decays, ran);
if(channel.products.size() != 2){
throw RHEJ::not_implemented{
"only decays into two particles are implemented"
};
}
std::vector<Particle> decay_products(channel.products.size());
for(size_t i = 0; i < channel.products.size(); ++i){
decay_products[i].type = channel.products[i];
}
// choose polar and azimuth angle in parent rest frame
const double E = parent.m()/2;
const double theta = 2.*M_PI*ran.flat();
const double cos_phi = 2.*ran.flat()-1.;
const double sin_phi = sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
const double px = E*cos(theta)*sin_phi;
const double py = E*sin(theta)*sin_phi;
const double pz = E*cos_phi;
decay_products[0].p.reset(px, py, pz, E);
decay_products[1].p.reset(-px, -py, -pz, E);
for(auto & particle: decay_products) particle.p.boost(parent.p);
return decay_products;
}
}
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 637a262..fff1079 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,761 +1,761 @@
#include "RHEJ/MatrixElement.hh"
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
#include "RHEJ/Constants.hh"
#include "RHEJ/currents.hh"
#include "RHEJ/PDG_codes.hh"
#include "RHEJ/uno.hh"
#include "RHEJ/utility.hh"
namespace RHEJ{
//cf. last line of eq. (22) in \ref Andersen:2011hs
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j, double lambda
) const {
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;
}
double MatrixElement::virtual_corrections(
double mur,
std::array<Particle, 2> const & in,
std::vector<Particle> const & out
) const{
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(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 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 || has_unob_gluon(in, out)){
q -= out[1].p;
++first_idx;
}
if(out.back().type == pid::Higgs || has_unof_gluon(in, out)){
--last_idx;
}
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, CLAMBDA)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| out.back().type == pid::Higgs
|| has_unof_gluon(in, out)
);
return exp(exponent);
}
} // namespace RHEJ
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
// cout << "#Fadin qa : "<<qav<<endl;
// cout << "#Fadin qb : "<<qbv<<endl;
// cout << "#Fadin p1 : "<<p1<<endl;
// cout << "#Fadin p2 : "<<p2<<endl;
// cout << "#Fadin p5 : "<<p5<<endl;
// cout << "#Fadin Gauge Check : "<< CL.dot(p5)<<endl;
// cout << "#Fadin C2L : "<< -CL.dot(CL)<<" "<<-CL.dot(CL)/(qav.m2()*qbv.m2())/(4./p5.perp2())<<endl;
// TODO can this dead test go?
// if (-CL.dot(CL)<0.)
// if (fabs(CL.dot(p5))>fabs(CL.dot(CL))) // not sufficient!
// return 0.;
// else
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
double kperp=(qav-qbv).perp();
if (kperp>RHEJ::CLAMBDA)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
return Cls-4./(kperp*kperp);
}
}
//! Lipatov vertex
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip,
CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
double kperp=(qav-qbv).perp();
if (kperp>RHEJ::CLAMBDA)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
double temp=Cls-4./(kperp*kperp);
return temp;
}
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
if (aptype==21&&bptype==21) {
return jM2gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2qg(pn,pb,p1,pa);
else
return jM2qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return jM2qg(p1,pa,pn,pb);
else
return jM2qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2qQ(pn,pb,p1,pa);
else
return jM2qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return jM2qQbar(p1,pa,pn,pb);
else
return jM2qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Matrix element squared for tree-level current-current scattering with Higgs
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
*/
double ME_Higgs_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype==21) // gg initial state
return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered forward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param punof Unordered Particle Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered forward emission
*/
double ME_Higgs_current_unof(
int aptype, int bptype,
CLHEP::HepLorentzVector const & punof,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (aptype>0)
return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param punob Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*/
double ME_Higgs_current_unob(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & punob,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (bptype==21&&aptype!=21) {
if (aptype > 0)
return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (bptype>0)
return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
CLHEP::HepLorentzVector to_HepLorentzVector(RHEJ::Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
} // namespace anonymous
namespace RHEJ{
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{}
double MatrixElement::operator()(
double mur,
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
return tree(
mur,
incoming, outgoing,
check_momenta
)*virtual_corrections(
mur,
incoming, outgoing
);
}
double MatrixElement::tree_kin(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
assert(
std::is_sorted(
incoming.begin(), incoming.end(),
[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
)
);
assert(std::is_sorted(outgoing.begin(), outgoing.end(), rapidity_less{}));
auto AWZH_boson = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(outgoing)){
return tree_kin_jets(incoming, outgoing, check_momenta);
}
switch(AWZH_boson->type){
case pid::Higgs: {
static constexpr double mH = 125.;
const double alpha_s_mH = alpha_s_(mH);
- return alpha_s_mH*alpha_s_mH/(256.*pow(M_PI, 5))*tree_kin_Higgs(
+ return alpha_s_mH*alpha_s_mH*tree_kin_Higgs(
incoming, outgoing, check_momenta
);
}
// TODO
case pid::Wp:
case pid::Wm:
case pid::photon:
case pid::Z:
default:
throw std::logic_error("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 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)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
std::array<Particle, 2> const & incoming,
std::vector<Particle> out_partons, bool check_momenta
) const{
if(!check_momenta){
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
fastjet::ClusterSequence cs(to_PseudoJet(out_partons), param_.jet_param.def);
const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission
if(has_unob_gluon(incoming, out_partons)){
assert(jets.size() >= 3);
++most_backward;
}
else if(has_unof_gluon(incoming, out_partons)){
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(RHEJ::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(
std::array<Particle, 2> const & incoming,
std::vector<Particle> partons,
bool check_momenta
) const {
partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
if(has_unob_gluon(incoming, partons) || has_unof_gluon(incoming, partons)){
throw std::logic_error("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
);
}
double MatrixElement::tree_kin_Higgs(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
if(has_uno_gluon(incoming, outgoing)){
return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
}
if(outgoing.front().type == pid::Higgs){
return tree_kin_Higgs_first(incoming, outgoing, check_momenta);
}
if(outgoing.back().type == pid::Higgs){
return tree_kin_Higgs_last(incoming, outgoing, check_momenta);
}
return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
}
double MatrixElement::MH2_forwardH(
RHEJ::ParticleID id,
CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
assert(RHEJ::is_parton(id));
if(id != RHEJ::pid::gluon){
return 9./2.*shat*shat*C2qHqm(p1in,p1out,pH)/(t1*t2);
}
// gluon case
#ifdef RHEJ_BUILD_WITH_QCDLOOP
if(!param_.Higgs_coupling.use_impact_factors){
return C_A/C_F*1./(16*M_PI*M_PI)*t1/t2*MH2gq_outsideH(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb
);
}
#endif
return 9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
assert(outgoing.front().type == pid::Higgs);
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
incoming,
std::vector<Particle>(begin(outgoing) + 1, end(outgoing)),
check_momenta
);
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();
double wt = MH2_forwardH(
incoming[0].type, p1, pa, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn
);
for(auto const & inc: incoming){
if(inc.type != pid::gluon) wt *= C_F/C_A;
}
return wt;
}
double MatrixElement::tree_kin_Higgs_last(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
assert(outgoing.back().type == pid::Higgs);
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
incoming,
std::vector<Particle>(begin(outgoing), end(outgoing) - 1),
check_momenta
);
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();
double wt = MH2_forwardH(
incoming[1].type, pn, pb, p1, pa, pH,
t2, t1
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn
);
for(auto const & inc: incoming){
if(inc.type != pid::gluon) wt *= C_F/C_A;
}
return wt;
}
double MatrixElement::tree_kin_Higgs_between(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
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);
std::vector<Particle> partons(begin(outgoing), the_Higgs);
partons.insert(end(partons), the_Higgs + 1, end(outgoing));
partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[has_unob_gluon(incoming, outgoing)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - (has_unof_gluon(incoming, outgoing)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && has_unob_gluon(incoming, outgoing))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && has_unof_gluon(incoming, outgoing))
|| (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(has_unob_gluon(incoming, outgoing)){
current_factor = 9./2.*ME_Higgs_current_unob(
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(has_unof_gluon(incoming, outgoing)){
current_factor = 9./2.*ME_Higgs_current_unof(
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
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn
);
return current_factor*9./8.*ladder_factor;
}
double MatrixElement::tree_param_partons(
double alpha_s, double mur,
std::vector<Particle> const & partons
) const{
const double gs2 = 4.*M_PI*alpha_s;
double wt = std::pow(gs2, partons.size());
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(partons.size() >= 2);
for(size_t i = 1; i < partons.size()-1; ++i){
wt *= 1 + alpha_s/(2*M_PI)*beta0*log(mur/partons[i].p.perp());
}
}
return wt;
}
double MatrixElement::tree_param(
double mur,
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing
) const{
const double alpha_s = alpha_s_(mur);
if(has_unob_gluon(incoming, outgoing)){
return 4*M_PI*alpha_s*tree_param_partons(
alpha_s, mur, filter_partons({begin(outgoing) + 1, end(outgoing)})
);
}
if(has_unof_gluon(incoming, outgoing)){
return 4*M_PI*alpha_s*tree_param_partons(
alpha_s, mur, filter_partons({begin(outgoing), end(outgoing) - 1})
);
}
return tree_param_partons(alpha_s, mur, filter_partons(outgoing));
}
double MatrixElement::tree(
double mur,
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
return tree_param(mur, incoming, outgoing)*tree_kin(
incoming, outgoing, check_momenta
);
}
} // namespace RHEJ
diff --git a/t/test_ME_h_3j.cc b/t/test_ME_h_3j.cc
index fe7bc17..a5fd45d 100644
--- a/t/test_ME_h_3j.cc
+++ b/t/test_ME_h_3j.cc
@@ -1,99 +1,104 @@
#include "LHEF/LHEF.h"
#include "RHEJ/config.hh"
#include "RHEJ/MatrixElement.hh"
#include "RHEJ/Event.hh"
struct Event{
std::array<RHEJ::Particle, 2> incoming;
std::vector<RHEJ::Particle> outgoing;
Event(
RHEJ::Particle in1, RHEJ::Particle in2,
std::initializer_list<RHEJ::Particle> out
):
incoming{std::move(in1), std::move(in2)},
outgoing{out}
{}
};
constexpr
RHEJ::pid::ParticleID PDG(int id){
return static_cast<RHEJ::pid::ParticleID>(id);
}
constexpr double alpha_s = 0.113559;
constexpr double mu = 125.;
constexpr double R = 0.4;
constexpr double min_jet_pt = 50.;
constexpr auto jet_def = fastjet::antikt_algorithm;
constexpr double ep = 1e-4;
void dump(Event const & ev){
LHEF::Writer writer{std::cerr};
std::cerr << std::setprecision(15);
RHEJ::UnclusteredEvent tmp;
tmp.incoming = ev.incoming;
tmp.outgoing = ev.outgoing;
tmp.central = {mu, mu, 0.};
RHEJ::Event out_ev{std::move(tmp), {jet_def, R}, min_jet_pt};
writer.hepeup = to_HEPEUP(std::move(out_ev), nullptr);
writer.writeEvent();
}
int main(){
std::vector<Event> events = {
/*
* reference matrix elements obtained from traditional HEJ (svn r3355)
* weights correspond to shat2*wt/(pdfPtr->pdfta*pdfPtr->pdftb)
* at the end of MEt2 in HJets.cxx
*/
#include "ME_test_events.dat"
};
- const std::vector<double> ref_weights{
+ std::vector<double> ref_weights{
1.20840366448943e-08,3.41482815561959e-18,4.09113550770465e-12,4.62704699562129e-15,4.48439126356004e-15,1.87972263060254e-12,1.44306589954128e-14,9.10660553817556e-13,6.78185181417658e-13,2.2313219279205e-14,5.04225136522914e-10,1.15445795400542e-13,4.25305992011228e-15,6.34464495221583e-11,8.1504882635896e-09,8.13050604628519e-10,1.16718768515638e-12,2.10641257506331e-08,6.62448031524987e-23,4.3319311109536e-07,1.20944494921349e-15,5.92539274973405e-05,4.4796169513938e-17,1.49273221119393e-15,2.43311569643869e-18,2.402863950606e-12,1.0620402696538e-11,5.35907540952987e-14,1.23031328706173e-12,2.79984692016006e-15,9.15681008462036e-19,6.82192590709917e-19,7.89237067022333e-13,1.86125787460036e-14,6.46979689455398e-07,3.04911830916725e-15,8.13163686285418e-15,1.53766124466997e-14,8.32276518653951e-14,8.88400099785712e-11,8.91102728936822e-06,1.80107644544759e-15,2.10582866365414e-09,3.10947796493188e-17,3.91243813845236e-15,3.26654480787734e-17,9.5447871679842e-14,2.59793669360488e-15,4.96134344054012e-13,5.81162371703576e-14,1.111167877034e-09,2.79109572058797e-16,4.46160513513727e-16,6.75218397576417e-15,2.68260561114305e-11,4.28989454505788e-16,5.8329868340234e-16,2.2024897389957e-09,2.57397955688743e-15,2.44654345522128e-14,2.30866752912176e-14,1.41827747056589e-07,1.08893147410797e-14,4.03382910318587e-11,3.1389869623674e-10,3.90505557430021e-13,1.08295237971538e-12,7.454591488958e-15,7.54483306433453e-14,1.02908259302562e-15,5.63531574394432e-14,6.30249429292698e-11,4.53762691808614e-13,1.89024041595359e-16,6.8158270448436e-15,1.14399844221906e-12,8.90256332512462e-11,6.06548092223796e-17,9.98362270779565e-12,4.59859716748231e-15,2.15099414771066e-10,2.11591109618363e-13,4.9742178019335e-09,1.33845681137089e-10,2.71186491627684e-16,8.86739836718398e-13,1.8249669781697e-14,1.92215499935624e-13,5.63538105112205e-14,4.14171568963499e-12,2.2717979957039e-15,5.55611156772101e-15,1.28364738716082e-16,5.3851134765245e-12,1.24211103260246e-08,4.83232530709888e-17,0.000521569751168655,1.6806219324588e-15,6.59977348355555e-11,9.1627970683281e-14,2.76820118228537e-12,4.0198982918336e-11,1.48426883756103e-18,1.40114232815435e-18,1.7965991343524e-11,7.6977628175923e-10,4.62196685302415e-14,2.69246069544232e-08,1.51516950568948e-13,6.99496684443045e-13,5.18573546306991e-16,0.0434630460883225,1.9837330600509e-23,2.67017601462925e-12,9.8878856049861e-16,8.72580364190897e-11,9.78094461016429e-06,8.71781003862465e-17,1.1559862791241e-08,4.84083438123941e-11,0.0139742563571754,1.99847601444777e-13,3.04915596353023e-17,5.93151845056274e-16,2.17435618150856e-14,5.33225221812805e-08,3.87773807827851e-14,4.50810202343348e-08,1.21607462589437e-08,3.1594598104258e-13,3.28656273761685e-16,2.50210924939855e-21,2.38126296013289e-13,1.25433376653286e-11,2.14322725557091e-14,1.49146107213853e-09,1.242585653879e-15,3.55954759815053e-17,3.58323190858796e-10,4.51037144634267e-14,7.55931651143187e-09,5.37698824199559e-14,7.38868781004529e-14,1.69357650230295e-11,5.80639629823124e-10,7.4837350747944e-17,6.96956083393151e-16,1.81834207579059e-10,2.82540028744915e-12,4.0395726561383e-13,1.79694539029015e-12,1.31041890677814e-16,1.30438189219159e-07,1.43221371664219e-10,9.114047729275e-16,7.34577355123687e-09,3.22593068239571e-15,4.15681922353836e-15,3.36408379142654e-16,1.44480854726148e-11,1.01363611458694e-13,1.00219358430676e-12,0.00338137229726986,4.29548017717404e-15,2.60395354258443e-14,4.1476296627427e-14,1.06855911096978e-13,8.22439039784712e-11,2.73629652512136e-10,5.15659349602634e-11,1.01584718933351e-13,6.94380398977803e-16,2.18931000347143e-18,3.07933496617785e-09,1.66390942837878e-16,3.52525907142114e-13,3.78080647432571e-14,1.17556440172277e-12,6.59751630435329e-17,1.21755902521401e-08,2.55450721249666e-11,1.53256550428156e-14,1.08380893991568e-12,4.16098986326189e-13,2.90205083341855e-14,5.511460594056e-18,1.57699398911427e-16,7.00104401760928e-13,3.58046798793545e-12,2.01229361268661e-12,3.73614518690553e-07,2.14710426845508e-14,3.70105600162546e-13,7.98071394035332e-15,7.37834963791502e-16,1.60529982053429e-09,1.08256973923745e-12,1.0334353499219e-13,9.49543099778791e-11,1.9115591364306e-12,4.0448657562838e-16,8.71422438179998e-15,1.25631778664749e-16,3.14479096095934e-19,2.64428535684163e-16,1.12107615349188e-11,1.23625298812219e-13,4.01428912903254e-13,2.70797308155832e-11,6.56279497820654e-09,7.13734059806893e-16,0.000782868291989126,3.76614026148548e-16,1.26257000484765e-13,9.58631950869197e-13,8.31708875806124e-11,2.34133301418612e-13,9.56863619261736e-14,4.64712032868663e-12,2.86185199150703e-12,1.33576465609689e-11,1.10079981208014e-15,1.0127254462106e-12,1.1462054990501e-12,6.6945137657857e-11,5.58956247731427e-17,9.55921087729154e-12,2.34217482968982e-13,9.06915368555449e-19,5.71005344841003e-15,7.30755001164554e-09,6.21879071656432e-07,3.18226338142145e-17,7.85187718829052e-15,1.12784019414768e-13,1.26620458266236e-12,3.69177699852588e-14,5.15075992359208e-08,4.92072565553936e-15,1.85781482577606e-18,4.45444190278706e-16,1.26266755594178e-13,2.4288203970545e-16,3.05176309462854e-13,1.10490541122231e-16,8.47749453725764e-16,5.15877353528789e-11,5.24254149207253e-17,1.91355494316452e-10,5.00556344079427e-17,2.17154626918891e-15,8.84358943960976e-15,3.61414423657153e-11,4.82735747825927e-09,1.29064865560832e-13,8.42127739585314e-10,3.92906216459766e-16,6.04446459041361e-17,5.89524177334148e-12,5.46194519536606e-19,1.34088171433487e-11,5.52860041827642e-15,1.75477788501097e-15,1.73667838135065e-16,1.65397026290053e-13,1.12185521855334e-09,6.08386662697057e-17,4.59141650497559e-07,1.5209325219688e-11,2.23037744763897e-13,2.01952877105987e-13,6.8597135743949e-17,1.21998302502671e-11,3.21171910404898e-17,5.26021980482437e-15,7.19705632175135e-19,1.20963419308456e-09,6.88528237531901e-13,1.98848279147693e-11,4.01262000581642e-12,5.20247739003617e-13,1.30428205622828e-14,4.11298291694491e-13,5.52043003142038e-10,1.23093415613032e-16,7.36616975809533e-17,3.45364902177034e-16,4.34985181424009e-08,5.38804276057675e-14,9.48300760869968e-13,1.22337104642822e-09,8.35787040119371e-14,1.55423397182344e-15,1.42203125672685e-13,3.18903277676011e-13,3.88652741986258e-15,7.04896809011702e-18,1.32646740227354e-21,7.99231856556843e-13,5.95802234870006e-09,1.65149306543534e-12,1.33064116256099e-11,1.56832702886549e-12,1.10299934631458e-08,3.67898989067657e-16,7.84682358054423e-17,2.57971628718788e-09,8.33182695406546e-17,1.18568992167241e-10,1.84355449893326e-12,2.55066898328171e-13,2.83673294248314e-14,1.78606190612042e-12,3.58348766614931e-12,4.59943073492537e-16,9.51522624162865e-13,8.0286195911246e-14,1.36695977212813e-16,2.66996900743536e-13,1.18874419459462e-15,4.20697505402443e-07,1.44157427471053e-12,3.88352207261609e-12,1.10352682706254e-15,9.16382540998255e-16,1.6918487768822e-14,6.27001020277801e-13,1.19982339734628e-15,2.08524585760556e-15,3.98635216491517e-11,2.29164410091261e-16,6.59993812570474e-14,9.09720299660866e-19,1.00791143935241e-13,3.48282488442434e-15,3.82462701684145e-11,1.82912920125976e-14,2.35934504263123e-12,3.07876467756239e-11,3.01958860729077e-05,1.34442938440306e-12,1.98172170171169e-11,2.54803366934403e-14,4.63525528427102e-13,1.98623090357032e-11,1.17532063246858e-12,1.62645374532443e-16,2.23093836031326e-11,5.23718977936426e-13,6.30254161980384e-15,5.64457070654369e-16,5.8598582297477e-15,1.54121478101166e-11,2.54757492800786e-11,7.24784320481052e-13,7.10590291356503e-11,1.90012029508776e-12,5.38421849395334e-05,9.2027530054952e-13,0.00802523961864568,1.69971085315301e-12,4.17128853468654e-15,2.58531014487319e-15,6.72871249369096e-12,3.68392430578061e-16,1.56192266889718e-14,8.09758960836169e-16,1.33670176619002e-07,3.40402526854793e-10,7.14029639915786e-15,1.15907463202726e-14
};
const std::vector<double> ref_virt{
0.000407437434887122,0.0616303680115081,0.0619338272462471,0.00113600815259533,0.0611805205486735,0.000565698082480416,0.516976611025519,0.0361140019339901,0.00145895360905874,0.134071980986426,0.000248975936203567,0.199294625870235,0.24358305091454,4.74361172738652e-05,5.89666196640687e-05,0.00531027816302915,0.000200243383490462,0.0805261711344379,0.00120027603203917,0.00201705178179895,0.339181495419129,0.000229863490455131,0.0523929260778993,0.000263057002724373,0.00244241108460007,0.0261467924657004,5.9657372131903e-05,0.200512468910312,0.00810099141722299,0.013973033333956,0.0108154079325531,0.158574821336736,0.00070998197888614,0.000498608284914822,3.73054997210496e-05,0.0589374030166269,0.00230622721548646,0.0339477696166811,0.000912463632412847,0.000736710687806361,0.00409015383244036,0.245951477936874,0.0837392661445892,0.222078993106222,0.00532590926455614,0.625630979237125,0.0496003388518398,0.10116468200925,0.0661786669164834,0.0210544045965245,0.0244625178289293,0.0146510511459362,0.101257018174142,0.00316763508016474,0.00143243489499069,0.333592621341247,0.00370907687545401,0.00213872056006003,0.01900030889465,0.167144065343861,2.6100961817368e-05,4.91550277446228e-05,0.0166499896061438,0.00011398434602661,0.000359768751250008,0.00225147558163185,0.0317147645964709,0.0842595514695342,0.35407612826878,0.00010526694382132,0.0430608610020923,3.5715013356257e-05,0.000838215519554185,0.00192018265054966,0.0713712284182796,0.000485232329561933,0.00152471141515176,0.000311862603200811,0.166619420201183,0.0199021050517686,1.41779057397457e-05,0.0017318628596011,0.00053760812759845,0.0686125677803084,0.0532620439395034,0.00166666380089786,0.0113533421860016,0.000109349820296255,0.0115667016173654,0.00173679164546385,0.0256607191277514,0.00120155627964714,0.0111815333349668,0.00349052851306921,0.00531735228893806,0.00152975285273546,4.66037248418355e-05,0.0425850611539657,0.00117697776465851,1.23387751381228e-05,0.0038237470220905,0.0905644410585002,0.02820127090728,0.00024367136851271,0.147636957219027,0.000593136877076053,0.0618471045067722,9.07687985864349e-06,0.000158960279301719,0.0887358461047651,0.0332353914977947,7.4341716367417e-06,0.000892008177580303,0.000741022113249675,0.00237741099974918,0.00354813689304541,0.00074173961384706,0.00467141577324516,0.00842296482521057,0.00821661392007147,4.83792120557372e-05,0.00425252516494632,0.000432204419963046,0.00586736518573752,0.000956530707046858,0.000236851602228008,0.0892427215701044,1.38830962250345e-05,0.0080267826500828,0.284615916615109,0.000228095910528196,0.000135342214864163,7.55341699407512e-06,0.000571884401074523,0.0101768016413147,0.0171951107578717,0.00189669876455143,0.00394930354845142,0.00328294862282636,0.00123193340562429,0.00331783483858544,0.00319051646342221,0.00155129659005506,8.87358953340508e-05,0.27623170363077,0.000223827470760779,0.0568015505265356,0.000286151801691539,0.0193226057191209,0.00123755083758374,0.000465248541258071,0.0760142560624272,0.00112490607994365,0.000144842352216183,0.197405681090406,0.00129883979140284,0.0930965339402953,0.0010308437485118,0.670013359364299,0.00432667098624779,0.0454915871138651,0.00251703872973342,0.00423441033466031,0.000345157767746845,0.00161415356368555,0.00416069958683314,0.00164939718639514,0.0111406419016027,0.00655421252542946,0.00851532238518754,0.00348709627275328,0.496059905302338,0.0371917367183897,0.000211531938236732,0.000258942974591912,0.17759337052515,0.00418622722335427,0.166120697713642,0.00360502919950296,0.000407414278464457,0.0360738279541873,0.0288893174873495,0.00047259015965458,1.56547248986137e-05,0.834956154090444,0.00732703650642937,0.00474657800944094,0.00273096554923287,0.257207154892208,0.00469808491585188,0.000270576144489275,0.000136075839471612,0.097205919234904,0.000607126031374257,0.000136619984636834,0.000160087925098193,0.0073858227644151,0.538982397271589,0.00336284303209108,0.00152656630103797,0.367113302445716,7.48038626150403e-05,2.60037367946894e-05,4.68016275421979e-05,0.000158257654442151,0.000139074963089943,0.000665698879408135,0.000343142682314476,0.0686327476602143,0.000242720484154181,0.0343708294206813,0.00114889968298977,0.0586454267607763,0.0681515891243611,0.00825429790468652,0.023992046337532,0.0176721744671367,0.000118812692493978,0.152442537631382,0.017153140472422,0.0064172004743378,3.15311961394686e-05,0.0176719977165195,0.00523267370956132,0.00696013339412051,0.0283129688131134,0.000387153677372906,0.00691199037138631,0.000146275451696295,0.701230156661611,0.000192846626056869,0.000164329434653947,0.00161262788918379,0.138649850457724,0.0298984433917698,0.0318122273046375,3.04827880239274e-05,0.000447390422746664,0.00113026494701542,0.00113285595566244,2.89098339981399e-05,0.00113856357833898,0.0539144571449436,0.0708152617112735,0.000146823107858066,0.0226524231026193,0.00337461891944858,0.000951319767163979,0.00173528759879105,0.027553837556441,0.000644790987241654,0.00202675597720372,0.0904812026317863,0.000121859990107361,0.0778263959001868,0.0078676495423855,0.00106041104696191,0.0460181415526085,0.111418229601818,0.00485578450515022,0.000301615731939099,0.0991074788208953,0.000545043001929007,0.00699096247420379,0.203553031240112,0.00129391506623192,0.0163036481873179,9.84784360850182e-06,0.0011930992876943,0.000233138513444116,0.00327134762224621,0.00934510369593952,0.0257690022271133,0.0586620958023085,0.0044269686967688,0.00458204111385776,0.000575173870997244,0.0196091953126155,0.0172446254918132,0.0023874923763825,0.00203444177442346,0.00590189106361814,0.160798676736245,0.000218465854877618,0.0449130791723244,0.0221907455403962,0.00577930006864582,4.68129848124974e-05,4.39088532565896e-05,0.0957291068496414,0.000222196152405895,0.00244492299475677,0.0649283175892222,0.00246575641284538,0.00395895850035338,0.452609135498808,2.18247701431579e-05,0.000103691495491759,0.0300254483476907,0.000997825085990438,0.00120168803041095,0.000488310402876861,5.02375648302248e-05,0.0889218155230037,0.886196353445691,0.537472442381394,0.000989158503397272,0.00585452117968649,0.000166646041974894,0.0559913557856292,0.271624614149917,1.9604193016588e-05,0.00400455993655262,0.000719953658022318,0.00769221521719457,0.0122222590022294,0.00869539856191144,0.00665334452394116,0.00269456106141176,0.0203144580365219,1.06269256412071e-05,0.00266305904353431,0.000866750257077464,0.0118678506617259,0.00774140571086026,0.472668086097713,0.00449314655439557,0.000363054262255041,0.000123958362353204,0.0593128481039698,6.49758497533422e-05,0.0147733538718899,0.000181422709422017,0.124358363237549,0.0149418567595413,0.000496790567694211,0.00319609174497007,0.171371757877329,0.000132314177013318,0.000118421517677613,0.0384411149022143,0.0560023977555578,0.000545951922843868,0.00923849055867683,7.42750022847054e-05,0.0034657141222543,0.0814951536826408,0.000247436489650107,0.00279135741337402,0.128217103654628,0.000239656955359697,0.0286081920788208,0.0264861644616376,0.000700223039958121,0.000224261722502375,0.206999869456071,0.0245477336638465,0.00020080171595025,0.000534370708560753,0.00114853261836608,0.0140588914482069,0.0491695416950957,0.0219392429988939,0.000187573794241001,2.76073824918907e-05,0.0177905351398601,0.0227102533101537,0.0156338659659345,0.0045130866637909,0.0118789100672146,0.0684121072806861
};
+
+ // reference weights contain an additional phase space factor
+ // remove it
+ for(auto & wt: ref_weights) wt *= 256.*pow(M_PI, 5);
+
RHEJ::MatrixElementConfig config;
config.jet_param.def = fastjet::JetDefinition(jet_def, R);
config.jet_param.min_pt = min_jet_pt;
config.log_correction = false;
config.Higgs_coupling = RHEJ::HiggsCouplingSettings{};
RHEJ::MatrixElement ME{
[](double){ return alpha_s; },
config
};
assert(events.size() == ref_weights.size());
for(size_t i = 0; i < events.size(); ++i){
std::sort(
begin(events[i].outgoing), end(events[i].outgoing),
RHEJ::rapidity_less{}
);
const double our_ME = ME.tree(
mu, events[i].incoming, events[i].outgoing, true
);
const double deviation = our_ME/ref_weights[i] - 1.;
if(std::abs(deviation) > ep){
std::cerr
<< "Matrix element deviates by factor of " << 1+deviation << ":\n"
"Is " << our_ME << ", should be " << ref_weights[i] << "\n"
"Event:\n";
dump(events[i]);
return EXIT_FAILURE;
}
const double our_virt = ME.virtual_corrections(
mu, events[i].incoming, events[i].outgoing
);
const double virt_deviation = our_virt/ref_virt[i] - 1.;
if(std::abs(virt_deviation) > ep){
std::cerr
<< "Virtual corrections deviate by factor of " << 1+virt_deviation << ":\n"
"Is " << our_virt << ", should be " << ref_virt[i] << "\n"
"Event:\n";
dump(events[i]);
return EXIT_FAILURE;
}
}
}
diff --git a/t/test_ME_hjets_mt174.cc b/t/test_ME_hjets_mt174.cc
index fe1cb0c..0016c51 100644
--- a/t/test_ME_hjets_mt174.cc
+++ b/t/test_ME_hjets_mt174.cc
@@ -1,90 +1,96 @@
#include "LHEF/LHEF.h"
#include "RHEJ/MatrixElement.hh"
#include "RHEJ/Event.hh"
struct Event{
double weight;
std::array<RHEJ::Particle, 2> incoming;
std::vector<RHEJ::Particle> outgoing;
Event(
double wt,
RHEJ::Particle in1, RHEJ::Particle in2,
std::initializer_list<RHEJ::Particle> out
):
weight{wt},
incoming{std::move(in1), std::move(in2)},
outgoing{out}
{}
};
constexpr
RHEJ::pid::ParticleID PDG(int id){
return static_cast<RHEJ::pid::ParticleID>(id);
}
constexpr double alpha_s = 0.112654;
constexpr double mu = 125.;
constexpr double R = 0.4;
constexpr double min_jet_pt = 35.;
constexpr auto jet_def = fastjet::antikt_algorithm;
constexpr double ep = 1e-4;
constexpr double alpha_s_mH_RHEJ = 0.113559;
void dump(Event const & ev){
LHEF::Writer writer{std::cerr};
std::cerr << std::setprecision(15);
RHEJ::UnclusteredEvent tmp;
tmp.incoming = ev.incoming;
tmp.outgoing = ev.outgoing;
tmp.central = {mu, mu, 0.};
RHEJ::Event out_ev{std::move(tmp), {jet_def, R}, min_jet_pt};
writer.hepeup = to_HEPEUP(std::move(out_ev), nullptr);
writer.writeEvent();
}
int main(){
std::vector<Event> events = {
/*
* reference matrix elements obtained from traditional HEJ (svn r3355)
* weights correspond to shat2*wt/(pdfPtr->pdfta*pdfPtr->pdftb)
* at the end of MEt2 in HJets.cxx
*/
#include "ME_hjets_mt174.dat"
};
+ // reference weights contain an additional phase space factor
+ // remove it
+ for(auto & event: events) {
+ event.weight *= 256.*pow(M_PI, 5);
+ }
+
RHEJ::MatrixElementConfig config;
config.jet_param.def = fastjet::JetDefinition(jet_def, R);
config.jet_param.min_pt = min_jet_pt;
config.log_correction = false;
RHEJ::HiggsCouplingSettings settings;
settings.mt = 174.;
settings.use_impact_factors = false;
settings.mb = 5.;
settings.include_bottom = true;
config.Higgs_coupling = settings;
RHEJ::MatrixElement ME{
[](double){ return alpha_s; },
config
};
for(size_t i = 0; i < events.size(); ++i){
std::sort(
begin(events[i].outgoing), end(events[i].outgoing),
RHEJ::rapidity_less{}
);
const double our_ME = ME.tree(
mu, events[i].incoming, events[i].outgoing, true
);
const double deviation = our_ME/events[i].weight - 1.;
if(std::abs(deviation) > ep){
std::cerr
<< "Matrix element deviates by factor of " << 1+deviation << ":\n"
"Is " << our_ME << ", should be " << events[i].weight << "\n"
"Event:\n";
dump(events[i]);
return EXIT_FAILURE;
}
}
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Jan 21, 1:44 AM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243484
Default Alt Text
(61 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment