Page MenuHomeHEPForge

No OneTemporary

This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
index 19cff80..b309881 100644
--- a/FixedOrderGen/src/PhaseSpacePoint.cc
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -1,710 +1,709 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "PhaseSpacePoint.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <iterator>
#include <limits>
#include <tuple>
#include <type_traits>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/kinematics.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PDF.hh"
#include "HEJ/RNG.hh"
#include "HEJ/utility.hh"
#include "JetParameters.hh"
#include "Process.hh"
#include "Subleading.hh"
namespace {
using namespace HEJ;
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
} // namespace anonymous
namespace HEJFOG{
Event::EventData to_EventData(PhaseSpacePoint psp){
//! @TODO Same function already in HEJ
Event::EventData result;
result.incoming = std::move(psp).incoming_;
assert(result.incoming.size() == 2);
result.outgoing = std::move(psp).outgoing_;
// technically Event::EventData doesn't have to be sorted,
// but PhaseSpacePoint should be anyway
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.decays = std::move(psp).decays_;
result.parameters.central = {NaN, NaN, psp.weight()};
return result;
}
namespace {
bool can_swap_to_uno(
Particle const & p1, Particle const & p2
){
return is_parton(p1)
&& p1.type != pid::gluon
&& p2.type == pid::gluon;
}
size_t count_gluons(std::vector<Particle>::const_iterator first,
std::vector<Particle>::const_iterator last
){
return std::count_if(first, last, [](Particle const & p)
{return p.type == pid::gluon;});
}
/** assumes FKL configurations between first and last,
* else there can be a quark in a non-extreme position
* e.g. uno configuration gqg would pass
*/
bool can_change_to_qqx(
std::vector<Particle>::const_iterator first,
std::vector<Particle>::const_iterator last
){
return 1 < count_gluons(first,last);
}
bool is_AWZ_proccess(Process const & proc){
return proc.boson && is_AWZ_boson(*proc.boson);
}
bool is_up_type(Particle const & part){
return is_anyquark(part) && !(std::abs(part.type)%2);
}
bool is_down_type(Particle const & part){
return is_anyquark(part) && std::abs(part.type)%2;
}
bool can_couple_to_W(Particle const & part, pid::ParticleID const W_id){
const int W_charge = W_id>0?1:-1;
return std::abs(part.type)<5
&& ( (W_charge*part.type > 0 && is_up_type(part))
|| (W_charge*part.type < 0 && is_down_type(part)) );
}
}
void PhaseSpacePoint::maybe_turn_to_subl(
double chance,
unsigned int const channels,
Process const & proc,
RNG & ran
){
if(proc.njets <= 2) return;
assert(outgoing_.size() >= 2);
// decide what kind of subleading process is allowed
bool allow_uno = false;
bool allow_strange = true;
const size_t nout = outgoing_.size();
const bool can_be_uno_backward = (channels&Subleading::uno)
&& can_swap_to_uno(outgoing_[0], outgoing_[1]);
const bool can_be_uno_forward = (channels&Subleading::uno)
&& can_swap_to_uno(outgoing_[nout-1], outgoing_[nout-2]);
allow_uno = can_be_uno_backward || can_be_uno_forward;
bool allow_qqx = (channels&Subleading::qqx)
&& can_change_to_qqx(outgoing_.cbegin(), outgoing_.cend());
if(is_AWZ_proccess(proc)) {
if(std::none_of(outgoing_.cbegin(), outgoing_.cend(),
[&proc](Particle const & p){
return can_couple_to_W(p, *proc.boson);})) {
// enforce qqx if A/W/Z can't couple somewhere else
assert(allow_qqx);
allow_uno = false;
chance = 1.;
// strange not allowed for W
if(std::abs(*proc.boson)== pid::Wp) allow_strange = false;
}
}
if(!allow_uno && !allow_qqx) return;
if(ran.flat() < chance){
weight_ /= chance;
if(allow_uno && !allow_qqx){
turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
} else if (!allow_uno && allow_qqx) {
turn_to_qqx(allow_strange, ran);
} else {
assert( allow_uno && allow_qqx);
if(ran.flat() < 0.5) turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
else turn_to_qqx(allow_strange, ran);
weight_ *= 2.;
}
} else weight_ /= 1 - chance;
}
void PhaseSpacePoint::turn_to_uno(
const bool can_be_uno_backward, const bool can_be_uno_forward,
RNG & ran
){
if(!can_be_uno_backward && !can_be_uno_forward) return;
const size_t nout = outgoing_.size();
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);
}
}
void PhaseSpacePoint::turn_to_qqx(const bool allow_strange, RNG & ran){
/// find first and last gluon in FKL chain
auto first = std::find_if(outgoing_.begin(), outgoing_.end(),
[](Particle const & p){return p.type == pid::gluon;});
std::vector<Particle*> FKL_gluons;
for(auto p = first; p!=outgoing_.end(); ++p){
if(p->type == pid::gluon) FKL_gluons.push_back(&*p);
else if(is_anyquark(*p)) break;
}
const size_t ng = FKL_gluons.size();
if(ng < 2)
throw std::logic_error("not enough gluons to create qqx");
// select flavour of quark
const double r1 = 2.*ran.flat()-1.;
const double max_flavour = allow_strange?n_f:n_f-1;
weight_ *= max_flavour*2;
int flavour = pid::down + std::floor(std::abs(r1)*max_flavour);
flavour*=r1<0.?-1:1;
// select gluon for switch
const size_t idx = std::floor((ng-1) * ran.flat());
weight_ *= (ng-1);
FKL_gluons[idx]->type = ParticleID(flavour);
FKL_gluons[idx+1]->type = ParticleID(-flavour);
}
template<class ParticleMomenta>
fastjet::PseudoJet PhaseSpacePoint::gen_last_momentum(
ParticleMomenta const & other_momenta,
const double mass_square, const double y
) const {
std::array<double,2> pt{0.,0.};
for (auto const & p: other_momenta) {
pt[0]-= p.px();
pt[1]-= p.py();
}
const double mperp = std::sqrt(pt[0]*pt[0]+pt[1]*pt[1]+mass_square);
const double pz=mperp*std::sinh(y);
const double E=mperp*std::cosh(y);
return {pt[0], pt[1], pz, E};
}
namespace {
//! adds a particle to target (in correct rapidity ordering)
//! @returns positon of insertion
auto insert_particle(std::vector<Particle> & target,
Particle && particle
){
const auto pos = std::upper_bound(
begin(target),end(target),particle,rapidity_less{}
);
target.insert(pos, std::move(particle));
return pos;
}
}
PhaseSpacePoint::PhaseSpacePoint(
Process const & proc,
JetParameters const & jet_param,
PDF & pdf, double E_beam,
double const subl_chance,
unsigned int const subl_channels,
ParticlesDecayMap const & particle_decays,
EWConstants const & ew_parameters,
RNG & ran
- )
- {
+ ){
assert(proc.njets >= 2);
status_ = good;
weight_ = 1;
const int nout = proc.njets + (proc.boson?1:0) + proc.boson_decay.size();
outgoing_.reserve(nout);
// generate parton momenta
const bool is_pure_jets = (nout == proc.njets);
auto partons = gen_LO_partons(
proc.njets, is_pure_jets, jet_param, E_beam, ran
);
// pre fill flavour with gluons
for(auto&& p_out: partons) {
outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out), {}});
}
if(status_ != good) return;
if(proc.boson){ // decay boson
const auto & boson_prop = ew_parameters.prop(*proc.boson) ;
auto boson{ gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran) };
const auto pos{insert_particle(outgoing_, std::move(boson))};
const size_t boson_idx = std::distance(begin(outgoing_), pos);
const auto & boson_decay = particle_decays.find(*proc.boson);
if( !proc.boson_decay.empty() ){ // decay given in proc
decays_.emplace(
boson_idx,
decay_boson(outgoing_[boson_idx], proc.boson_decay, ran)
);
} else if( boson_decay != particle_decays.end()
&& !boson_decay->second.empty() ){ // decay given explicitly
decays_.emplace(
boson_idx,
decay_boson(outgoing_[boson_idx], boson_decay->second, ran)
);
}
}
// normalisation of momentum-conserving delta function
weight_ *= std::pow(2*M_PI, 4);
/** @TODO
* uf (jet_param.min_pt) doesn't correspond to our final scale choice.
* The HEJ scale generators currently expect a full event as input,
* so fixing this is not completely trivial
*/
reconstruct_incoming(proc, subl_channels, 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;
maybe_turn_to_subl(subl_chance, subl_channels, proc, ran);
if(proc.boson) couple_boson(*proc.boson, ran);
}
double PhaseSpacePoint::gen_hard_pt(
int np , double ptmin, double ptmax, double y,
RNG & ran
) {
// heuristic parameters for pt sampling
const double ptpar = ptmin + np/5.;
const double arg_small_y = std::atan((ptmax - ptmin)/ptpar);
const double y_cut = 3.;
const double r1 = ran.flat();
if(y < y_cut){
const double pt = ptmin + ptpar*std::tan(r1*arg_small_y);
const double temp = std::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, 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,
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, bool is_pure_jets,
JetParameters const & jet_param,
double max_pt,
RNG & ran
){
if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"};
weight_ /= std::pow(16.*std::pow(M_PI,3),np);
weight_ /= std::tgamma(np+1); //remove rapidity ordering
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 bool is_last_parton = i+1 == np;
if(is_pure_jets && is_last_parton) {
constexpr double parton_mass_sq = 0.;
partons.emplace_back(gen_last_momentum(partons, parton_mass_sq, y));
break;
}
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(
ParticleID bosonid, double mass, double width,
RNG & ran
){
// Usual phase space measure
weight_ /= 16.*std::pow(M_PI, 3);
// Generate a y Gaussian distributed around 0
/// @TODO check magic numbers for different boson Higgs
/// @TODO better sampling for W
const double stddev_y = 1.6;
const double y = random_normal(stddev_y, ran);
const double r1 = ran.flat();
const double s_boson = mass*(
mass + width*std::tan(M_PI/2.*r1 + (r1-1.)*std::atan(mass/width))
);
// off-shell s_boson sampling, compensates for Breit-Wigner
/// @TODO use a flag instead
if(std::abs(bosonid) == pid::Wp){
weight_/=M_PI*M_PI*8.;
weight_*= mass*width*( M_PI+2.*std::atan(mass/width) )
/ ( 1. + std::cos( M_PI*r1 + 2.*(r1-1.)*std::atan(mass/width) ) );
}
auto p = gen_last_momentum(outgoing_, s_boson, y);
return Particle{bosonid, std::move(p), {}};
}
Particle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> const & partons
) const{
if(!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(!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(!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(!is_parton(partons[last_idx])) return partons[last_idx-1];
return partons[last_idx];
}
namespace {
/// partons are ordered: even = anti, 0 = gluon
ParticleID index_to_pid(size_t i){
if(!i) return pid::gluon;
return static_cast<ParticleID>( i%2 ? (i+1)/2 : -i/2 );
}
/// partons are ordered: even = anti, 0 = gluon
size_t pid_to_index(ParticleID id){
if(id==pid::gluon) return 0;
return id>0 ? id*2-1 : std::abs(id)*2;
}
std::bitset<11> init_allowed(ParticleID const id){
if(std::abs(id) == pid::proton)
return ~0;
std::bitset<11> out{0};
if(is_parton(id))
out[pid_to_index(id)] = 1;
return out;
}
/// decides which "index" (see index_to_pid) are allowed for process
std::bitset<11> allowed_quarks(ParticleID const boson){
std::bitset<11> allowed = ~0;
if(std::abs(boson) == pid::Wp){
// special case W:
// Wp: anti-down or up-type quark, no b/t
// Wm: down or anti-up-type quark, no b/t
allowed = boson>0? 0b00011001101
:0b00100110011;
}
return allowed;
}
}
/**
* checks which partons are allowed as initial state:
* 1. only allow what is given in the Runcard (p -> all)
* 2. A/W/Z require something to couple to
* a) no qqx => no incoming gluon
* b) 2j => no incoming gluon
* c) 3j => can couple OR is gluon => 2 gluons become qqx later
*/
std::array<std::bitset<11>,2> PhaseSpacePoint::filter_partons(
Process const & proc, unsigned int const subl_channels, RNG & ran
){
std::array<std::bitset<11>,2> allowed_partons{
init_allowed(proc.incoming[0]),
init_allowed(proc.incoming[1])
};
bool const allow_qqx = subl_channels&Subleading::qqx;
// special case A/W/Z
if(is_AWZ_proccess(proc) && ((proc.njets < 4) || !allow_qqx)){
// all possible incoming states
auto allowed(allowed_quarks(*proc.boson));
if(proc.njets == 2 || !allow_qqx) allowed[0]=0;
// possible states per leg
std::array<std::bitset<11>,2> const maybe_partons{
allowed_partons[0]&allowed, allowed_partons[1]&allowed};
if(maybe_partons[0].any() && maybe_partons[1].any()){
// two options to get allowed initial state => choose one at random
const size_t idx = ran.flat() < 0.5;
allowed_partons[idx] = maybe_partons[idx];
// else choose the possible
} else if(maybe_partons[0].any()) {
allowed_partons[0] = maybe_partons[0];
} else if(maybe_partons[1].any()) {
allowed_partons[1] = maybe_partons[1];
} else{
throw std::invalid_argument{"Incoming state not allowed."};
}
}
return allowed_partons;
}
void PhaseSpacePoint::reconstruct_incoming(
Process const & proc, unsigned int const subl_channels,
PDF & pdf, double E_beam,
double uf,
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].E()-incoming_[0].pz())/sqrts;
const double xb=(incoming_[1].E()+incoming_[1].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;
}
auto const & ids = proc.incoming;
std::array<std::bitset<11>,2> allowed_partons(
filter_partons(proc, subl_channels, ran));
for(size_t i = 0; i < 2; ++i){
if(ids[i] == pid::proton || ids[i] == pid::p_bar){
// pick ids according to pdfs
incoming_[i].type =
generate_incoming_id(i, i?xb:xa, uf, pdf, allowed_partons[i], ran);
} else {
assert(allowed_partons[i][pid_to_index(ids[i])]);
incoming_[i].type = ids[i];
}
}
assert(momentum_conserved(1e-7));
}
ParticleID PhaseSpacePoint::generate_incoming_id(
size_t const beam_idx, double const x, double const uf,
PDF & pdf, std::bitset<11> allowed_partons, RNG & ran
){
std::array<double,11> pdf_wt;
pdf_wt[0] = allowed_partons[0]?
std::fabs(pdf.pdfpt(beam_idx,x,uf,pid::gluon)):0.;
double pdftot = pdf_wt[0];
for(size_t i = 1; i < pdf_wt.size(); ++i){
pdf_wt[i] = allowed_partons[i]?
4./9.*std::fabs(pdf.pdfpt(beam_idx,x,uf,index_to_pid(i))):0;
pdftot += pdf_wt[i];
}
const double r1 = pdftot * ran.flat();
double sum = 0;
for(size_t i=0; i < pdf_wt.size(); ++i){
if (r1 < (sum+=pdf_wt[i])){
weight_*= pdftot/pdf_wt[i];
return index_to_pid(i);
}
}
std::cerr << "Error in choosing incoming parton: "<<x<<" "<<uf<<" "
<<sum<<" "<<pdftot<<" "<<r1<<std::endl;
throw std::logic_error{"Failed to choose parton flavour"};
}
void PhaseSpacePoint::couple_boson(
ParticleID const boson, RNG & ran
){
if(std::abs(boson) != pid::Wp) return; // only matters for W
/// @TODO this could be use to sanity check gamma and Z
// find all possible quarks
std::vector<Particle*> allowed_parts;
for(auto & part: outgoing_){
// Wp -> up OR anti-down, Wm -> anti-up OR down, no bottom
if ( can_couple_to_W(part, boson) )
allowed_parts.push_back(&part);
}
if(allowed_parts.size() == 0){
throw std::logic_error{"Found no parton for coupling with boson"};
}
// select one and flip it
size_t idx = 0;
if(allowed_parts.size() > 1){
/// @TODO more efficient sampling
/// old code: probability[i] = exp(parton[i].y - W.y)
idx = std::floor(ran.flat()*allowed_parts.size());
weight_ *= allowed_parts.size();
}
const int W_charge = boson>0?1:-1;
allowed_parts[idx]->type =
static_cast<ParticleID>( allowed_parts[idx]->type - W_charge );
}
double PhaseSpacePoint::random_normal(
double stddev,
RNG & ran
){
const double r1 = ran.flat();
const double r2 = ran.flat();
const double lninvr1 = -std::log(r1);
const double result = stddev*std::sqrt(2.*lninvr1)*std::cos(2.*M_PI*r2);
weight_ *= exp(result*result/(2*stddev*stddev))*std::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,
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;
if(decays.size()==1) return decays.front();
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(
Particle const & parent,
std::vector<Decay> const & decays,
RNG & ran
){
const auto channel = select_decay_channel(decays, ran);
return decay_boson(parent, channel.products, ran);
}
std::vector<Particle> PhaseSpacePoint::decay_boson(
Particle const & parent,
std::vector<ParticleID> const & decays,
RNG & ran
){
if(decays.size() != 2){
throw not_implemented{
"only decays into two particles are implemented"
};
}
std::vector<Particle> decay_products(decays.size());
for(size_t i = 0; i < decays.size(); ++i){
decay_products[i].type = decays[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 = std::sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
const double px = E*std::cos(theta)*sin_phi;
const double py = E*std::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/include/HEJ/BufferedEventReader.hh b/include/HEJ/BufferedEventReader.hh
index 77a6231..3b70a0f 100644
--- a/include/HEJ/BufferedEventReader.hh
+++ b/include/HEJ/BufferedEventReader.hh
@@ -1,69 +1,69 @@
/** \file
* \brief Event reader with internal buffer
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <stack>
#include <cstddef>
#include <string>
#include <utility>
#include <vector>
#include "LHEF/LHEF.h"
#include "HEJ/EventReader.hh"
#include "HEJ/optional.hh"
namespace HEJ {
//! Event reader with internal buffer
/**
* Having a buffer makes it possible to put events back into the reader
* so that they are read again
*/
class BufferedEventReader: public EventReader {
public:
explicit BufferedEventReader(std::unique_ptr<EventReader> reader):
reader_{std::move(reader)}
{}
//! Read an event
bool read_event() override;
//! Access header text
std::string const & header() const override {
return reader_->header();
}
//! Access run information
LHEF::HEPRUP const & heprup() const override {
return reader_->heprup();
}
//! Access last read event
LHEF::HEPEUP const & hepeup() const override {
return cur_event_;
}
//! Guess number of events from header
- HEJ::optional<std::size_t> number_events() const override {
+ optional<std::size_t> number_events() const override {
return reader_->number_events();
}
//! Push event back into reader
template< class... T>
void emplace(T&&... args) {
buffer_.emplace(std::forward<T>(args)...);
}
private:
std::stack<LHEF::HEPEUP, std::vector<LHEF::HEPEUP>> buffer_;
std::unique_ptr<EventReader> reader_;
LHEF::HEPEUP cur_event_;
};
} // namespace HEJ
diff --git a/include/HEJ/Config.hh b/include/HEJ/Config.hh
index 70d64de..99b8224 100644
--- a/include/HEJ/Config.hh
+++ b/include/HEJ/Config.hh
@@ -1,234 +1,234 @@
/** \file
* \brief HEJ 2 configuration parameters
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <map>
#include <string>
#include <utility>
#include <vector>
#include "fastjet/JetDefinition.hh"
#include "yaml-cpp/yaml.h"
#include "HEJ/Constants.hh"
#include "HEJ/event_types.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/Fraction.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/optional.hh"
#include "HEJ/output_formats.hh"
#include "HEJ/ScaleFunction.hh"
namespace HEJ{
//! Jet parameters
struct JetParameters{
fastjet::JetDefinition def; /**< Jet Definition */
double min_pt; /**< Minimum Jet Transverse Momentum */
};
//! Settings for scale variation
struct ScaleConfig{
//! Base scale choices
std::vector<ScaleFunction> base;
//! Factors for multiplicative scale variation
std::vector<double> factors;
//! Maximum ratio between renormalisation and factorisation scale
double max_ratio;
};
//! Settings for random number generator
struct RNGConfig {
//! Random number generator name
std::string name;
//! Optional initial seed
optional<std::string> seed;
};
//! Settings for partial unweighting
struct PartialUnweightConfig {
//! Number of trials for training
size_t trials;
//! Maximum distance in standard deviations from mean logarithmic weight
double max_dev;
};
/**! Possible treatments for fixed-order input events.
*
* The program will decide on how to treat an event based on
* the value of this enumeration.
*/
enum class EventTreatment{
reweight, /**< Perform resummation */
keep, /**< Keep the event */
discard, /**< Discard the event */
};
//! Container to store the treatments for various event types
using EventTreatMap = std::map<event_type::EventType, EventTreatment>;
//! Possible setting for the event weight
enum class WeightType{
weighted, //!< weighted events
unweighted_resum, //!< unweighted only resummation part
partially_unweighted //!< mixed weighted and unweighted
};
/**! Input parameters.
*
* This struct handles stores all configuration parameters
* needed in a HEJ 2 run.
*
* \internal To add a new option:
* 1. Add a member to the Config struct.
* 2. Inside "src/YAMLreader.cc":
* - Add the option name to the "supported" Node in
* get_supported_options.
* - Initialise the new Config member in to_Config.
* The functions set_from_yaml (for mandatory options) and
* set_from_yaml_if_defined (non-mandatory) may be helpful.
* 3. Add a new entry (with short description) to config.yaml
* 4. Update the user documentation in "doc/Sphinx/"
*/
struct Config {
//! %Parameters for scale variation
ScaleConfig scales;
//! Resummation jet properties
JetParameters resummation_jets;
//! Fixed-order jet properties
JetParameters fixed_order_jets;
//! Minimum transverse momentum for extremal partons
//! \deprecated This will be removed in future versions.
//! Use \ref max_ext_soft_pt_fraction instead.
double min_extparton_pt = 0;
//! Maximum transverse momentum fraction from soft radiation in extremal jets
Fraction<double> max_ext_soft_pt_fraction;
//! The regulator lambda for the subtraction terms
double regulator_lambda = CLAMBDA;
//! Number of resummation configurations to generate per fixed-order event
size_t trials;
//! Maximal number of events
optional<size_t> max_events;
//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
bool log_correction;
//! Event output files names and formats
std::vector<OutputFile> output;
//! Parameters for random number generation
RNGConfig rng;
//! Map to decide what to do for different event types
EventTreatMap treat;
//! %Parameters for custom analysis
//! @deprecated use analyses_parameters instead
YAML::Node analysis_parameters;
//! %Parameters for custom analyses
std::vector<YAML::Node> analyses_parameters;
//! Settings for effective Higgs-gluon coupling
HiggsCouplingSettings Higgs_coupling;
//! elector weak parameters
EWConstants ew_parameters;
//! Type of event weight e.g. (un)weighted
WeightType weight_type;
//! Settings for partial unweighting
- HEJ::optional<PartialUnweightConfig> unweight_config;
+ optional<PartialUnweightConfig> unweight_config;
};
//! Configuration options for the PhaseSpacePoint class
struct PhaseSpacePointConfig {
//! Properties of resummation jets
JetParameters jet_param;
//! Minimum transverse momentum for extremal partons
//! \deprecated This will be removed in future versions.
//! Use \ref max_ext_soft_pt_fraction instead.
double min_extparton_pt = 0;
//! Maximum transverse momentum fraction from soft radiation in extremal jets
Fraction<double> max_ext_soft_pt_fraction;
};
//! Configuration options for the MatrixElement class
struct MatrixElementConfig {
MatrixElementConfig() = default;
MatrixElementConfig(
bool log_correction,
HiggsCouplingSettings Higgs_coupling,
EWConstants ew_parameters,
double regulator_lambda = CLAMBDA
):
log_correction{log_correction},
Higgs_coupling{std::move(Higgs_coupling)},
ew_parameters{std::move(ew_parameters)},
regulator_lambda{regulator_lambda}
{}
//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
bool log_correction;
//! Settings for effective Higgs-gluon coupling
HiggsCouplingSettings Higgs_coupling;
//! elector weak parameters
EWConstants ew_parameters;
//! The regulator lambda for the subtraction terms
double regulator_lambda = CLAMBDA;
};
//! Configuration options for the EventReweighter class
struct EventReweighterConfig {
//! Settings for phase space point generation
PhaseSpacePointConfig psp_config;
//! Settings for matrix element calculation
MatrixElementConfig ME_config;
//! Access properties of resummation jets
JetParameters & jet_param() {
return psp_config.jet_param;}
//! Access properties of resummation jets (const version)
JetParameters const & jet_param() const {
return psp_config.jet_param;}
//! Treatment of the various event types
EventTreatMap treat;
};
/**! Extract PhaseSpacePointConfig from Config
*
* \internal We do not provide a PhaseSpacePointConfig constructor from Config
* so that PhaseSpacePointConfig remains an aggregate.
* This faciliates writing client code (e.g. the HEJ fixed-order generator)
* that creates a PhaseSpacePointConfig *without* a Config object.
*
* @see to_MatrixElementConfig, to_EventReweighterConfig
*/
inline
PhaseSpacePointConfig to_PhaseSpacePointConfig(Config const & conf) {
return {
conf.resummation_jets,
conf.min_extparton_pt,
conf.max_ext_soft_pt_fraction
};
}
/**! Extract MatrixElementConfig from Config
*
* @see to_PhaseSpacePointConfig, to_EventReweighterConfig
*/
inline
MatrixElementConfig to_MatrixElementConfig(Config const & conf) {
return {conf.log_correction, conf.Higgs_coupling,
conf.ew_parameters, conf.regulator_lambda};
}
/**! Extract EventReweighterConfig from Config
*
* @see to_PhaseSpacePointConfig, to_MatrixElementConfig
*/
inline
EventReweighterConfig to_EventReweighterConfig(Config const & conf) {
return {
to_PhaseSpacePointConfig(conf),
to_MatrixElementConfig(conf),
conf.treat
};
}
} // namespace HEJ
diff --git a/include/HEJ/CrossSectionAccumulator.hh b/include/HEJ/CrossSectionAccumulator.hh
index c5ca18f..a846c06 100644
--- a/include/HEJ/CrossSectionAccumulator.hh
+++ b/include/HEJ/CrossSectionAccumulator.hh
@@ -1,99 +1,99 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <iterator>
#include <map>
#include <ostream>
#include <vector>
#include "HEJ/event_types.hh"
namespace HEJ {
class Event;
//! Collection of Cross Section with its uncertainty
template<typename T>
struct XSWithError {
T value = T{}; //!< Cross Section
T error = T{}; //!< Error
};
/**
* @brief Sum of Cross Section for different subproccess
*/
class CrossSectionAccumulator {
public:
//! Fill with single event
//! @note for multiple resummation events use fill_correlated() instead
void fill(Event const & ev);
//! Fill by weight and type
//! @note for multiple resummation events use fill_correlated() instead
void fill(double weight, event_type::EventType type);
//! Fill by weight, error and type
//! @note The error will be _added_ to the current error
void fill(double weight, double delta_error, event_type::EventType type);
/**
* @brief Fill with multiple correlated weights
* @details This should be used to fill multiple reweighted events,
* coming from the same fixed order point.
* Total error for \f$N\f$ fixed order points each giving \f$M_i\f$
* resummation events is:
* \f[
* \delta^2=\sum_i \left(\sum_j w_{i,j}\right)^2
* +\sum_{i,j} \left(w_{i,j}\right)^2,
* \f]
* @note This is equivalent to fill() for only one reweighted event
* coming from each fixed order point (\f$M_i=1\f$)
*/
void fill_correlated(std::vector<Event> const & evts);
//! iterator implementation of fill_correlated()
template<class ConstIt>
void fill_correlated(ConstIt begin, ConstIt end);
//! explicit version of fill_correlated() by giving sum(wt) and sum(wt^2)
void fill_correlated(double sum, double sum2, event_type::EventType type);
//! begin of Cross Section and error for subprocesses
auto begin() const {
return std::begin(xs_);
}
//! end of Cross Section and error for subprocesses
auto end() const {
return std::end(xs_);
}
//! total Cross Section and error
XSWithError<double> const & total() const {
return total_;
}
private:
- std::map<HEJ::event_type::EventType, XSWithError<double>> xs_;
+ std::map<event_type::EventType, XSWithError<double>> xs_;
XSWithError<double> total_;
};
//! Print CrossSectionAccumulator to stream
std::ostream& operator<<(std::ostream& os, const CrossSectionAccumulator& xs);
// ------------ Implementation ------------
template<class ConstIt>
void CrossSectionAccumulator::fill_correlated(ConstIt begin, ConstIt end){
if(std::distance(begin, end) < 2){ // only one event
fill(*begin);
return;
}
double sum = 0.;
double sum2 = 0.;
const auto type = begin->type();
for(; begin != end; ++begin){
double const wt = begin->central().weight;
sum += wt;
sum2 += wt*wt;
}
fill_correlated(sum, sum2, type);
}
} // namespace HEJ
diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
index 730b622..6b01b9e 100644
--- a/include/HEJ/Event.hh
+++ b/include/HEJ/Event.hh
@@ -1,372 +1,372 @@
/** \file
* \brief Declares the Event class and helpers
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <cstddef>
#include <iosfwd>
#include <iterator>
#include <unordered_map>
#include <utility>
#include <vector>
#include "boost/iterator/filter_iterator.hpp"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/PseudoJet.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Parameters.hh"
#include "HEJ/Particle.hh"
namespace LHEF {
class HEPEUP;
class HEPRUP;
}
namespace fastjet {
class JetDefinition;
}
namespace HEJ {
struct RNG;
struct UnclusteredEvent;
/** @brief An event with clustered jets
*
* This is the main HEJ 2 event class.
* It contains kinematic information including jet clustering,
* parameter (e.g. scale) settings and the event weight.
*/
class Event {
public:
class EventData;
//! Iterator over partons
using ConstPartonIterator = boost::filter_iterator<
bool (*)(Particle const &),
std::vector<Particle>::const_iterator
>;
//! Reverse Iterator over partons
using ConstReversePartonIterator = std::reverse_iterator<
ConstPartonIterator>;
//! No default Constructor
Event() = delete;
//! Event Constructor adding jet clustering to an unclustered event
//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0
[[deprecated("UnclusteredEvent will be replaced by EventData")]]
Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
);
//! @name Particle Access
//! @{
//! Incoming particles
std::array<Particle, 2> const & incoming() const{
return incoming_;
}
//! Outgoing particles
std::vector<Particle> const & outgoing() const{
return outgoing_;
}
//! Iterator to the first outgoing parton
ConstPartonIterator begin_partons() const;
//! Iterator to the first outgoing parton
ConstPartonIterator cbegin_partons() const;
//! Iterator to the end of the outgoing partons
ConstPartonIterator end_partons() const;
//! Iterator to the end of the outgoing partons
ConstPartonIterator cend_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator rbegin_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator crbegin_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator rend_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator crend_partons() const;
//! Particle decays
/**
* The key in the returned map corresponds to the index in the
* vector returned by outgoing()
*/
std::unordered_map<std::size_t, std::vector<Particle>> const & decays()
const {
return decays_;
}
//! The jets formed by the outgoing partons, sorted in rapidity
std::vector<fastjet::PseudoJet> const & jets() const{
return jets_;
}
//! @}
//! @name Weight variations
//! @{
//! All chosen parameter, i.e. scale choices (const version)
Parameters<EventParameters> const & parameters() const{
return parameters_;
}
//! All chosen parameter, i.e. scale choices
Parameters<EventParameters> & parameters(){
return parameters_;
}
//! Central parameter choice (const version)
EventParameters const & central() const{
return parameters_.central;
}
//! Central parameter choice
EventParameters & central(){
return parameters_.central;
}
//! Parameter (scale) variations (const version)
std::vector<EventParameters> const & variations() const{
return parameters_.variations;
}
//! Parameter (scale) variations
std::vector<EventParameters> & variations(){
return parameters_.variations;
}
//! Parameter (scale) variation (const version)
/**
* @param i Index of the requested variation
*/
EventParameters const & variations(std::size_t i) const{
return parameters_.variations.at(i);
}
//! Parameter (scale) variation
/**
* @param i Index of the requested variation
*/
EventParameters & variations(std::size_t i){
return parameters_.variations.at(i);
}
//! @}
//! Indices of the jets the outgoing partons belong to
/**
* @param jets Jets to be tested
* @returns A vector containing, for each outgoing parton,
* the index in the vector of jets the considered parton
* belongs to. If the parton is not inside any of the
* passed jets, the corresponding index is set to -1.
*/
std::vector<int> particle_jet_indices(
std::vector<fastjet::PseudoJet> const & jets
) const {
return cs_.particle_jet_indices(jets);
}
//! particle_jet_indices() of the Event jets()
std::vector<int> particle_jet_indices() const {
return particle_jet_indices(jets());
}
//! Jet definition used for clustering
fastjet::JetDefinition const & jet_def() const{
return cs_.jet_def();
}
//! Minimum jet transverse momentum
double min_jet_pt() const{
return min_jet_pt_;
}
//! Event type
event_type::EventType type() const{
return type_;
}
//! Give colours to each particle
/**
* @returns true if new colours are generated, i.e. same as is_resummable()
* @details Colour ordering is done according to leading colour in the MRK
* limit, see \cite Andersen:2011zd. This only affects \ref
* is_resummable() "HEJ" configurations, all other \ref event_type
* "EventTypes" will be ignored.
* @note This overwrites all previously set colours.
*/
- bool generate_colours(HEJ::RNG &);
+ bool generate_colours(RNG &);
//! Check that current colours are leading in the high energy limit
/**
* @details Checks that the colour configuration can be split up in
* multiple, rapidity ordered, non-overlapping ladders. Such
* configurations are leading in the MRK limit, see
* \cite Andersen:2011zd
*
* @note This is _not_ to be confused with \ref is_resummable(), however
* for all resummable states it is possible to create a leading colour
* configuration, see generate_colours()
*/
bool is_leading_colour() const;
/**
* @brief Check if given event could have been produced by HEJ
* @details A HEJ state has to fulfil:
* 1. type() has to be \ref is_resummable() "resummable"
* 2. Soft radiation in the tagging jets contributes at most to
* `max_ext_soft_pt_fraction` of the total jet \f$ p_\perp \f$
*
* @note This is true for any resummed stated produced by the
* EventReweighter or any \ref is_resummable() "resummable" Leading
* Order state.
*
* @param max_ext_soft_pt_fraction Maximum transverse momentum fraction from
* soft radiation in extremal jets
* @param min_extparton_pt Absolute minimal \f$ p_\perp \f$,
* \b deprecated use max_ext_soft_pt_fraction
* instead
* @return True if this state could have been produced by HEJ
*/
bool valid_hej_state(
double max_ext_soft_pt_fraction, double min_extparton_pt = 0.) const;
private:
//! \internal
//! @brief Construct Event explicitly from input.
/** This is only intended to be called from EventData.
*
* \warning The input is taken _as is_, sorting and classification has to be
* done externally, i.e. by EventData
*/
Event(
std::array<Particle, 2> && incoming,
std::vector<Particle> && outgoing,
std::unordered_map<std::size_t, std::vector<Particle>> && decays,
Parameters<EventParameters> && parameters,
fastjet::JetDefinition const & jet_def,
double const min_jet_pt
);
//! Iterator over partons (non-const)
using PartonIterator = boost::filter_iterator<
bool (*)(Particle const &),
std::vector<Particle>::iterator
>;
//! Reverse Iterator over partons (non-const)
using ReversePartonIterator = std::reverse_iterator<PartonIterator>;
//! Iterator to the first outgoing parton (non-const)
PartonIterator begin_partons();
//! Iterator to the end of the outgoing partons (non-const)
PartonIterator end_partons();
//! Reverse Iterator to the first outgoing parton (non-const)
ReversePartonIterator rbegin_partons();
//! Reverse Iterator to the first outgoing parton (non-const)
ReversePartonIterator rend_partons();
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
std::unordered_map<std::size_t, std::vector<Particle>> decays_;
std::vector<fastjet::PseudoJet> jets_;
Parameters<EventParameters> parameters_;
fastjet::ClusterSequence cs_;
double min_jet_pt_;
event_type::EventType type_;
}; // end class Event
//! Class to store general Event setup, i.e. Phase space and weights
class Event::EventData {
public:
//! Default Constructor
EventData() = default;
//! Constructor from LesHouches event information
EventData(LHEF::HEPEUP const & hepeup);
//! Constructor with all values given
EventData(
std::array<Particle, 2> incoming,
std::vector<Particle> outgoing,
std::unordered_map<std::size_t, std::vector<Particle>> decays,
Parameters<EventParameters> parameters
):
incoming(std::move(incoming)), outgoing(std::move(outgoing)),
decays(std::move(decays)), parameters(std::move(parameters))
{}
//! Generate an Event from the stored EventData.
/**
* @details Do jet clustering and classification.
* Use this to generate an Event.
*
* @note Calling this function destroys EventData
*
* @param jet_def Jet definition
* @param min_jet_pt minimal \f$p_T\f$ for each jet
*
* @returns Full clustered and classified event.
*/
Event cluster(
fastjet::JetDefinition const & jet_def, double const min_jet_pt);
//! Alias for cluster()
Event operator()(
fastjet::JetDefinition const & jet_def, double const min_jet_pt){
return cluster(jet_def, min_jet_pt);
}
//! Sort particles in rapidity
void sort();
//! Reconstruct intermediate particles from final-state leptons
/**
* Final-state leptons are created from virtual photons, W, or Z bosons.
* This function tries to reconstruct such intermediate bosons if they
* are not part of the event record.
*/
void reconstruct_intermediate();
//! Incoming particles
std::array<Particle, 2> incoming;
//! Outcoing particles
std::vector<Particle> outgoing;
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<std::size_t, std::vector<Particle>> decays;
//! Parameters, e.g. scale or inital weight
Parameters<EventParameters> parameters;
}; // end class EventData
//! Print Event
std::ostream& operator<<(std::ostream & os, Event const & ev);
//! Square of the partonic centre-of-mass energy \f$\hat{s}\f$
double shat(Event const & ev);
//! Convert an event to a LHEF::HEPEUP
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP *);
// put deprecated warning at the end, so don't get the warning inside Event.hh,
// additionally doxygen can not identify [[deprecated]] correctly
struct [[deprecated("UnclusteredEvent will be replaced by EventData")]]
UnclusteredEvent;
//! An event before jet clustering
//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0
struct UnclusteredEvent{
//! Default Constructor
UnclusteredEvent() = default;
//! Constructor from LesHouches event information
UnclusteredEvent(LHEF::HEPEUP const & hepeup);
std::array<Particle, 2> incoming; /**< Incoming Particles */
std::vector<Particle> outgoing; /**< Outgoing Particles */
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<std::size_t, std::vector<Particle>> decays;
//! Central parameter (e.g. scale) choice
EventParameters central;
std::vector<EventParameters> variations; /**< For parameter variation */
};
} // namespace HEJ
diff --git a/include/HEJ/EventReader.hh b/include/HEJ/EventReader.hh
index f4952bf..61614ae 100644
--- a/include/HEJ/EventReader.hh
+++ b/include/HEJ/EventReader.hh
@@ -1,60 +1,60 @@
/** \file
* \brief Header file for event reader interface
*
* This header defines an abstract base class for reading events from files.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <cstddef>
#include <memory>
#include <string>
#include "HEJ/optional.hh"
namespace LHEF {
class HEPEUP;
class HEPRUP;
}
namespace HEJ {
//! Abstract base class for reading events from files
struct EventReader {
//! Read an event
virtual bool read_event() = 0;
//! Access header text
virtual std::string const & header() const = 0;
//! Access run information
virtual LHEF::HEPRUP const & heprup() const = 0;
//! Access last read event
virtual LHEF::HEPEUP const & hepeup() const = 0;
//! Guess number of events from header
- virtual HEJ::optional<std::size_t> number_events() const {
+ virtual optional<std::size_t> number_events() const {
std::size_t start = header().rfind("Number of Events");
start = header().find_first_of("123456789", start);
if(start == std::string::npos) {
return {};
}
const std::size_t end = header().find_first_not_of("0123456789", start);
return std::stoi(header().substr(start, end - start));
}
virtual ~EventReader() = default;
};
//! Factory function for event readers
/**
* @param filename The name of the input file
* @returns A pointer to an instance of an EventReader
* for the input file
*/
std::unique_ptr<EventReader> make_reader(std::string const & filename);
} // namespace HEJ
diff --git a/include/HEJ/HDF5Reader.hh b/include/HEJ/HDF5Reader.hh
index 25bcc3b..0706548 100644
--- a/include/HEJ/HDF5Reader.hh
+++ b/include/HEJ/HDF5Reader.hh
@@ -1,51 +1,51 @@
/** \file
* \brief Header file for reading events in the HDF5 event format.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <string>
#include <memory>
#include "HEJ/EventReader.hh"
namespace HEJ{
//! Class for reading events from a file in the HDF5 file format
/**
* @details This format is specified in \cite Hoeche:2019rti.
*/
class HDF5Reader : public EventReader{
public:
HDF5Reader() = delete;
//! Contruct object reading from the given file
explicit HDF5Reader(std::string const & filename);
//! Read an event
bool read_event() override;
//! Access header text
std::string const & header() const override;
//! Access run information
LHEF::HEPRUP const & heprup() const override;
//! Access last read event
LHEF::HEPEUP const & hepeup() const override;
//! Get number of events
- HEJ::optional<size_t> number_events() const override;
+ optional<size_t> number_events() const override;
~HDF5Reader() override;
private:
struct HDF5ReaderImpl;
std::unique_ptr<HDF5ReaderImpl> impl_;
};
} // namespace HEJ
diff --git a/include/HEJ/LesHouchesReader.hh b/include/HEJ/LesHouchesReader.hh
index 8f11749..f83e9c8 100644
--- a/include/HEJ/LesHouchesReader.hh
+++ b/include/HEJ/LesHouchesReader.hh
@@ -1,82 +1,82 @@
/** \file
* \brief Header file for reading events in the Les Houches Event File format.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <string>
#include "LHEF/LHEF.h"
#include "HEJ/EventReader.hh"
#include "HEJ/stream.hh"
namespace HEJ{
//! Class for reading events from a file in the Les Houches Event File format
class LesHouchesReader : public EventReader{
public:
//! Contruct object reading from the given file
explicit LesHouchesReader(std::string const & filename):
stream_{filename},
reader_{stream_}
{
// always use the newest LHE version
reader_.heprup.version = LHEF::HEPRUP().version;
}
//! Read an event
bool read_event() override {
return reader_.readEvent();
}
//! Access header text
std::string const & header() const override {
return reader_.headerBlock;
}
//! Access run information
LHEF::HEPRUP const & heprup() const override {
return reader_.heprup;
}
//! Access last read event
LHEF::HEPEUP const & hepeup() const override {
return reader_.hepeup;
}
private:
- HEJ::istream stream_;
+ istream stream_;
protected:
//! Underlying reader
LHEF::Reader reader_;
};
/**
* @brief Les Houches Event file reader for LHE files created by Sherpa
* @details In Sherpa the cross section is given by
* sum(weights)/(number of trials). This EventReader converts the
* weights such that cross section=sum(weights)
* @note Reading from a pipe is not possible!
*/
class SherpaLHEReader : public LesHouchesReader{
public:
//! Inialise Reader for a Sherpa LHE file
explicit SherpaLHEReader(std::string const & filename);
bool read_event() override;
- HEJ::optional<size_t> number_events() const override {
+ optional<size_t> number_events() const override {
return num_events;
}
private:
double num_trials;
size_t num_events;
};
} // namespace HEJ
diff --git a/include/HEJ/RivetAnalysis.hh b/include/HEJ/RivetAnalysis.hh
index 9f28c78..a32748e 100644
--- a/include/HEJ/RivetAnalysis.hh
+++ b/include/HEJ/RivetAnalysis.hh
@@ -1,70 +1,70 @@
/** \file
* \brief HEJ 2 interface to rivet analyses
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#include <vector>
#include "LHEF/LHEF.h"
#include "HEJ/Analysis.hh"
#include "HEJ/optional.hh"
namespace YAML {
class Node;
}
namespace HEJ {
/**
* @brief Class representing a Rivet analysis
*
* This class inherits from Analysis and can therefore be used
* like any other HEJ 2 analysis.
*/
- class RivetAnalysis: public HEJ::Analysis {
+ class RivetAnalysis: public Analysis {
public:
//! Create RivetAnalysis
static std::unique_ptr<Analysis> create(
YAML::Node const & config, LHEF::HEPRUP const & heprup);
//! Constructor
/**
* @param config Configuration parameters
* @param heprup General run informations
*
* config["rivet"] should be the name of a single Rivet analysis or
* a list of Rivet analyses. config["output"] is the prefix for
* the .yoda output files.
*/
RivetAnalysis(YAML::Node const & config, LHEF::HEPRUP const & heprup);
~RivetAnalysis() override;
//! Pass an event to the underlying Rivet analysis
- void fill(HEJ::Event const & event, HEJ::Event const &) override;
- bool pass_cuts(HEJ::Event const &, HEJ::Event const &) override
+ void fill(Event const & event, Event const &) override;
+ bool pass_cuts(Event const &, Event const &) override
{return true;} //!< no additional cuts are applied
void finalise() override;
private:
std::vector<std::string> analyses_names_;
std::string output_name_;
LHEF::HEPRUP heprup_;
/// struct to organise the infos per rivet run/scale setting
struct rivet_info;
std::vector<rivet_info> rivet_runs_;
/**
* \internal
* @brief Calculates the scale variation from the first event for the output
* file
*/
- void init(HEJ::Event const & event);
+ void init(Event const & event);
bool first_event_;
};
} // namespace HEJ
diff --git a/include/HEJ/Wjets.hh b/include/HEJ/Wjets.hh
index dbe850e..021f3c1 100644
--- a/include/HEJ/Wjets.hh
+++ b/include/HEJ/Wjets.hh
@@ -1,509 +1,509 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in W+Jets.
*
* This file contains all the W+Jet specific components to compute
* the current contractions for valid HEJ processes, to form a full
* W+Jets ME, currently one would have to use functions from the
* jets.hh header also. We can calculate all subleading processes for
* W+Jets.
*/
#pragma once
#include <vector>
#include "CLHEP/Vector/LorentzVector.h"
namespace HEJ {
class ParticleProperties;
namespace currents {
typedef CLHEP::HepLorentzVector HLV;
//! Square of qQ->qenuQ W+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qQ->qenuQ
*
* This returns the square of the current contractions in qQ->qenuQ scattering
* with an emission of a W Boson.
*/
double ME_W_qQ(HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! Square of qbarQ->qbarenuQ W+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qbarQ->qbarenuQ
*
* This returns the square of the current contractions in qbarQ->qbarenuQ
* scattering with an emission of a W Boson.
*/
double ME_W_qbarQ(HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! Square of qQbar->qenuQbar W+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qQbar->qenuQbar
*
* This returns the square of the current contractions in qQbar->qenuQbar
* scattering with an emission of a W Boson.
*/
double ME_W_qQbar(HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! Square of qbarQbar->qbarenuQbar W+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for
* qbarQbar->qbarenuQbar
*
* This returns the square of the current contractions in
* qbarQbar->qbarenuQbar scattering with an emission of a W Boson.
*/
double ME_W_qbarQbar(HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in, HLV const & p2out,
HLV const & p2in,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! Square of qg->qenug W+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qg->qenug
*
* This returns the square of the current contractions in qg->qenug scattering
* with an emission of a W Boson.
*/
double ME_W_qg(HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! Square of qbarg->qbarenug W+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qbarg->qbarenug
*
* This returns the square of the current contractions in qbarg->qbarenug
* scattering with an emission of a W Boson.
*/
double ME_W_qbarg(HLV const & p1out, HLV
const & plbar, HLV const & pl,
HLV const & p1in, HLV const & p2out, HLV const & p2in,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! qQg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qQ->qQg Scattering
*
* This returns the square of the current contractions in qQg->qQg scattering
* with an emission of a W Boson.
*/
double ME_W_unob_qQ(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pg,
HLV const & plbar, HLV const & pl,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! qbarQg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qbarQ->qbarQg
*
* This returns the square of the current contractions in qbarQg->qbarQg
* scattering with an emission of a W Boson.
*/
double ME_W_unob_qbarQ(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pg,
HLV const & plbar, HLV const & pl,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! qQbarg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qQbar->qQbarg
*
* This returns the square of the current contractions in qQbarg->qQbarg
* scattering with an emission of a W Boson.
*/
double ME_W_unob_qQbar(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pg,
HLV const & plbar, HLV const & pl,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! qbarQbarg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qbarQbar->qbarQbarg
*
* This returns the square of the current contractions in qbarQbarg->qbarQbarg
* scattering with an emission of a W Boson.
*/
double ME_W_unob_qbarQbar(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in,
HLV const & pg,
HLV const & plbar, HLV const & pl,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+uno same leg
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qQ->qQg
*
* This returns the square of the current contractions in gqQ->gqQ scattering
* with an emission of a W Boson.
*/
double ME_Wuno_qQ(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pg,
HLV const & plbar, HLV const & pl,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+uno same leg. quark anti-quark
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qQbar->gqQbar
*
* This returns the square of the current contractions in gqQbar->gqQbar
* scattering with an emission of a W Boson. (Unordered Same Leg)
*/
double ME_Wuno_qQbar(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pg,
HLV const & plbar, HLV const & pl,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+uno same leg. quark gluon
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state gluon b
* @param p2in Momentum of intial state gluon b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qg->gqg
*
* This returns the square of the current contractions in qg->gqg scattering
* with an emission of a W Boson.
*/
double ME_Wuno_qg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pg,
HLV const & plbar, HLV const & pl,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+uno same leg. anti-quark quark
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qbarQ->gqbarQ
*
* This returns the square of the current contractions in qbarQ->gqbarQ
* scattering with an emission of a W Boson.
*/
double ME_Wuno_qbarQ(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pg,
HLV const & plbar, HLV const & pl,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+uno same leg. anti-quark anti-quark
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qbarQbar->gqbarQbar
*
* This returns the square of the current contractions in gqbarQbar->qbarQbar
* scattering with an emission of a W Boson.
*/
double ME_Wuno_qbarQbar(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in,
HLV const & pg,
HLV const & plbar, HLV const & pl,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+uno same leg. anti-quark gluon
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state gluon b
* @param p2in Momentum of intial state gluon b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for ->gqbarg
*
* This returns the square of the current contractions in qbarg->gqbarg
* scattering with an emission of a W Boson.
*/
double ME_Wuno_qbarg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pg,
HLV const & plbar, HLV const & pl,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+Extremal qqx. qxqQ
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state anti-quark b
* @param p2in Momentum of final state gluon b
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for ->qbarqQ
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WExqqx_qbarqQ(HLV const & pgin, HLV const & pqout,
HLV const & plbar, HLV const & pl,
HLV const & pqbarout, HLV const & p2out,
HLV const & p2in,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+Extremal qqx. qqxQ
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state anti-quark b
* @param p2in Momentum of final state gluon b
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for ->qqbarQ
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WExqqx_qqbarQ(HLV const & pgin, HLV const & pqout,
HLV const & plbar, HLV const & pl,
HLV const & pqbarout, HLV const & p2out,
HLV const & p2in,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+Extremal qqx. gg->qxqg
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state gluon b
* @param p2in Momentum of final state gluon b
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for gg->qbarqg
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WExqqx_qbarqg(HLV const & pgin, HLV const & pqout,
HLV const & plbar, HLV const & pl,
HLV const & pqbarout, HLV const & p2out,
HLV const & p2in,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+Extremal qqx. gg->qqxg
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state gluon a
* @param p2in Momentum of final state gluon b
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for gg->qqbarg
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WExqqx_qqbarg(HLV const & pgin, HLV const & pqbarout,
HLV const & plbar, HLV const & pl,
HLV const & pqout, HLV const & p2out,
HLV const & p2in,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+Extremal qqx. gg->qqxg. qqx on forwards leg, W emission backwards leg.
/**
* @param pa Momentum of initial state (anti-)quark
* @param pb Momentum of initial state gluon
* @param p1 Momentum of final state (anti-)quark (after W emission)
* @param p2 Momentum of final state anti-quark
* @param p3 Momentum of final state quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param aqlinepa Is opposite extremal leg to qqx a quark or antiquark line
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for gq->qqbarqW
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated via current contraction of existing
* currents. Assumes qqx split from forwards leg, W emission from backwards
* leg. Switch input (pa<->pb, p1<->pn) if calculating forwards qqx.
*/
double ME_W_Exqqx_QQq(HLV const & pa, HLV const & pb, HLV const & p1,
HLV const & p2, HLV const & p3,
HLV const & plbar, HLV const & pl, bool aqlinepa,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+Jets qqxCentral. qqx W emission.
/**
* @param pa Momentum of initial state particle a
* @param pb Momentum of initial state particle b
* @param pl Momentum of final state lepton
* @param plbar Momentum of final state anti-lepton
* @param partons Vector of outgoing parton momenta
* @param aqlinepa True= pa is anti-quark
* @param aqlinepb True= pb is anti-quark
* @param qqxmarker Ordering of the qqbar pair produced (qqx vs qxq)
* @param nabove Number of lipatov vertices "above" qqbar pair
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qq>qQQbarWq
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WCenqqx_qq(HLV const & pa, HLV const & pb,
HLV const & pl, HLV const & plbar,
std::vector<HLV> const & partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
//! W+Jets qqxCentral. W emission from backwards leg.
/**
* @param ka Momentum of initial state particle a
* @param kb Momentum of initial state particle b
* @param pl Momentum of final state lepton
* @param plbar Momentum of final state anti-lepton
* @param partons outgoing parton momenta
* @param aqlinepa True= pa is anti-quark
* @param aqlinepb True= pb is anti-quark
* @param qqxmarker Ordering of the qqbar pair produced (qqx vs qxq)
* @param nabove Number of lipatov vertices "above" qqbar pair
* @param nbelow Number of lipatov vertices "below" qqbar pair
* @param forwards Swap to emission off front leg
* TODO: remove so args can be const
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qq>qQQbarWq
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_W_Cenqqx_qq(HLV ka, HLV kb,
HLV pl, HLV plbar,
std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker,
int nabove, int nbelow, bool forwards,
- HEJ::ParticleProperties const & wprop);
+ ParticleProperties const & wprop);
} // namespace currents
} // namespace HEJ
diff --git a/include/HEJ/Zjets.hh b/include/HEJ/Zjets.hh
index 4ebcd48..0292b0a 100644
--- a/include/HEJ/Zjets.hh
+++ b/include/HEJ/Zjets.hh
@@ -1,127 +1,127 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in Z+Jets.
*/
#pragma once
#include <vector>
#include "CLHEP/Vector/LorentzVector.h"
#include "HEJ/PDG_codes.hh"
namespace HEJ {
class ParticleProperties;
namespace currents {
typedef CLHEP::HepLorentzVector HLV;
//! Square of qQ->qQe+e- Z+Jets Scattering Current
/**
* @param pa Momentum of initial state quark
* @param pb Momentum of initial state quark
* @param p1 Momentum of final state quark
* @param p2 Momentum of final state quark
* @param pep Momentum of final state positron
* @param pem Momentum of final state electron
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
* @param zprop Mass and width of the Z boson
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns Square of the current contractions for qQ->qQe+e- Scattering
*
* This returns the square of the current contractions in qQ->qQ scattering
* with an emission of a Z Boson.
*/
std::vector <double> ME_Z_qQ(const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem,
- const HEJ::ParticleID aptype, const HEJ::ParticleID bptype,
- HEJ::ParticleProperties const & zprop,
+ const ParticleID aptype, const ParticleID bptype,
+ ParticleProperties const & zprop,
const double stw2, const double ctw);
//! Square of qg->qge+e- Z+Jets Scattering Current
/**
* @param pa Momentum of initial state quark
* @param pb Momentum of initial state gluon
* @param p1 Momentum of final state quark
* @param p2 Momentum of final state gluon
* @param pep Momentum of final state positron
* @param pem Momentum of final state electron
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
* @param zprop Mass and width of the Z boson
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns Square of the current contractions for qg->qge+e- Scattering
*
* This returns the square of the current contractions in qg->qg scattering
* with an emission of a Z Boson.
*/
double ME_Z_qg(const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem,
- const HEJ::ParticleID aptype, const HEJ::ParticleID bptype,
- HEJ::ParticleProperties const & zprop,
+ const ParticleID aptype, const ParticleID bptype,
+ ParticleProperties const & zprop,
const double stw2, const double ctw);
//! Square of qQ->gqQe+e- Z+Jets Unordered Current
/**
* @param pa Momentum of initial state quark a
* @param pb Momentum of initial state quark b
* @param pg Momentum of final state unordered gluon
* @param p1 Momentum of final state quark a
* @param p2 Momentum of final state quark b
* @param pep Momentum of final state positron
* @param pem Momentum of final state electron
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
* @param zprop Mass and width of the Z boson
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns Square of the current contractions for qQ->gqQe+e- Scattering
*
* This returns the square of the current contractions in qQ->gqQ scattering
* with an emission of a Z Boson.
*/
std::vector <double> ME_Zuno_qQ(const HLV & pa, const HLV & pb,
const HLV & pg, const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem,
- const HEJ::ParticleID aptype, const HEJ::ParticleID bptype,
- HEJ::ParticleProperties const & zprop,
+ const ParticleID aptype, const ParticleID bptype,
+ ParticleProperties const & zprop,
const double stw2, const double ctw);
//! Square of qg->gqge+e- Z+Jets Unordered Current
/**
* @param pa Momentum of initial state quark
* @param pb Momentum of initial state gluon
* @param pg Momentum of final state unordered gluon
* @param p1 Momentum of final state quark
* @param p2 Momentum of final state gluon
* @param pep Momentum of final state positron
* @param pem Momentum of final state electron
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
* @param zprop Mass and width of the Z boson
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns Square of the current contractions for qg->gqge+e- Scattering
*
* This returns the square of the current contractions in qg->gqg scattering
* with an emission of a Z Boson.
*/
double ME_Zuno_qg(const HLV & pa, const HLV & pb,
const HLV & pg, const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem,
- const HEJ::ParticleID aptype, const HEJ::ParticleID bptype,
- HEJ::ParticleProperties const & zprop,
+ const ParticleID aptype, const ParticleID bptype,
+ ParticleProperties const & zprop,
const double stw2, const double ctw);
} // namespace currents
} // namespace HEJ
diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index aa4c057..e311c28 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,257 +1,257 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/EventReweighter.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <functional>
#include <string>
#include <unordered_map>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/PseudoJet.hh"
#include "LHEF/LHEF.h"
#include "HEJ/Fraction.hh"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/PhaseSpacePoint.hh"
namespace HEJ{
EventReweighter::EventReweighter(
LHEF::HEPRUP const & heprup,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
std::shared_ptr<RNG> ran
):
EventReweighter{
- HEJ::Beam{
+ Beam{
heprup.EBMUP.first,
{{
- static_cast<HEJ::ParticleID>(heprup.IDBMUP.first),
- static_cast<HEJ::ParticleID>(heprup.IDBMUP.second)
+ static_cast<ParticleID>(heprup.IDBMUP.first),
+ static_cast<ParticleID>(heprup.IDBMUP.second)
}}
},
heprup.PDFSUP.first,
std::move(scale_gen),
std::move(conf),
std::move(ran)
}
{
if(heprup.EBMUP.second != E_beam_){
throw std::invalid_argument(
"asymmetric beam: " + std::to_string(E_beam_)
+ " ---> <--- " + std::to_string(heprup.EBMUP.second)
);
}
if(heprup.PDFSUP.second != pdf_.id()){
throw std::invalid_argument(
"conflicting PDF ids: " + std::to_string(pdf_.id())
+ " vs. " + std::to_string(heprup.PDFSUP.second)
);
}
}
EventReweighter::EventReweighter(
Beam const & beam,
int pdf_id,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
std::shared_ptr<RNG> ran
):
param_{std::move(conf)},
E_beam_{beam.E},
pdf_{pdf_id, beam.type.front(), beam.type.back()},
MEt2_{
[this](double mu){ return pdf_.Halphas(mu); },
param_.ME_config
},
scale_gen_{std::move(scale_gen)},
ran_{std::move(ran)}
{
assert(ran_);
}
PDF const & EventReweighter::pdf() const{
return pdf_;
}
std::vector<Event> EventReweighter::reweight(
Event const & input_ev, std::size_t num_events
){
auto res_events{ gen_res_events(input_ev, num_events) };
if(res_events.empty()) return {};
for(auto & event: res_events) event = scale_gen_(std::move(event));
return rescale(input_ev, std::move(res_events));
}
EventTreatment EventReweighter::treatment(EventType type) const {
return param_.treat.at(type);
}
std::vector<Event> EventReweighter::gen_res_events(
Event const & ev,
std::size_t phase_space_points
){
assert(ev.variations().empty());
status_.clear();
switch(treatment(ev.type())){
case EventTreatment::discard: {
status_.emplace_back(StatusCode::discard);
return {};
}
case EventTreatment::keep:
if(! jets_pass_resummation_cuts(ev)) {
status_.emplace_back(StatusCode::failed_resummation_cuts);
return {};
}
else {
status_.emplace_back(StatusCode::good);
return {ev};
}
default:;
}
const double Born_shat = shat(ev);
std::vector<Event> resummation_events;
status_.reserve(phase_space_points);
for(std::size_t psp_number = 0; psp_number < phase_space_points; ++psp_number){
PhaseSpacePoint psp{ev, param_.psp_config, *ran_};
status_.emplace_back(psp.status());
assert(psp.status() != StatusCode::unspecified);
if(psp.status() != StatusCode::good) continue;
assert(psp.weight() != 0.);
if(psp.incoming()[0].E() > E_beam_ || psp.incoming()[1].E() > E_beam_) {
status_.back() = StatusCode::too_much_energy;
continue;
}
resummation_events.emplace_back(
to_EventData( std::move(psp) ).cluster(
param_.jet_param().def, param_.jet_param().min_pt
)
);
auto & new_event = resummation_events.back();
assert( new_event.valid_hej_state(
param_.psp_config.max_ext_soft_pt_fraction,
param_.psp_config.min_extparton_pt ) );
if( new_event.type() != ev.type() )
throw std::logic_error{"Resummation Event does not match Born event"};
new_event.generate_colours(*ran_);
assert(new_event.variations().empty());
new_event.central().mur = ev.central().mur;
new_event.central().muf = ev.central().muf;
const double resum_shat = shat(new_event);
new_event.central().weight *= ev.central().weight*Born_shat*Born_shat/
(phase_space_points*resum_shat*resum_shat);
}
return resummation_events;
}
std::vector<Event> EventReweighter::rescale(
Event const & Born_ev,
std::vector<Event> events
) const{
const double Born_pdf = pdf_factors(Born_ev).central;
const double Born_ME = tree_matrix_element(Born_ev);
for(auto & cur_event: events){
const auto pdf = pdf_factors(cur_event);
assert(pdf.variations.size() == cur_event.variations().size());
const auto ME = matrix_elements(cur_event);
assert(ME.variations.size() == cur_event.variations().size());
cur_event.parameters() *= pdf*ME/(Born_pdf*Born_ME);
}
return events;
}
bool EventReweighter::jets_pass_resummation_cuts(
Event const & ev
) const{
const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
fastjet::ClusterSequence cs{out_as_PseudoJet, param_.jet_param().def};
return cs.inclusive_jets(param_.jet_param().min_pt).size() == ev.jets().size();
}
Weights EventReweighter::pdf_factors(Event const & ev) const{
auto const & a = ev.incoming().front();
auto const & b = ev.incoming().back();
const double xa = a.p.e()/E_beam_;
const double xb = b.p.e()/E_beam_;
Weights result;
std::unordered_map<double, double> known_pdf;
result.central =
pdf_.pdfpt(0,xa,ev.central().muf,a.type)*
pdf_.pdfpt(1,xb,ev.central().muf,b.type);
known_pdf.emplace(ev.central().muf, result.central);
result.variations.reserve(ev.variations().size());
for(auto const & ev_param: ev.variations()){
const double muf = ev_param.muf;
auto cur_pdf = known_pdf.find(muf);
if(cur_pdf == known_pdf.end()){
cur_pdf = known_pdf.emplace(
muf,
pdf_.pdfpt(0,xa,muf,a.type)*pdf_.pdfpt(1,xb,muf,b.type)
).first;
}
result.variations.emplace_back(cur_pdf->second);
}
assert(result.variations.size() == ev.variations().size());
return result;
}
Weights
EventReweighter::matrix_elements(Event const & ev) const{
assert(param_.treat.count(ev.type()) > 0);
if(param_.treat.find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev);
}
return MEt2_(ev);
}
double EventReweighter::tree_matrix_element(Event const & ev) const{
assert(ev.variations().empty());
assert(param_.treat.count(ev.type()) > 0);
if(param_.treat.find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev).central;
}
return MEt2_.tree(ev).central;
}
Weights
EventReweighter::fixed_order_scale_ME(Event const & ev) const{
int alpha_s_power = 0;
for(auto const & part: ev.outgoing()){
if(is_parton(part))
++alpha_s_power;
else if(part.type == pid::Higgs) {
alpha_s_power += 2;
}
// nothing to do for other uncoloured particles
}
Weights result;
result.central = std::pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
for(auto const & var: ev.variations()){
result.variations.emplace_back(
std::pow(pdf_.Halphas(var.mur), alpha_s_power)
);
}
return result;
}
} // namespace HEJ
diff --git a/src/HDF5Reader.cc b/src/HDF5Reader.cc
index 0c3459f..e040a78 100644
--- a/src/HDF5Reader.cc
+++ b/src/HDF5Reader.cc
@@ -1,322 +1,322 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/HDF5Reader.hh"
#include "HEJ/ConfigFlags.hh"
#ifdef HEJ_BUILD_WITH_HDF5
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <functional>
#include <iterator>
#include <numeric>
#include <utility>
#include <vector>
#include "highfive/H5File.hpp"
#include "LHEF/LHEF.h"
#else
#include "HEJ/exceptions.hh"
#endif
#ifdef HEJ_BUILD_WITH_HDF5
namespace HEJ {
namespace {
// buffer size for reader
// each "reading from disk" reads "chunk_size" many event at once
constexpr std::size_t chunk_size = 10000;
struct ParticleData {
std::vector<int> id;
std::vector<int> status;
std::vector<int> mother1;
std::vector<int> mother2;
std::vector<int> color1;
std::vector<int> color2;
std::vector<double> px;
std::vector<double> py;
std::vector<double> pz;
std::vector<double> e;
std::vector<double> m;
std::vector<double> lifetime;
std::vector<double> spin;
};
struct EventRecords {
std::vector<int> particle_start;
std::vector<int> nparticles;
std::vector<int> pid;
std::vector<double> weight;
std::vector<double> scale;
std::vector<double> fscale;
std::vector<double> rscale;
std::vector<double> aqed;
std::vector<double> aqcd;
double trials;
ParticleData particles;
};
class ConstEventIterator {
public:
// iterator traits
using iterator_category = std::bidirectional_iterator_tag;
using value_type = LHEF::HEPEUP;
using difference_type = std::ptrdiff_t;
using pointer = const LHEF::HEPEUP*;
using reference = LHEF::HEPEUP const &;
using iterator = ConstEventIterator;
friend iterator cbegin(EventRecords const & records) noexcept;
friend iterator cend(EventRecords const & records) noexcept;
iterator& operator++() {
particle_offset_ += records_.get().nparticles[idx_];
++idx_;
return *this;
}
iterator& operator--() {
--idx_;
particle_offset_ -= records_.get().nparticles[idx_];
return *this;
}
iterator operator--(int) {
auto res = *this;
--(*this);
return res;
}
bool operator==(iterator const & other) const {
return idx_ == other.idx_;
}
bool operator!=(iterator other) const {
return !(*this == other);
}
value_type operator*() const {
value_type hepeup{};
auto const & r = records_.get();
hepeup.NUP = r.nparticles[idx_];
hepeup.IDPRUP = r.pid[idx_];
hepeup.XWGTUP = r.weight[idx_]/r.trials;
hepeup.weights.emplace_back(hepeup.XWGTUP, nullptr);
hepeup.SCALUP = r.scale[idx_];
hepeup.scales.muf = r.fscale[idx_];
hepeup.scales.mur = r.rscale[idx_];
hepeup.AQEDUP = r.aqed[idx_];
hepeup.AQCDUP = r.aqcd[idx_];
const size_t start = particle_offset_;
const size_t end = start + hepeup.NUP;
auto const & p = r.particles;
hepeup.IDUP = std::vector<long>( begin(p.id)+start, begin(p.id)+end );
hepeup.ISTUP = std::vector<int>( begin(p.status)+start, begin(p.status)+end );
hepeup.VTIMUP = std::vector<double>( begin(p.lifetime)+start, begin(p.lifetime)+end );
hepeup.SPINUP = std::vector<double>( begin(p.spin)+start, begin(p.spin)+end );
hepeup.MOTHUP.resize(hepeup.NUP);
hepeup.ICOLUP.resize(hepeup.NUP);
hepeup.PUP.resize(hepeup.NUP);
for(size_t i = 0; i < hepeup.MOTHUP.size(); ++i) {
const size_t idx = start + i;
assert(idx < end);
hepeup.MOTHUP[i] = std::make_pair(p.mother1[idx], p.mother2[idx]);
hepeup.ICOLUP[i] = std::make_pair(p.color1[idx], p.color2[idx]);
hepeup.PUP[i] = std::vector<double>{
p.px[idx], p.py[idx], p.pz[idx], p.e[idx], p.m[idx]
};
}
return hepeup;
}
private:
explicit ConstEventIterator(EventRecords const & records):
records_{records} {}
std::reference_wrapper<const EventRecords> records_;
size_t idx_ = 0;
size_t particle_offset_ = 0;
}; // end ConstEventIterator
ConstEventIterator cbegin(EventRecords const & records) noexcept {
return ConstEventIterator{records};
}
ConstEventIterator cend(EventRecords const & records) noexcept {
auto it =ConstEventIterator{records};
it.idx_ = records.aqcd.size(); // or size of any other records member
return it;
}
} // end anonymous namespace
struct HDF5Reader::HDF5ReaderImpl{
HighFive::File file;
std::size_t event_idx;
std::size_t particle_idx;
std::size_t nevents;
EventRecords records;
ConstEventIterator cur_event;
LHEF::HEPRUP heprup;
LHEF::HEPEUP hepeup;
explicit HDF5ReaderImpl(std::string const & filename):
file{filename},
event_idx{0},
particle_idx{0},
nevents{
file.getGroup("event")
.getDataSet("nparticles") // or any other dataset
.getSpace().getDimensions().front()
},
records{},
cur_event{cbegin(records)},
heprup{},
hepeup{}
{
read_heprup();
read_event_records(chunk_size);
}
void read_heprup() {
const auto init = file.getGroup("init");
init.getDataSet( "PDFgroupA" ).read(heprup.PDFGUP.first);
init.getDataSet( "PDFgroupB" ).read(heprup.PDFGUP.second);
init.getDataSet( "PDFsetA" ).read(heprup.PDFSUP.first);
init.getDataSet( "PDFsetB" ).read(heprup.PDFSUP.second);
init.getDataSet( "beamA" ).read(heprup.IDBMUP.first);
init.getDataSet( "beamB" ).read(heprup.IDBMUP.second);
init.getDataSet( "energyA" ).read(heprup.EBMUP.first);
init.getDataSet( "energyB" ).read(heprup.EBMUP.second);
init.getDataSet( "numProcesses" ).read(heprup.NPRUP);
init.getDataSet( "weightingStrategy" ).read(heprup.IDWTUP);
const auto proc_info = file.getGroup("procInfo");
proc_info.getDataSet( "procId" ).read(heprup.LPRUP);
proc_info.getDataSet( "xSection" ).read(heprup.XSECUP);
proc_info.getDataSet( "error" ).read(heprup.XERRUP);
// TODO: is this identification correct?
proc_info.getDataSet( "unitWeight" ).read(heprup.XMAXUP);
std::vector<double> trials;
file.getGroup("event").getDataSet("trials").read(trials);
records.trials = std::accumulate(begin(trials), end(trials), 0.);
}
std::size_t read_event_records(std::size_t count) {
count = std::min(count, nevents-event_idx);
auto events = file.getGroup("event");
events.getDataSet("nparticles").select({event_idx}, {count}).read(records.nparticles);
assert(records.nparticles.size() == count);
events.getDataSet("pid").select( {event_idx}, {count} ).read( records.pid );
events.getDataSet("weight").select( {event_idx}, {count} ).read( records.weight );
events.getDataSet("scale").select( {event_idx}, {count} ).read( records.scale );
events.getDataSet("fscale").select( {event_idx}, {count} ).read( records.fscale );
events.getDataSet("rscale").select( {event_idx}, {count} ).read( records.rscale );
events.getDataSet("aqed").select( {event_idx}, {count} ).read( records.aqed );
events.getDataSet("aqcd").select( {event_idx}, {count} ).read( records.aqcd );
const std::size_t particle_count = std::accumulate(
begin(records.nparticles), end(records.nparticles), 0
);
auto pdata = file.getGroup("particle");
auto & particles = records.particles;
pdata.getDataSet("id").select( {particle_idx}, {particle_count} ).read( particles.id );
pdata.getDataSet("status").select( {particle_idx}, {particle_count} ).read( particles.status );
pdata.getDataSet("mother1").select( {particle_idx}, {particle_count} ).read( particles.mother1 );
pdata.getDataSet("mother2").select( {particle_idx}, {particle_count} ).read( particles.mother2 );
pdata.getDataSet("color1").select( {particle_idx}, {particle_count} ).read( particles.color1 );
pdata.getDataSet("color2").select( {particle_idx}, {particle_count} ).read( particles.color2 );
pdata.getDataSet("px").select( {particle_idx}, {particle_count} ).read( particles.px );
pdata.getDataSet("py").select( {particle_idx}, {particle_count} ).read( particles.py );
pdata.getDataSet("pz").select( {particle_idx}, {particle_count} ).read( particles.pz );
pdata.getDataSet("e").select( {particle_idx}, {particle_count} ).read( particles.e );
pdata.getDataSet("m").select( {particle_idx}, {particle_count} ).read( particles.m );
pdata.getDataSet("lifetime").select( {particle_idx}, {particle_count} ).read( particles.lifetime );
pdata.getDataSet("spin").select( {particle_idx}, {particle_count} ).read( particles.spin );
event_idx += count;
particle_idx += particle_count;
return count;
}
};
HDF5Reader::HDF5Reader(std::string const & filename):
impl_{std::make_unique<HDF5ReaderImpl>(filename)}
{}
bool HDF5Reader::read_event() {
if(impl_->cur_event == cend(impl_->records)) {
// end of active chunk, read new events from file
const auto events_read = impl_->read_event_records(chunk_size);
impl_->cur_event = cbegin(impl_->records);
if(events_read == 0) return false;
}
impl_->hepeup = *impl_->cur_event;
++impl_->cur_event;
return true;
}
namespace {
static const std::string nothing = "";
}
std::string const & HDF5Reader::header() const {
return nothing;
}
LHEF::HEPRUP const & HDF5Reader::heprup() const {
return impl_->heprup;
}
LHEF::HEPEUP const & HDF5Reader::hepeup() const {
return impl_->hepeup;
}
- HEJ::optional<size_t> HDF5Reader::number_events() const {
+ optional<size_t> HDF5Reader::number_events() const {
return impl_->nevents;
}
}
#else // no HDF5 support
namespace HEJ {
class HDF5Reader::HDF5ReaderImpl{};
HDF5Reader::HDF5Reader(std::string const &){
throw std::invalid_argument{
"Failed to create HDF5 reader: "
"HEJ 2 was built without HDF5 support"
};
}
bool HDF5Reader::read_event() {
throw std::logic_error{"unreachable"};
}
std::string const & HDF5Reader::header() const {
throw std::logic_error{"unreachable"};
}
LHEF::HEPRUP const & HDF5Reader::heprup() const {
throw std::logic_error{"unreachable"};
}
LHEF::HEPEUP const & HDF5Reader::hepeup() const {
throw std::logic_error{"unreachable"};
}
- HEJ::optional<size_t> HDF5Reader::number_events() const {
+ optional<size_t> HDF5Reader::number_events() const {
throw std::logic_error{"unreachable"};
}
}
#endif
namespace HEJ {
HDF5Reader::~HDF5Reader() = default;
}
diff --git a/src/HDF5Writer.cc b/src/HDF5Writer.cc
index 5bd43b2..1a3f96d 100644
--- a/src/HDF5Writer.cc
+++ b/src/HDF5Writer.cc
@@ -1,417 +1,417 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/HDF5Writer.hh"
#include <cassert>
#include "LHEF/LHEF.h"
#include "HEJ/ConfigFlags.hh"
#ifdef HEJ_BUILD_WITH_HDF5
#include <cmath>
#include <cstddef>
#include <iterator>
#include <type_traits>
#include <utility>
#include <vector>
#include "highfive/H5File.hpp"
#include "HEJ/Event.hh"
#include "HEJ/event_types.hh"
#else
#include "HEJ/exceptions.hh"
#endif
#ifdef HEJ_BUILD_WITH_HDF5
namespace HEJ{
using HighFive::File;
using HighFive::DataSpace;
namespace {
using std::size_t;
constexpr size_t chunk_size = 1000;
constexpr unsigned compression_level = 3;
size_t to_index(event_type::EventType const type){
return type==0?0:std::floor(std::log2(static_cast<size_t>(type)))+1;
}
template<typename T>
void write_dataset(HighFive::Group & group, std::string const & name, T val) {
using data_t = std::decay_t<T>;
group.createDataSet<data_t>(name, DataSpace::From(val)).write(val);
}
template<typename T>
void write_dataset(
HighFive::Group & group, std::string const & name,
std::vector<T> const & val
) {
using data_t = std::decay_t<T>;
group.createDataSet<data_t>(name, DataSpace::From(val)).write(val);
}
struct Cache {
explicit Cache(size_t capacity):
capacity{capacity}
{
nparticles.reserve(capacity);
start.reserve(capacity);
pid.reserve(capacity);
weight.reserve(capacity);
scale.reserve(capacity);
fscale.reserve(capacity);
rscale.reserve(capacity);
aqed.reserve(capacity);
aqcd.reserve(capacity);
trials.reserve(capacity);
npLO.reserve(capacity);
npNLO.reserve(capacity);
}
- void fill(HEJ::Event const & ev) {
+ void fill(Event const & ev) {
const auto hepeup = to_HEPEUP(ev, nullptr);
// HEJ event to get nice wrapper
const auto num_partons = std::distance(ev.cbegin_partons(),
ev.cend_partons());
assert(num_partons>0);
// Number of partons for LO matching, HEJ requires at least 2 partons
npLO.emplace_back(num_partons>1?num_partons-2:num_partons);
// Number of real emissions in NLO, HEJ is LO -> -1
npNLO.emplace_back(-1);
fill_event_params(hepeup);
fill_event_particles(hepeup);
}
void fill_event_params(LHEF::HEPEUP const & ev) {
nparticles.emplace_back(ev.NUP);
start.emplace_back(particle_pos);
pid.emplace_back(ev.IDPRUP);
weight.emplace_back(ev.XWGTUP);
scale.emplace_back(ev.SCALUP);
fscale.emplace_back(ev.scales.muf);
rscale.emplace_back(ev.scales.mur);
aqed.emplace_back(ev.AQEDUP);
aqcd.emplace_back(ev.AQCDUP);
// set first trial=1 for first event
// -> sum(trials) = 1 -> xs=sum(weights)/sum(trials) as in Sherpa
if(particle_pos == 0){
trials.emplace_back(1.);
} else {
trials.emplace_back(0.);
}
particle_pos += ev.NUP;
}
void fill_event_particles(LHEF::HEPEUP const & ev) {
id.insert(end(id), begin(ev.IDUP), end(ev.IDUP));
status.insert(end(status), begin(ev.ISTUP), end(ev.ISTUP));
lifetime.insert(end(lifetime), begin(ev.VTIMUP), end(ev.VTIMUP));
spin.insert(end(spin), begin(ev.SPINUP), end(ev.SPINUP));
for(int i = 0; i < ev.NUP; ++i) {
mother1.emplace_back(ev.MOTHUP[i].first);
mother2.emplace_back(ev.MOTHUP[i].second);
color1.emplace_back(ev.ICOLUP[i].first);
color2.emplace_back(ev.ICOLUP[i].second);
px.emplace_back(ev.PUP[i][0]);
py.emplace_back(ev.PUP[i][1]);
pz.emplace_back(ev.PUP[i][2]);
e.emplace_back(ev.PUP[i][3]);
m.emplace_back(ev.PUP[i][4]);
}
}
bool is_full() const {
return nparticles.size() >= capacity;
}
void clear() {
nparticles.clear();
start.clear();
pid.clear();
id.clear();
status.clear();
mother1.clear();
mother2.clear();
color1.clear();
color2.clear();
weight.clear();
scale.clear();
fscale.clear();
rscale.clear();
aqed.clear();
aqcd.clear();
trials.clear();
npLO.clear();
npNLO.clear();
px.clear();
py.clear();
pz.clear();
e.clear();
m.clear();
lifetime.clear();
spin.clear();
}
size_t capacity;
std::vector<int> nparticles, start, pid, id, status,
mother1, mother2, color1, color2, npLO, npNLO;
std::vector<double> weight, scale, fscale, rscale, aqed,
aqcd, trials, px, py, pz, e, m, lifetime, spin;
private:
size_t particle_pos = 0;
};
}
struct HDF5Writer::HDF5WriterImpl{
File file;
LHEF::HEPRUP heprup;
Cache cache{chunk_size};
size_t event_idx = 0;
size_t particle_idx = 0;
HDF5WriterImpl(std::string const & filename, LHEF::HEPRUP && hepr):
file{filename, File::ReadWrite | File::Create | File::Truncate},
heprup{std::move(hepr)}
{
// TODO: code duplication with Les Houches Writer
const int max_number_types = to_index(event_type::last_type)+1;
heprup.NPRUP = max_number_types;
// ids of event types
heprup.LPRUP.clear();
heprup.LPRUP.reserve(max_number_types);
heprup.LPRUP.emplace_back(0);
for(size_t i=event_type::first_type+1; i<=event_type::last_type; i*=2) {
heprup.LPRUP.emplace_back(i);
}
heprup.XSECUP = std::vector<double>(max_number_types, 0.);
heprup.XERRUP = std::vector<double>(max_number_types, 0.);
heprup.XMAXUP = std::vector<double>(max_number_types, 0.);
write_init();
create_event_group();
create_particle_group();
}
void write_init() {
auto init = file.createGroup("init");
write_dataset(init, "PDFgroupA" , heprup.PDFGUP.first);
write_dataset(init, "PDFgroupB" , heprup.PDFGUP.second);
write_dataset(init, "PDFsetA" , heprup.PDFSUP.first);
write_dataset(init, "PDFsetB" , heprup.PDFSUP.second);
write_dataset(init, "beamA" , heprup.IDBMUP.first);
write_dataset(init, "beamB" , heprup.IDBMUP.second);
write_dataset(init, "energyA" , heprup.EBMUP.first);
write_dataset(init, "energyB" , heprup.EBMUP.second);
write_dataset(init, "numProcesses" , heprup.NPRUP);
write_dataset(init, "weightingStrategy", heprup.IDWTUP);
auto proc_info = file.createGroup("procInfo");
write_dataset(proc_info, "procId", heprup.LPRUP);
}
static HighFive::DataSetCreateProps const & hdf5_chunk() {
static const auto props = [](){
HighFive::DataSetCreateProps props;
props.add(HighFive::Chunking({chunk_size}));
props.add(HighFive::Deflate(compression_level));
return props;
}();
return props;
}
void create_event_group() {
static const auto dim = DataSpace({0}, {DataSpace::UNLIMITED});
auto events = file.createGroup("event");
events.createDataSet<int>("nparticles", dim, hdf5_chunk());
events.createDataSet<int>("start", dim, hdf5_chunk());
events.createDataSet<int>("pid", dim, hdf5_chunk());
events.createDataSet<double>("weight", dim, hdf5_chunk());
events.createDataSet<double>("scale", dim, hdf5_chunk());
events.createDataSet<double>("fscale", dim, hdf5_chunk());
events.createDataSet<double>("rscale", dim, hdf5_chunk());
events.createDataSet<double>("aqed", dim, hdf5_chunk());
events.createDataSet<double>("aqcd", dim, hdf5_chunk());
events.createDataSet<double>("trials", dim, hdf5_chunk());
events.createDataSet<int>("npLO", dim, hdf5_chunk());
events.createDataSet<int>("npNLO", dim, hdf5_chunk());
}
void resize_event_group(size_t new_size) {
auto events = file.getGroup("event");
events.getDataSet("nparticles").resize({new_size});
events.getDataSet("start").resize({new_size});
events.getDataSet("pid").resize({new_size});
events.getDataSet("weight").resize({new_size});
events.getDataSet("scale").resize({new_size});
events.getDataSet("fscale").resize({new_size});
events.getDataSet("rscale").resize({new_size});
events.getDataSet("aqed").resize({new_size});
events.getDataSet("aqcd").resize({new_size});
events.getDataSet("trials").resize({new_size});
events.getDataSet("npLO").resize({new_size});
events.getDataSet("npNLO").resize({new_size});
}
void create_particle_group() {
static const auto dim = DataSpace({0}, {DataSpace::UNLIMITED});
auto particles = file.createGroup("particle");
particles.createDataSet<int>("id", dim, hdf5_chunk());
particles.createDataSet<int>("status", dim, hdf5_chunk());
particles.createDataSet<int>("mother1", dim, hdf5_chunk());
particles.createDataSet<int>("mother2", dim, hdf5_chunk());
particles.createDataSet<int>("color1", dim, hdf5_chunk());
particles.createDataSet<int>("color2", dim, hdf5_chunk());
particles.createDataSet<double>("px", dim, hdf5_chunk());
particles.createDataSet<double>("py", dim, hdf5_chunk());
particles.createDataSet<double>("pz", dim, hdf5_chunk());
particles.createDataSet<double>("e", dim, hdf5_chunk());
particles.createDataSet<double>("m", dim, hdf5_chunk());
particles.createDataSet<double>("lifetime", dim, hdf5_chunk());
particles.createDataSet<double>("spin", dim, hdf5_chunk());
}
void resize_particle_group(size_t new_size) {
auto particles = file.getGroup("particle");
particles.getDataSet("id").resize({new_size});
particles.getDataSet("status").resize({new_size});
particles.getDataSet("mother1").resize({new_size});
particles.getDataSet("mother2").resize({new_size});
particles.getDataSet("color1").resize({new_size});
particles.getDataSet("color2").resize({new_size});
particles.getDataSet("px").resize({new_size});
particles.getDataSet("py").resize({new_size});
particles.getDataSet("pz").resize({new_size});
particles.getDataSet("e").resize({new_size});
particles.getDataSet("m").resize({new_size});
particles.getDataSet("lifetime").resize({new_size});
particles.getDataSet("spin").resize({new_size});
}
void write(Event const & ev){
cache.fill(ev);
if(cache.is_full()) {
dump_cache();
}
const double wt = ev.central().weight;
const size_t idx = to_index(ev.type());
heprup.XSECUP[idx] += wt;
heprup.XERRUP[idx] += wt*wt;
if(wt > heprup.XMAXUP[idx]){
heprup.XMAXUP[idx] = wt;
}
}
void dump_cache() {
write_event_params();
write_event_particles();
cache.clear();
}
void write_event_params() {
auto events = file.getGroup("event");
// choose arbitrary dataset to find size
const auto dataset = events.getDataSet("nparticles");
const size_t size = dataset.getSpace().getDimensions().front();
resize_event_group(size + cache.nparticles.size());
#define WRITE_FROM_CACHE(GROUP, PROPERTY) \
GROUP.getDataSet(#PROPERTY).select({size}, {cache.PROPERTY.size()}).write(cache.PROPERTY)
WRITE_FROM_CACHE(events, nparticles);
WRITE_FROM_CACHE(events, start);
WRITE_FROM_CACHE(events, pid);
WRITE_FROM_CACHE(events, weight);
WRITE_FROM_CACHE(events, scale);
WRITE_FROM_CACHE(events, fscale);
WRITE_FROM_CACHE(events, rscale);
WRITE_FROM_CACHE(events, aqed);
WRITE_FROM_CACHE(events, aqcd);
WRITE_FROM_CACHE(events, trials);
WRITE_FROM_CACHE(events, npLO);
WRITE_FROM_CACHE(events, npNLO);
}
void write_event_particles() {
auto particles = file.getGroup("particle");
// choose arbitrary dataset to find size
const auto dataset = particles.getDataSet("id");
const size_t size = dataset.getSpace().getDimensions().front();
resize_particle_group(size + cache.id.size());
WRITE_FROM_CACHE(particles, id);
WRITE_FROM_CACHE(particles, status);
WRITE_FROM_CACHE(particles, lifetime);
WRITE_FROM_CACHE(particles, spin);
WRITE_FROM_CACHE(particles, mother1);
WRITE_FROM_CACHE(particles, mother2);
WRITE_FROM_CACHE(particles, color1);
WRITE_FROM_CACHE(particles, color2);
WRITE_FROM_CACHE(particles, px);
WRITE_FROM_CACHE(particles, py);
WRITE_FROM_CACHE(particles, pz);
WRITE_FROM_CACHE(particles, e);
WRITE_FROM_CACHE(particles, m);
}
#undef WRITE_FROM_CACHE
~HDF5WriterImpl(){
dump_cache();
auto proc_info = file.getGroup("procInfo");
write_dataset(proc_info, "xSection", heprup.XSECUP);
write_dataset(proc_info, "error", heprup.XERRUP);
write_dataset(proc_info, "unitWeight", heprup.XMAXUP);
}
};
HDF5Writer::HDF5Writer(std::string const & filename, LHEF::HEPRUP heprup):
impl_{std::make_unique<HDF5Writer::HDF5WriterImpl>(
filename, std::move(heprup))}
{}
void HDF5Writer::write(Event const & ev){
impl_->write(ev);
}
}
#else // no HDF5 support
namespace HEJ{
class HDF5Writer::HDF5WriterImpl{};
HDF5Writer::HDF5Writer(std::string const &, LHEF::HEPRUP){
throw std::invalid_argument{
"Failed to create HDF5 writer: "
"HEJ 2 was built without HDF5 support"
};
}
void HDF5Writer::write(Event const &){
assert(false);
}
}
#endif
namespace HEJ {
HDF5Writer::~HDF5Writer() = default;
}
diff --git a/src/JetSplitter.cc b/src/JetSplitter.cc
index ade6731..e1ac9c2 100644
--- a/src/JetSplitter.cc
+++ b/src/JetSplitter.cc
@@ -1,192 +1,191 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/JetSplitter.hh"
#include <array>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <numeric>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/PseudoJet.hh"
#include "HEJ/Constants.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/RNG.hh"
namespace HEJ {
- namespace{
- constexpr double ccut=HEJ::CMINPT; // min parton pt
+ namespace {
template<class Iterator>
bool same_pt_and_rapidity(
Iterator begin, Iterator end,
fastjet::PseudoJet const & jet
){
constexpr double ep = 1e-2;
const fastjet::PseudoJet reconstructed_jet = std::accumulate(
begin, end, fastjet::PseudoJet{}
);
return
(std::abs(reconstructed_jet.pt() - jet.pt()) < ep)
&& (std::abs(reconstructed_jet.rapidity() - jet.rapidity()) < ep)
;
}
bool all_in_one_jet(
std::vector<fastjet::PseudoJet> const & partons,
fastjet::JetDefinition jet_def, double min_jet_pt
){
fastjet::ClusterSequence ev(partons, jet_def);
const std::vector<fastjet::PseudoJet> testjet = ev.inclusive_jets(min_jet_pt);
return testjet.size() == 1u
&& testjet[0].constituents().size() == partons.size();
}
}
JetSplitter::JetSplitter(
fastjet::JetDefinition jet_def, double min_pt
):
min_jet_pt_{min_pt},
jet_def_{std::move(jet_def)}
{}
using SplitResult = JetSplitter::SplitResult;
SplitResult JetSplitter::split(
fastjet::PseudoJet const & j2split, int ncons, RNG & ran
) const{
if(ncons <= 0) {
throw std::invalid_argument{
"number of requested jet constituents less than 1"
};
}
double swt = 1.;
std::vector<fastjet::PseudoJet> jcons;
if(ncons == 1){
jcons.emplace_back(j2split);
jcons.back().set_user_index(0);
return {jcons, swt};
}
if(ncons == 2){
return Split2(j2split, ran);
}
const double R_max = R_factor*jet_def_.R();
assert(R_max < M_PI);
double pt_remaining = j2split.pt();
const double phi_jet = j2split.phi();
const double y_jet = j2split.rapidity();
for(int i = 0; i < ncons - 1; ++i){
/**
* Generate rapidity and azimuthal angle with a distance
* R = sqrt(delta_y^2 + delta_phi^2) < R_max
* from the jet centre
*/
const double R = R_max*ran.flat();
const double theta = 2*M_PI*ran.flat();
const double delta_phi = R*std::cos(theta);
const double delta_y = R*std::sin(theta);
/**
* Generate pt such that the total contribution of all partons
* along the jet pt axis does not exceed the jet pt
*/
const double pt_max = pt_remaining/std::cos(delta_phi);
assert(pt_max > 0);
- if(pt_max < ccut) return {}; // no pt remaining for this parton
- const double pt = (pt_max - ccut)*ran.flat() + ccut;
+ if(pt_max < CMINPT) return {}; // no pt remaining for this parton
+ const double pt = (pt_max - CMINPT)*ran.flat() + CMINPT;
pt_remaining -= pt*std::cos(delta_phi);
jcons.emplace_back(
pt*std::cos(phi_jet + delta_phi), pt*std::sin(phi_jet + delta_phi),
pt*std::sinh(y_jet + delta_y), pt*std::cosh(y_jet + delta_y)
);
jcons.back().set_user_index(i);
- swt *= 2*M_PI*R*R_max*pt*(pt_max - ccut);
+ swt *= 2*M_PI*R*R_max*pt*(pt_max - CMINPT);
}
const fastjet::PseudoJet p_total = std::accumulate(
jcons.cbegin(), jcons.cend(), fastjet::PseudoJet{}
);
// Calculate the pt of the last parton
const double last_px = j2split.px() - p_total.px();
const double last_py = j2split.py() - p_total.py();
const double last_pt = std::sqrt(last_px*last_px + last_py*last_py);
- if(last_pt < ccut) return {};
+ if(last_pt < CMINPT) return {};
// Calculate the rapidity of the last parton using the requirement that the
// new jet must have the same rapidity as the LO jet.
const double exp_2y_jet = (j2split.e() + j2split.pz())/(j2split.e() - j2split.pz());
const double bb = (p_total.e()+p_total.pz()) - exp_2y_jet*(p_total.e()-p_total.pz());
const double lasty = std::log(
( -bb+std::sqrt(bb*bb+4.*exp_2y_jet*last_pt*last_pt) ) / ( 2.*last_pt ));
jcons.emplace_back(
last_px, last_py, last_pt*std::sinh(lasty), last_pt*std::cosh(lasty)
);
jcons.back().set_user_index(ncons-1);
assert(same_pt_and_rapidity(jcons.cbegin(), jcons.cend(), j2split));
// Test that the last parton is not too far away from the jet centre.
if (jcons.back().delta_R(j2split) > R_max) return {};
if(! all_in_one_jet(jcons, jet_def_, min_jet_pt_)) return {};
return {jcons, swt};
}
double JetSplitter::sample_distance_2p(double & wt, RNG & ran) const{
static constexpr double x_small = 0.1;
static constexpr double p_small = 0.4;
const double pR = ran.flat();
if(pR < p_small){
wt *= x_small/p_small;
return x_small/p_small*pR;
}
wt *= (1-x_small)/(1-p_small);
return (1-x_small)/(1-p_small)*(pR-p_small) + x_small;
}
SplitResult JetSplitter::Split2(
fastjet::PseudoJet const & j2split, RNG & ran
) const {
static constexpr std::size_t ncons = 2;
std::vector<fastjet::PseudoJet> jcons(ncons);
std::array<double, ncons> R, phi, y, pt;
double wt = 1;
const double theta = 2*M_PI*ran.flat(); // angle in y-phi plane
// empiric observation: we are always within the jet radius
R[0] = sample_distance_2p(wt, ran)*jet_def_.R();
R[1] = -sample_distance_2p(wt, ran)*jet_def_.R();
for(std::size_t i = 0; i <= 1; ++i){
phi[i] = j2split.phi() + R[i]*std::cos(theta);
y[i] = j2split.rapidity() + R[i]*std::sin(theta);
}
for(std::size_t i = 0; i <= 1; ++i){
pt[i] = (j2split.py() - std::tan(phi[1-i])*j2split.px())/
(std::sin(phi[i]) - std::tan(phi[1-i])*std::cos(phi[i]));
- if(pt[i] < ccut) return {};
+ if(pt[i] < CMINPT) return {};
jcons[i].reset_PtYPhiM(pt[i], y[i], phi[i]);
jcons[i].set_user_index(i);
}
assert(same_pt_and_rapidity(jcons.cbegin(), jcons.cend(), j2split));
if(! all_in_one_jet(jcons, jet_def_, min_jet_pt_)) return {};
wt *= 2*M_PI*pt[0]*R[0]*jet_def_.R()*jet_def_.R();
// from transformation of delta(R[1] - ...) to delta(pt[0] - ...)
const double dphi0 = phi[0] - j2split.phi();
const double ptJ = j2split.pt();
const double jacobian = std::cos(theta)*pt[1]*pt[1]/(ptJ*std::sin(dphi0));
return {jcons, jacobian*wt};
}
}
diff --git a/src/LesHouchesWriter.cc b/src/LesHouchesWriter.cc
index a8b8726..84f1dff 100644
--- a/src/LesHouchesWriter.cc
+++ b/src/LesHouchesWriter.cc
@@ -1,125 +1,125 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/LesHouchesWriter.hh"
#include <cassert>
#include <cmath>
#include <cstddef>
#include <utility>
#include <vector>
#include "HEJ/Event.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Parameters.hh"
-namespace HEJ{
+namespace HEJ {
namespace {
std::size_t to_index(event_type::EventType const type){
return type==0?0:1+std::floor(std::log2(static_cast<std::size_t>(type)));
}
}
LesHouchesWriter::LesHouchesWriter(
std::string const & file, LHEF::HEPRUP heprup
):
out_{file, std::fstream::in | std::fstream::out | std::fstream::trunc},
writer_{std::make_unique<LHEF::Writer>(out_)}
{
if(! out_.is_open()){
throw std::ios_base::failure("Failed to open " + file);
};
// scientific style is needed to allow rewriting the init block
out_ << std::scientific;
writer_->heprup = std::move(heprup);
// lhe Standard: IDWTUP (negative => weights = +/-)
// IDWTUP: HEJ -> SHG/Pythia/next program
// 1: weighted->unweighted, xs = mean(weight), XMAXUP given
// 2: weighted->unweighted, xs = XSECUP, XMAXUP given
// 3: unweighted (weight=+1)->unweighted, no additional information
// 4: weighted->weighted, xs = mean(weight)
//
// None of these codes actually match what we want:
// 1 and 4 require xs = mean(weight), which is impossible until after generation
// 2 tells the SHG to unweight our events, which is wasteful
// 3 claims we produce unweighted events, which is both wasteful _and_
// impossible until after generation (we don't know the maximum weight before)
//
// For the time being, we choose -3. If the consumer (like Pythia) assumes
// weight=+1, the final weights have to be corrected by multiplying with
// the original weight we provided. We are also often use NLO-PDFs which can
// give negative weights, hence the native IDWTUP.
//
writer_->heprup.IDWTUP = -3;
// always use the newest LHE version
// Pythia only saves weights (hepeup.XWGTUP) for version >=2
writer_->heprup.version = LHEF::HEPRUP().version;
const int max_number_types = to_index(event_type::last_type)+1;
writer_->heprup.NPRUP = max_number_types;
// ids of event types
writer_->heprup.LPRUP.clear();
writer_->heprup.LPRUP.reserve(max_number_types);
writer_->heprup.LPRUP.emplace_back(0);
for(auto i=event_type::first_type+1; i<=event_type::last_type; i*=2)
writer_->heprup.LPRUP.emplace_back(i);
// use placeholders for unknown init block values
// we can overwrite them after processing all events
writer_->heprup.XSECUP = std::vector<double>(max_number_types, 0.);
writer_->heprup.XERRUP = std::vector<double>(max_number_types, 0.);
writer_->heprup.XMAXUP = std::vector<double>(max_number_types, 0.);
write_init();
}
void LesHouchesWriter::write(Event const & ev){
assert(writer_ && out_.is_open());
const double wt = ev.central().weight;
- writer_->hepeup = HEJ::to_HEPEUP(std::move(ev), &heprup());
+ writer_->hepeup = to_HEPEUP(std::move(ev), &heprup());
writer_->writeEvent();
assert(heprup().XSECUP.size() > to_index(ev.type()));
heprup().XSECUP[to_index(ev.type())] += wt;
heprup().XERRUP[to_index(ev.type())] += wt*wt;
if(wt > heprup().XMAXUP[to_index(ev.type())]){
heprup().XMAXUP[to_index(ev.type())] = wt;
}
}
// this function is called after overwritting the Les Houches init block
// assert that we have overwritten *exactly* the init block,
// i.e. we are at the end of the file or an intact event block is next
void assert_next_event_intact(std::istream & out){
(void) out; // suppress compiler warnings if not in debug mode
#ifndef NDEBUG
std::string line;
getline(out, line);
assert(out.eof() || line.rfind("<event", 0) == 0);
#endif
}
void LesHouchesWriter::rewrite_init(){
assert(writer_ && out_.is_open());
// replace placeholder entries
const auto pos = out_.tellp();
out_.seekp(0);
write_init();
assert_next_event_intact(out_);
out_.seekp(pos);
}
LesHouchesWriter::~LesHouchesWriter(){
assert(writer_ && out_.is_open());
for(auto & xs_err: heprup().XERRUP)
{
xs_err = std::sqrt(xs_err);
}
rewrite_init();
}
-}
+} // namespace HEJ
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 9144d95..245b1b0 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,2217 +1,2217 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <iterator>
#include <limits>
#include <unordered_map>
#include <utility>
#include "CLHEP/Vector/LorentzVector.h"
#include "fastjet/PseudoJet.hh"
#include "HEJ/ConfigFlags.hh"
#include "HEJ/Constants.hh"
#include "HEJ/Event.hh"
#include "HEJ/event_types.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/Hjets.hh"
#include "HEJ/jets.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/utility.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Zjets.hh"
-namespace HEJ{
+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*std::log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
return (
1. + alpha_s/(4.*M_PI)*beta0*std::log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(Event const & event) const {
std::vector <double> tree_kin_part=tree_kin(event);
std::vector <Weights> virtual_part=virtual_corrections(event);
if(tree_kin_part.size() != virtual_part.size()) {
throw std::logic_error("tree and virtuals have different sizes");
}
Weights sum = Weights{0., std::vector<double>(event.variations().size(), 0.)};
for(size_t i=0; i<tree_kin_part.size(); ++i) {
sum += tree_kin_part.at(i)*virtual_part.at(i);
}
return tree_param(event)*sum;
}
Weights MatrixElement::tree(Event const & event) const {
std::vector <double> tree_kin_part=tree_kin(event);
double sum = 0.;
for(size_t i=0; i<tree_kin_part.size(); ++i) {
sum += tree_kin_part.at(i);
}
return tree_param(event)*sum;
}
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;
}
std::vector<Weights> MatrixElement::virtual_corrections(Event const & event) const {
if(! is_resummable(event.type())) {
return {Weights{0., std::vector<double>(event.variations().size(), 0.)}};
}
// only compute once for each renormalisation scale
std::unordered_map<double, std::vector<double> > known_vec;
std::vector<double> central_vec=virtual_corrections(event, event.central().mur);
known_vec.emplace(event.central().mur, central_vec);
for(auto const & var: event.variations()) {
const auto ME_it = known_vec.find(var.mur);
if(ME_it == end(known_vec)) {
known_vec.emplace(var.mur, virtual_corrections(event, var.mur));
}
}
// At this stage known_vec contains one vector of virtual corrections for each mur value
// Now put this into a vector of Weights
std::vector<Weights> result_vec;
for(size_t i=0; i<central_vec.size(); ++i) {
Weights result;
result.central = central_vec.at(i);
for(auto const & var: event.variations()) {
const auto ME_it = known_vec.find(var.mur);
result.variations.emplace_back(ME_it->second.at(i));
}
result_vec.emplace_back(result);
}
return result_vec;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
const 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;
std::size_t first_idx = 0;
std::size_t last_idx = partons.size() - 1;
#ifndef NDEBUG
bool wc = true;
#endif
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::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else if (event.type() == event_type::qqxexb) {
q -= partons[1].p;
++first_idx;
if (std::abs(partons[0].type) != std::abs(partons[1].type)){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else {
if(event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
if (in[0].type != partons[0].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
std::size_t first_idx_qqx = last_idx;
std::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(std::abs(backquark->type) != std::abs((backquark+1)->type)) {
wqq=true;
#ifndef NDEBUG
wc=false;
#endif
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(std::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(std::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;
}
#ifndef NDEBUG
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
);
#endif
return std::exp(exponent);
}
std::vector <double> MatrixElement::virtual_corrections_Z_qq(
Event const & event,
const double mur,
Particle const & ZBoson
) 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;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q_t = pa - partons[0].p - ZBoson.p;
fastjet::PseudoJet q_b = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q_t -= partons[1].p;
q_b -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--last_idx;
}
double sum_top=0., sum_bot=0., sum_mix=0.;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
const double dy = partons[j+1].rapidity() - partons[j].rapidity();
const double tmp_top = omega0(alpha_s, mur, q_t)*dy;
const double tmp_bot = omega0(alpha_s, mur, q_b)*dy;
sum_top += tmp_top;
sum_bot += tmp_bot;
sum_mix += (tmp_top + tmp_bot) / 2.;
q_t -= partons[j+1].p;
q_b -= partons[j+1].p;
}
return {exp(sum_top), exp(sum_bot), exp(sum_mix)};
}
double MatrixElement::virtual_corrections_Z_qg(
Event const & event,
const double mur,
Particle const & ZBoson,
const bool is_gq_event
) 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;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
// If this is a gq event, don't subtract the Z momentum from first q
fastjet::PseudoJet q = (is_gq_event ? pa - partons[0].p : pa - partons[0].p - ZBoson.p);
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--last_idx;
}
double sum=0.;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
sum += omega0(alpha_s, mur, q)*(partons[j+1].rapidity()
- partons[j].rapidity());
q -= partons[j+1].p;
}
return exp(sum);
}
std::vector<double> MatrixElement::virtual_corrections(
Event const & event,
const 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) && std::abs(AWZH_boson->type) == pid::Wp){
return {virtual_corrections_W(event, mur, *AWZH_boson)};
}
if(AWZH_boson != end(out) && AWZH_boson->type == pid::Z_photon_mix){
if(is_gluon(in.back().type)){
// This is a qg event
return {virtual_corrections_Z_qg(event, mur, *AWZH_boson, false)};
}else if(is_gluon(in.front().type)){
// This is a gq event
return {virtual_corrections_Z_qg(event, mur, *AWZH_boson, true)};
}else{
// This is a qq event
return virtual_corrections_Z_qq(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;
std::size_t first_idx = 0;
std::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;
}
std::size_t first_idx_qqx = last_idx;
std::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(std::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(std::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 {std::exp(exponent)};
}
namespace {
//! Lipatov vertex for partons emitted into extremal jets
CLHEP::HepLorentzVector CLipatov(
CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
const CLHEP::HepLorentzVector p5 = qav-qbv;
const CLHEP::HepLorentzVector CL = -(qav+qbv)
+ 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;
}
double C2Lipatov(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2
){
const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, p1, p2);
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
const double kperp=(qav-qbv).perp();
if (kperp>lambda)
return Cls;
return Cls-4./(kperp*kperp);
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, p1, p2);
const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, p1, p2);
return -CL_t.dot(CL_b);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, p1, p2)
/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
const double kperp = (qav_t - qbv_t).perp();
if (kperp > lambda) {
return Cls;
} else {
return Cls - 4.0 / (kperp * kperp);
}
}
CLHEP::HepLorentzVector CLipatov(
CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
){
const CLHEP::HepLorentzVector p5 = qav-qbv;
const CLHEP::HepLorentzVector CL = -(qav+qbv)
+ 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;
}
//! Lipatov vertex
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
){
const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, pim, pip, pom, pop);
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
const double kperp=(qav-qbv).perp();
if (kperp>lambda)
return Cls;
return Cls-4./(kperp*kperp);
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
) {
const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, pim, pip, pom, pop);
const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, pim, pip, pom, pop);
return -CL_t.dot(CL_b);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, pa, pb, p1, p2)
/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
const double kperp = (qav_t - qbv_t).perp();
if (kperp > lambda) {
return Cls;
} else {
return Cls - 4.0 / (kperp * kperp);
}
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pg Unordered gluon momentum
* @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
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_uno_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
using namespace currents;
assert(aptype!=pid::gluon); // aptype cannot be gluon
if (bptype==pid::gluon) {
if (is_quark(aptype))
return ME_unob_qg(pg,p1,pa,pn,pb);
else
return ME_unob_qbarg(pg,p1,pa,pn,pb);
}
else if (is_antiquark(bptype)) {
if (is_quark(aptype))
return ME_unob_qQbar(pg,p1,pa,pn,pb);
else
return ME_unob_qbarQbar(pg,p1,pa,pn,pb);
}
else { //bptype == quark
if (is_quark(aptype))
return ME_unob_qQ(pg,p1,pa,pn,pb);
else
return ME_unob_qbarQ(pg,p1,pa,pn,pb);
}
throw std::logic_error("unreachable");
}
/** Matrix element squared for tree-level current-current scattering
* @param bptype Particle b PDG ID
* @param pgin Incoming gluon momentum
* @param pq Quark from splitting Momentum
* @param pqbar Anti-quark from splitting Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param swap_q_qx Boolean. Ordering of qqbar pair. False: pqbar extremal.
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The qqxf contribution can be calculated by reversing the argument ordering.
*/
double ME_qqx_current(
ParticleID bptype,
CLHEP::HepLorentzVector const & pgin,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
bool const swap_q_qx
){
using namespace currents;
if (bptype==pid::gluon) {
if (swap_q_qx) // pq extremal
return ME_Exqqx_qqbarg(pgin,pq,pqbar,pn,pb);
else // pqbar extremal
return ME_Exqqx_qbarqg(pgin,pq,pqbar,pn,pb);
}
else { // b leg quark line
if (swap_q_qx) //extremal pq
return ME_Exqqx_qqbarQ(pgin,pq,pqbar,pn,pb);
else
return ME_Exqqx_qbarqQ(pgin,pq,pqbar,pn,pb);
}
throw std::logic_error("unreachable");
}
/* \brief Matrix element squared for central qqx tree-level current-current
* scattering
*
* @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
* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
*/
double ME_qqxmid_current(
ParticleID aptype, ParticleID bptype, int nabove,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<CLHEP::HepLorentzVector> const & partons
){
using namespace currents;
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
const bool swap_q_qx=pqbar.rapidity() < pq.rapidity();
double wt=1.;
if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F;
return wt*ME_Cenqqx_qq(pa, pb, partons, is_antiquark(bptype),
is_antiquark(aptype), swap_q_qx, nabove);
}
/** 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(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
using namespace currents;
if (aptype==pid::gluon && bptype==pid::gluon) {
return ME_gg(pn,pb,p1,pa);
} else if (aptype==pid::gluon && bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_qg(pn,pb,p1,pa);
else
return ME_qbarg(pn,pb,p1,pa);
}
else if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_qg(p1,pa,pn,pb);
else
return ME_qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_qQ(pn,pb,p1,pa);
else
return ME_qQbar(pn,pb,p1,pa);
}
else {
if (is_quark(aptype))
return ME_qQbar(p1,pa,pn,pb);
else
return ME_qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unreachable");
}
/** 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(
ParticleID aptype, ParticleID 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, ParticleProperties const & Wprop
){
using namespace currents;
// We know it cannot be gg incoming.
assert(!(aptype==pid::gluon && bptype==pid::gluon));
if (aptype==pid::gluon && bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_W_qg(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qbarg(pn,plbar,pl,pb,p1,pa,Wprop);
}
else if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_W_qg(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qbarg(p1,plbar,pl,pa,pn,pb,Wprop);
}
else { // they are both quark
if (wc==true){ // emission off b, (first argument pbout)
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_W_qQ(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
else {
if (is_quark(aptype))
return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
}
else{ // emission off a, (first argument paout)
if (is_quark(aptype)) {
if (is_quark(bptype))
return ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
else { // a is anti-quark
if (is_quark(bptype))
return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
}
}
throw std::logic_error("unreachable");
}
/** 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
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_W_uno_current(
ParticleID aptype, ParticleID 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, ParticleProperties const & Wprop
){
using namespace currents;
// we know they are not both gluons
assert(bptype != pid::gluon || aptype != pid::gluon);
if (bptype == pid::gluon && aptype != pid::gluon) {
// b gluon => W emission off a
if (is_quark(aptype))
return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else { // they are both quark
if (wc) {// emission off b, i.e. b is first current
if (is_quark(bptype)){
if (is_quark(aptype))
return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else{
if (is_quark(aptype))
return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (is_quark(aptype)) {
if (is_quark(bptype)) //qq
return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else //qqbar
return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else { // a is anti-quark
if (is_quark(bptype)) //qbarq
return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else //qbarqbar
return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
}
}
throw std::logic_error("unreachable");
}
/** \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 swap_q_qx Boolean. Ordering of qqbar pair. False: pqbar extremal.
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxb Tree-Level Current-Current Scattering
*
* @note calculate forwards qqx contribution by reversing argument ordering.
*/
double ME_W_qqx_current(
ParticleID aptype, ParticleID 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 swap_q_qx, bool const wc,
ParticleProperties const & Wprop
){
using namespace currents;
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
- const double CFbackward = K_g( (swap_q_qx)?pq:pqbar ,pa)/HEJ::C_F;
+ const double CFbackward = K_g( (swap_q_qx)?pq:pqbar ,pa)/C_F;
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==pid::gluon && bptype==pid::gluon) {
//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission
// Site.
if (swap_q_qx)
return ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb, Wprop)
* CFbackward;
else
return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)
* CFbackward;
}
else {
assert(aptype==pid::gluon && bptype!=pid::gluon );
//a gluon => W emission off b leg or qqx
if (!wc){ // W Emitted from backwards qqx
if (swap_q_qx)
return ME_WExqqx_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop)
* CFbackward;
else
return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)
* CFbackward;
}
else { // W Must be emitted from forwards leg.
if (swap_q_qx)
return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, is_antiquark(bptype), Wprop)
* CFbackward;
else
return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop)
* CFbackward;
}
}
throw std::logic_error("unreachable");
}
/* \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(
ParticleID aptype, ParticleID bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<CLHEP::HepLorentzVector> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc,
ParticleProperties const & Wprop
){
using namespace currents;
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
const bool swap_q_qx=pqbar.rapidity() < pq.rapidity();
double wt=1.;
- if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/HEJ::C_F;
- if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/HEJ::C_F;
+ if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/C_F;
+ if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F;
if(wqq)
return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons,
is_antiquark(bptype),is_antiquark(aptype),
swap_q_qx, nabove, Wprop);
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons,
is_antiquark(bptype), is_antiquark(aptype),
swap_q_qx, nabove, nbelow, wc, Wprop);
}
/** Matrix element squared for tree-level current-current scattering With Z+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 plbar Final state positron momentum
* @param pl Final state electron momentum
* @param Zprop Z properties
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
std::vector<double> ME_Z_current(
const ParticleID aptype, const ParticleID 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,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
if(is_anyquark(aptype) && is_gluon(bptype)){
// This is a qg event
const double current_factor=ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
return { current_factor };
}else if(is_gluon(aptype) && is_anyquark(bptype)){
// This is a gq event
const double current_factor=ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw);
return { current_factor };
}else if(is_anyquark(aptype) && is_anyquark(bptype)){
// This is a qq event
return ME_Z_qQ(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
}else{
// It cannot be gg incoming
throw std::logic_error("Bad particle types");
}
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With Z+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 plbar Final state positron momentum
* @param pl Final state electron momentum
* @param Zprop Z properties
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
std::vector<double> ME_Z_uno_current(
const ParticleID aptype, const ParticleID 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,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
if (is_anyquark(aptype) && is_gluon(bptype)) {
// This is a qg event
const double current_factor=ME_Zuno_qg(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
return { current_factor };
} else if (is_gluon(aptype) && is_anyquark(bptype)) {
// This is a gq event
const double current_factor=ME_Zuno_qg(pb,pa,pg,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw);
return { current_factor };
} else if (is_anyquark(aptype) && is_anyquark(bptype)) {
// This is a qq event
return ME_Zuno_qQ(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
}
throw std::logic_error("unreachable");
}
/** \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(
ParticleID aptype, ParticleID 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, double vev
){
using namespace currents;
if (aptype==pid::gluon && bptype==pid::gluon)
// gg initial state
return ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev);
else if (aptype==pid::gluon&&bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_H_qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
else
return ME_H_qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
}
else if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_H_qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
else
return ME_H_qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
}
else { // they are both quark
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
else
return ME_H_qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
else {
if (is_quark(aptype))
return ME_H_qbarQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
else
return ME_H_qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
}
throw std::logic_error("unreachable");
}
/** \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 pg 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
*
* @note This function assumes unordered gluon backwards from pa-p1 current.
* For unof, reverse call order
*/
double ME_Higgs_current_uno(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pg,
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, double vev
){
using namespace currents;
if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_H_unob_gQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_gQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
else { // they are both quark
if (is_quark(aptype)) {
if (is_quark(bptype))
return ME_H_unob_qQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_qbarQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
else {
if (is_quark(bptype))
return ME_H_unob_qQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_qbarQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
}
throw std::logic_error("unreachable");
}
- CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){
+ CLHEP::HepLorentzVector to_HepLorentzVector(Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
- void validate(HEJ::MatrixElementConfig const & config) {
+ void validate(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
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
std::vector<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)};
case pid::Z_photon_mix:
return tree_kin_Z(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;
}
template<class InputIterator>
std::vector <double> FKL_ladder_weight_mix(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0_t, CLHEP::HepLorentzVector const & q0_b,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
const double lambda
){
double wt_top = 1;
double wt_bot = 1;
double wt_mix = 1;
auto qi_t = q0_t;
auto qi_b = q0_b;
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_t = qi_t - g;
const auto qip1_b = qi_b - g;
if(treat_as_extremal(*gluon_it)){
wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, lambda)*C_A;
wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, lambda)*C_A;
wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, lambda)*C_A;
} else{
wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, p1, pn, lambda)*C_A;
wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A;
wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A;
}
qi_t = qip1_t;
qi_b = qip1_b;
}
return {wt_top, wt_bot, wt_mix};
}
} // 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;
}
const auto & jets = ev.jets();
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 = ev.particle_jet_indices(
{*most_backward, *most_forward}
);
assert(extremal_jet_indices.size() == out_partons.size());
for(std::size_t i = 0; i < out_partons.size(); ++i){
- assert(HEJ::is_parton(out_partons[i]));
+ assert(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;
}
namespace {
double tree_kin_jets_qqxmid(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,
std::vector<Particle> const & partons,
double lambda
){
CLHEP::HepLorentzVector 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;
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);
}
const int nabove = std::distance(begin_ladder, backmidquark);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for (std::size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_qqxmid_current(
aptype, bptype, nabove, pa, pb,
pq, pqbar, partonsHLV
);
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*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_qqx(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda){
const bool swap_q_qx = is_quark(*BeginPart);
const auto pgin = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pq = to_HepLorentzVector(*(BeginPart+(swap_q_qx?0:1)));
const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_q_qx?1:0)));
const auto p1 = to_HepLorentzVector(*(BeginPart));
const auto pn = to_HepLorentzVector(*(EndPart-1));
assert((BeginIn)->type==pid::gluon); // Incoming a must be gluon.
const double current_factor = ME_qqx_current(
(EndIn-1)->type, pgin, pq, pqbar, pn, pb, swap_q_qx
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pgin-pq-pqbar, pgin, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const double current_factor = ME_uno_current(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pa-p1-pg, pa, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
}
double MatrixElement::tree_kin_jets(Event const & ev) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
- if (ev.type()==HEJ::event_type::FKL){
+ if (ev.type()==event_type::FKL){
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
);
}
- else if (ev.type()==HEJ::event_type::unordered_backward){
+ else if (ev.type()==event_type::unordered_backward){
return tree_kin_jets_uno(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
- else if (ev.type()==HEJ::event_type::unordered_forward){
+ else if (ev.type()==event_type::unordered_forward){
return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
- else if (ev.type()==HEJ::event_type::extremal_qqxb){
+ else if (ev.type()==event_type::extremal_qqxb){
return tree_kin_jets_qqx(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
- else if (ev.type()==HEJ::event_type::extremal_qqxf){
+ else if (ev.type()==event_type::extremal_qqxf){
return tree_kin_jets_qqx(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
- else if (ev.type()==HEJ::event_type::central_qqx){
+ else if (ev.type()==event_type::central_qqx){
return tree_kin_jets_qqxmid(incoming[0].type, incoming[1].type,
to_HepLorentzVector(incoming[0]),
to_HepLorentzVector(incoming[1]),
partons, param_.regulator_lambda);
}
else {
throw std::logic_error("Cannot reweight non-resummable processes in Pure Jets");
}
}
namespace{
double tree_kin_W_FKL(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector plbar, CLHEP::HepLorentzVector pl,
double lambda, ParticleProperties const & Wprop
){
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_uno(InIter BeginIn, partIter BeginPart,
partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
double lambda, ParticleProperties const & Wprop
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
bool wc = (BeginIn)->type==(BeginPart+1)->type; //leg b emits w
auto q0 = pa - p1 - pg;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_uno_current(
(BeginIn)->type, (BeginIn+1)->type, pn, pb,
p1, pa, pg, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_qqx(InIter BeginIn, partIter BeginPart,
partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
double lambda, ParticleProperties const & Wprop
){
const bool swap_q_qx=is_quark(*BeginPart);
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pq = to_HepLorentzVector(*(BeginPart+(swap_q_qx?0:1)));
const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_q_qx?1:0)));
const auto p1 = to_HepLorentzVector(*(BeginPart));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const bool wc = (BeginIn+1)->type!=(EndPart-1)->type; //leg b emits w
auto q0 = pa - pq - pqbar;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_qqx_current(
(BeginIn)->type, (BeginIn+1)->type, pa, pb,
pq, pqbar, pn, plbar, pl, swap_q_qx, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
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(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector pa,
CLHEP::HepLorentzVector pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector plbar, CLHEP::HepLorentzVector pl,
double lambda, ParticleProperties const & Wprop
){
CLHEP::HepLorentzVector 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.front());
auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqx emit W
bool wc = !wqq && (aptype==partons.front().type); //leg b emits w
assert(!wqq || !wc);
if(wqq){ // emission from qqx
qqxt -= pl + plbar;
} else if(!wc) { // emission from leg a
q0 -= pl + plbar;
qqxt -= pl + plbar;
}
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);
}
const int nabove = std::distance(begin_ladder, backmidquark);
const int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for (std::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, Wprop
);
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());
#ifndef NDEBUG
// assert that there is exactly one decay corresponding to the W
assert(ev.decays().size() == 1);
auto const & w_boson{
std::find_if(ev.outgoing().cbegin(), ev.outgoing().cend(),
[] (Particle const & p) -> bool {
return std::abs(p.type) == ParticleID::Wp;
}) };
assert(w_boson != ev.outgoing().cend());
assert( static_cast<long int>(ev.decays().cbegin()->first)
== std::distance(ev.outgoing().cbegin(), w_boson) );
#endif
// find decay products of W
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert( ( is_anylepton(decay.at(0)) && is_anyneutrino(decay.at(1)) )
|| ( is_anylepton(decay.at(1)) && is_anyneutrino(decay.at(0)) ) );
// get lepton & neutrino
CLHEP::HepLorentzVector plbar, pl;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
}
else{
pl = to_HepLorentzVector(decay.at(0));
plbar = to_HepLorentzVector(decay.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() == FKL){
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_backward){
return tree_kin_W_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_forward){
return tree_kin_W_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqxb){
return tree_kin_W_qqx(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqxf){
return tree_kin_W_qqx(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
assert(ev.type() == central_qqx);
return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
namespace{
std::vector <double> tree_kin_Z_FKL(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw
){
const auto p1 = to_HepLorentzVector(partons[0]);
const auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
const std::vector <double> current_factor = ME_Z_current(
aptype, bptype, pn, pb, p1, pa,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
q0, pa, pb, p1, pn, lambda));
}else if(is_gluon(aptype)){
// This is a gq event
const auto q0 = pa-p1;
ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a qq event
const auto q0 = pa-p1-plbar-pl;
const auto q1 = pa-p1;
ladder_factor=FKL_ladder_weight_mix(begin_ladder, end_ladder,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); i++){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
template<class InIter, class partIter>
std::vector <double> tree_kin_Z_uno(InIter BeginIn, partIter BeginPart, partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const ParticleID aptype = (BeginIn)->type;
const ParticleID bptype = (BeginIn+1)->type;
const std::vector <double> current_factor = ME_Z_uno_current(
aptype, bptype, pn, pb, p1, pa, pg,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-pg-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else if(is_gluon(aptype)){
// This is a gq event
const auto q0 = pa-pg-p1;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a qq event
const auto q0 = pa-pg-p1-plbar-pl;
const auto q1 = pa-pg-p1;
ladder_factor=FKL_ladder_weight_mix(BeginPart+2, EndPart-1,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); i++){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
} // namespace anonymous
std::vector<double> MatrixElement::tree_kin_Z(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
// find decay products of Z
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert(is_anylepton(decay.at(0)) && !is_anyneutrino(decay.at(0))
&& decay.at(0).type==-decay.at(1).type);
// get leptons
CLHEP::HepLorentzVector plbar, pl;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
}
else{
pl = to_HepLorentzVector(decay.at(0));
plbar = to_HepLorentzVector(decay.at(1));
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
const double stw2 = param_.ew_parameters.sin2_tw();
const double ctw = param_.ew_parameters.cos_tw();
if(ev.type() == FKL){
return tree_kin_Z_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_backward){
return tree_kin_Z_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_forward){
return tree_kin_Z_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
throw std::logic_error("Can only reweight FKL or uno processes in Z+Jets");
}
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
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(type == pid::gluon) return currents::K_g(pout, pin);
return C_F;
}
#endif
// Colour factor in strict MRK limit
double K_MRK(ParticleID type) {
return (type == pid::gluon)?C_A:C_F;
}
}
double MatrixElement::MH2_forwardH(
CLHEP::HepLorentzVector const & p1out,
CLHEP::HepLorentzVector const & p1in,
ParticleID type2,
CLHEP::HepLorentzVector const & p2out,
CLHEP::HepLorentzVector const & p2in,
CLHEP::HepLorentzVector const & pH,
double t1, double t2
) const{
using namespace currents;
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
const double vev = param_.ew_parameters.vev();
// 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, vev
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH,vev) + C2gHgm(p1in,p1out,pH,vev)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(outgoing[1].type != pid::gluon) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
const auto q0 = pa - p1 - pH;
const double t1 = q0.m2();
const double t2 = (pn - pb).m2();
return MH2_forwardH(
p1, pa, incoming[1].type, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_last(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.back().type == pid::Higgs);
if(outgoing[outgoing.size()-2].type != pid::gluon) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
const double t1 = q0.m2();
const double t2 = (pn + pH - pb).m2();
return MH2_forwardH(
pn, pb, incoming[0].type, p1, pa, pH,
t2, t1
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
namespace {
template<class InIter, class partIter>
double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart,
CLHEP::HepLorentzVector const & qH,
CLHEP::HepLorentzVector const & qHp1,
double mt, bool inc_bot, double mb, double vev
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
return ME_Higgs_current_uno(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa,
qH, qHp1, mt, inc_bot, mb, vev
);
}
}
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() == FKL){
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,
param_.ew_parameters.vev()
);
}
else if(ev.type() == unob){
- current_factor = HEJ::C_A*HEJ::C_A/2*tree_kin_Higgs_uno(
+ current_factor = C_A*C_A/2*tree_kin_Higgs_uno(
begin(incoming), end(incoming), begin(partons),
end(partons), qH, qH-pH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
- current_factor = HEJ::C_A*HEJ::C_A/2*tree_kin_Higgs_uno(
+ current_factor = C_A*C_A/2*tree_kin_Higgs_uno(
rbegin(incoming), rend(incoming), rbegin(partons),
rend(partons), qH-pH, qH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
throw std::logic_error("Can only reweight FKL or uno processes in H+Jets");
}
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, double alpha_w) {
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 alpha_w*alpha_w;
case pid::Z_photon_mix:
return alpha_w*alpha_w;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
}
double MatrixElement::tree_param(Event const & ev, double mur) const {
assert(is_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*std::log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res;
}
} // namespace HEJ
diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc
index 004bba3..adf7649 100644
--- a/src/PhaseSpacePoint.cc
+++ b/src/PhaseSpacePoint.cc
@@ -1,849 +1,849 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/PhaseSpacePoint.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdlib>
#include <functional>
#include <iterator>
#include <limits>
#include <numeric>
#include <random>
#include <tuple>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "HEJ/Constants.hh"
#include "HEJ/Event.hh"
#include "HEJ/event_types.hh"
#include "HEJ/JetSplitter.hh"
#include "HEJ/kinematics.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/resummation_jet.hh"
#include "HEJ/RNG.hh"
#include "HEJ/utility.hh"
-namespace HEJ{
+namespace HEJ {
namespace {
constexpr int max_jet_user_idx = PhaseSpacePoint::ng_max;
bool is_nonjet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() > max_jet_user_idx;
}
bool is_jet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() <= max_jet_user_idx;
}
// user indices for partons with extremal rapidity
constexpr int qqxmid1_uid = -9;
constexpr int qqxmid2_uid = -8;
constexpr int qqxb_uid = -7;
constexpr int qqxf_uid = -6;
constexpr int unob_uid = -5;
constexpr int unof_uid = -4;
constexpr int backward_FKL_uid = -3;
constexpr int forward_FKL_uid = -2;
double estimate_ng_mean(std::vector<fastjet::PseudoJet> const & Born_jets){
const double delta_y =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
assert(delta_y > 0);
// Formula derived from fit in arXiv:1805.04446 (see Fig. 2)
return 0.975052*delta_y;
}
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
} // namespace anonymous
Event::EventData to_EventData(PhaseSpacePoint psp){
Event::EventData result;
result.incoming = std::move(psp).incoming_;
assert(result.incoming.size() == 2);
result.outgoing = std::move(psp).outgoing_;
// technically Event::EventData doesn't have to be sorted,
// but PhaseSpacePoint should be anyway
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
- HEJ::rapidity_less{}
+ rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.decays = std::move(psp).decays_;
result.parameters.central = {NaN, NaN, psp.weight()};
return result;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
return sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
}
bool PhaseSpacePoint::pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const{
return cluster_jets(jets).size() == jets.size();
}
int PhaseSpacePoint::sample_ng(
std::vector<fastjet::PseudoJet> const & Born_jets, RNG & ran
){
const double ng_mean = estimate_ng_mean(Born_jets);
std::poisson_distribution<int> dist(ng_mean);
const int ng = dist(ran);
assert(ng >= 0);
assert(ng < ng_max);
weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng);
return ng;
}
void PhaseSpacePoint::copy_AWZH_boson_from(Event const & event){
auto const & from = event.outgoing();
const auto AWZH_boson = std::find_if(
begin(from), end(from),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson == end(from)) return;
auto insertion_point = std::lower_bound(
begin(outgoing_), end(outgoing_), *AWZH_boson, rapidity_less{}
);
outgoing_.insert(insertion_point, *AWZH_boson);
// copy decay products
const int idx = std::distance(begin(from), AWZH_boson);
assert(idx >= 0);
const auto decay_it = event.decays().find(idx);
if(decay_it != end(event.decays())){
const int new_idx = std::distance(begin(outgoing_), insertion_point);
assert(new_idx >= 0);
assert(outgoing_[new_idx].type == AWZH_boson->type);
decays_.emplace(new_idx, decay_it->second);
}
assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
}
namespace {
auto get_first_anyquark_emission(Event const & ev) {
// find born quarks (ignore extremal partons)
auto const firstquark = std::find_if(
std::next(ev.begin_partons()), std::prev(ev.end_partons(), 2),
[](Particle const & s){ return (is_anyquark(s)); }
);
// assert that it is a q-q_bar pair.
assert(std::distance(firstquark, ev.end_partons()) != 2);
assert(
( is_quark(*firstquark) && is_antiquark(*std::next(firstquark)) )
|| ( is_antiquark(*firstquark) && is_quark(*std::next(firstquark)) )
);
return firstquark;
}
//! returns index of most backward q-qbar jet
template<class Iterator>
int get_back_quark_jet(Event const & ev, Iterator firstquark){
// find jets at FO corresponding to the quarks
// technically this isn't necessary for LO
std::vector<int> const born_indices{ ev.particle_jet_indices() };
const auto firstquark_idx = std::distance(ev.begin_partons(), firstquark);
int const firstjet_idx = born_indices[firstquark_idx];
assert(firstjet_idx>0);
assert( born_indices[firstquark_idx+1] == firstjet_idx+1 );
return firstjet_idx;
}
//! returns index of most backward q-qbar jet
int getBackQuarkJet(Event const & ev){
const auto firstquark = get_first_anyquark_emission(ev);
return get_back_quark_jet(ev, firstquark);
}
template<class ConstIterator, class Iterator>
void label_extremal_qqx(
ConstIterator born_begin, ConstIterator born_end,
Iterator first_out
){
// find born quarks
const auto firstquark = std::find_if(
born_begin, born_end-1,
[](Particle const & s){ return (is_anyquark(s)); }
);
assert(firstquark != born_end-1);
const auto secondquark = std::find_if(
firstquark+1, born_end,
[](Particle const & s){ return (is_anyquark(s)); }
);
assert(secondquark != born_end);
assert( ( is_quark(*firstquark) && is_antiquark(*secondquark) )
|| ( is_antiquark(*firstquark) && is_quark(*secondquark) ));
assert(first_out->type == ParticleID::gluon);
assert((first_out+1)->type == ParticleID::gluon);
// copy type from born
first_out->type = firstquark->type;
(first_out+1)->type = secondquark->type;
}
}
void PhaseSpacePoint::label_qqx(Event const & event){
assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
assert(filter_partons(outgoing_).size() == outgoing_.size());
if(qqxb_){
label_extremal_qqx( event.outgoing().cbegin(), event.outgoing().cend(),
outgoing_.begin()
);
return;
}
if(qqxf_){ // same as qqxb with reversed order
label_extremal_qqx( event.outgoing().crbegin(), event.outgoing().crend(),
outgoing_.rbegin()
);
return;
}
// central qqx
const auto firstquark = get_first_anyquark_emission(event);
// find jets at FO corresponding to the quarks
// technically this isn't necessary for LO
const auto firstjet_idx = get_back_quark_jet(event, firstquark);
// find corresponding jets after resummation
fastjet::ClusterSequence cs{to_PseudoJet(outgoing_), param_.jet_param.def};
auto const jets = fastjet::sorted_by_rapidity(
cs.inclusive_jets( param_.jet_param.min_pt ));
std::vector<int> const resum_indices{ cs.particle_jet_indices({jets}) };
// assert that jets didn't move
assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx )->rapidity(),
jets[ firstjet_idx ].rapidity(), 1e-2) );
assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx+1 )->rapidity(),
jets[ firstjet_idx+1 ].rapidity(), 1e-2) );
// find last partons in first (central) jet
size_t idx_out = 0;
for(size_t i=resum_indices.size()-2; i>0; --i)
if(resum_indices[i] == firstjet_idx){
idx_out = i;
break;
}
assert(idx_out != 0);
// check that there is sufficient pt in jets from the quarks
const double minpartonjetpt = 1. - param_.max_ext_soft_pt_fraction;
if (outgoing_[idx_out].p.pt()<minpartonjetpt*( event.jets().cbegin()+firstjet_idx )->pt()){
weight_=0.;
status_ = StatusCode::wrong_jets;
return;
}
if (outgoing_[idx_out+1].p.pt()<minpartonjetpt*( event.jets().cbegin()+firstjet_idx+1 )->pt()){
weight_=0.;
status_ = StatusCode::wrong_jets;
return;
}
// check that no additional emission between jets
// such configurations are possible if we have an gluon gets generated
// inside the rapidities of the qqx chain, but clusted to a
// differnet/outside jet. Changing this is non trivial
if(resum_indices[idx_out+1] != resum_indices[idx_out]+1){
weight_=0.;
status_ = StatusCode::gluon_in_qqx;
return;
}
outgoing_[idx_out].type = firstquark->type;
outgoing_[idx_out+1].type = std::next(firstquark)->type;
}
void PhaseSpacePoint::label_quarks(Event const & ev){
const auto WZEmit = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & s){ return (std::abs(s.type) == pid::Wp || s.type == pid::Z_photon_mix); }
);
if (WZEmit != end(ev.outgoing())){
if(!qqxb_) {
const size_t backward_FKL_idx = unob_?1:0;
const auto backward_FKL = std::next(ev.begin_partons(), backward_FKL_idx);
outgoing_[backward_FKL_idx].type = backward_FKL->type;
}
if(!qqxf_) {
const size_t forward_FKL_idx = unof_?1:0;
const auto forward_FKL = std::prev(ev.end_partons(), 1+forward_FKL_idx);
outgoing_.rbegin()[unof_].type = forward_FKL->type;
}
} else {
most_backward_FKL(outgoing_).type = ev.incoming().front().type;
most_forward_FKL(outgoing_).type = ev.incoming().back().type;
}
if(qqxmid_||qqxb_||qqxf_){
label_qqx(ev);
}
}
PhaseSpacePoint::PhaseSpacePoint(
Event const & ev, PhaseSpacePointConfig conf, RNG & ran
):
unob_{ev.type() == event_type::unob},
unof_{ev.type() == event_type::unof},
qqxb_{ev.type() == event_type::qqxexb},
qqxf_{ev.type() == event_type::qqxexf},
qqxmid_{ev.type() == event_type::qqxmid},
param_{std::move(conf)},
status_{unspecified}
{
weight_ = 1;
const auto & Born_jets = ev.jets();
const int ng = sample_ng(Born_jets, ran);
weight_ /= std::tgamma(ng + 1);
const int ng_jets = sample_ng_jets(ng, Born_jets, ran);
std::vector<fastjet::PseudoJet> out_partons = gen_non_jet(
ng - ng_jets, CMINPT, param_.jet_param.min_pt, ran
);
int qqxbackjet(-1);
if(qqxmid_){
qqxbackjet = getBackQuarkJet(ev);
}
const auto qperp = std::accumulate(
begin(out_partons), end(out_partons),
fastjet::PseudoJet{}
);
const auto jets = reshuffle(Born_jets, qperp);
if(weight_ == 0.) {
status_ = failed_reshuffle;
return;
}
if(! pass_resummation_cuts(jets)){
status_ = failed_resummation_cuts;
weight_ = 0.;
return;
}
std::vector<fastjet::PseudoJet> jet_partons = split(
jets, ng_jets, qqxbackjet, ran
);
if(weight_ == 0.) {
status_ = StatusCode::failed_split;
return;
}
if(qqxmid_){
rescale_qqx_rapidities(
out_partons, jets,
most_backward_FKL(jet_partons).rapidity(),
most_forward_FKL(jet_partons).rapidity(),
qqxbackjet
);
}
else{
rescale_rapidities(
out_partons,
most_backward_FKL(jet_partons).rapidity(),
most_forward_FKL(jet_partons).rapidity()
);
}
if(! cluster_jets(out_partons).empty()){
weight_ = 0.;
status_ = StatusCode::empty_jets;
return;
}
std::sort(begin(out_partons), end(out_partons), rapidity_less{});
assert(
std::is_sorted(begin(jet_partons), end(jet_partons), rapidity_less{})
);
const auto first_jet_parton = out_partons.insert(
end(out_partons), begin(jet_partons), end(jet_partons)
);
std::inplace_merge(
begin(out_partons), first_jet_parton, end(out_partons), rapidity_less{}
);
if(! jets_ok(Born_jets, out_partons)){
weight_ = 0.;
status_ = StatusCode::wrong_jets;
return;
}
weight_ *= phase_space_normalisation(Born_jets.size(), out_partons.size());
outgoing_.reserve(out_partons.size() + 1); // one slot for possible A, W, Z, H
for(auto & p: out_partons){
outgoing_.emplace_back(Particle{pid::gluon, std::move(p), {}});
}
assert(!outgoing_.empty());
label_quarks(ev);
if(weight_ == 0.) {
//! @TODO optimise s.t. this is not possible
// status is handled internally
return;
}
copy_AWZH_boson_from(ev);
reconstruct_incoming(ev.incoming());
status_ = StatusCode::good;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_non_jet(
int count, double ptmin, double ptmax, RNG & ran
){
// heuristic parameters for pt sampling
const double ptpar = 1.3 + count/5.;
const double temp1 = std::atan((ptmax - ptmin)/ptpar);
std::vector<fastjet::PseudoJet> partons(count);
for(size_t i = 0; i < static_cast<size_t>(count); ++i){
const double r1 = ran.flat();
const double pt = ptmin + ptpar*std::tan(r1*temp1);
const double temp2 = std::cos(r1*temp1);
const double phi = 2*M_PI*ran.flat();
weight_ *= 2.0*M_PI*pt*ptpar*temp1/(temp2*temp2);
// we don't know the allowed rapidity span yet,
// set a random value to be rescaled later on
const double y = ran.flat();
partons[i].reset_PtYPhiM(pt, y, phi);
// Set user index higher than any jet-parton index
// in order to assert that these are not inside jets
partons[i].set_user_index(i + 1 + ng_max);
assert(ptmin-1e-5 <= partons[i].pt() && partons[i].pt() <= ptmax+1e-5);
}
assert(std::all_of(partons.cbegin(), partons.cend(), is_nonjet_parton));
return sorted_by_rapidity(partons);
}
void PhaseSpacePoint::rescale_qqx_rapidities(
std::vector<fastjet::PseudoJet> & out_partons,
std::vector<fastjet::PseudoJet> const & jets,
const double ymin1, const double ymax2,
const int qqxbackjet
){
const double ymax1 = jets[qqxbackjet].rapidity();
const double ymin2 = jets[qqxbackjet+1].rapidity();
constexpr double ep = 1e-7;
const double tot_y = ymax1 - ymin1 + ymax2 - ymin2;
std::vector<std::reference_wrapper<fastjet::PseudoJet>> refpart(
out_partons.begin(), out_partons.end());
double ratio = (ymax1 - ymin1)/tot_y;
const auto gap{ std::find_if(refpart.begin(), refpart.end(),
[ratio](fastjet::PseudoJet p){
return (p.rapidity()>=ratio);} ) };
double ymin = ymin1;
double ymax = ymax1;
double dy = ymax - ymin - 2*ep;
double offset = 0.;
for(auto it_part=refpart.begin(); it_part<refpart.end(); ++it_part){
if(it_part == gap){
ymin = ymin2;
ymax = ymax2;
dy = ymax - ymin - 2*ep;
offset = ratio;
ratio = 1-ratio;
}
fastjet::PseudoJet & part = *it_part;
assert(offset <= part.rapidity() && part.rapidity() < ratio+offset);
const double y = ymin + ep + dy*((part.rapidity()-offset)/ratio);
part.reset_momentum_PtYPhiM(part.pt(), y, part.phi());
weight_ *= tot_y-4.*ep;
assert(ymin <= part.rapidity() && part.rapidity() <= ymax);
}
assert(is_sorted(begin(out_partons), end(out_partons), rapidity_less{}));
}
void PhaseSpacePoint::rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
){
constexpr double ep = 1e-7;
for(auto & parton: partons){
assert(0 <= parton.rapidity() && parton.rapidity() <= 1);
const double dy = ymax - ymin - 2*ep;
const double y = ymin + ep + dy*parton.rapidity();
parton.reset_momentum_PtYPhiM(parton.pt(), y, parton.phi());
weight_ *= dy;
assert(ymin <= parton.rapidity() && parton.rapidity() <= ymax);
}
}
namespace {
template<typename T, typename... Rest>
auto min(T const & a, T const & b, Rest&&... r) {
using std::min;
return min(a, min(b, std::forward<Rest>(r)...));
}
}
double PhaseSpacePoint::probability_in_jet(
std::vector<fastjet::PseudoJet> const & Born_jets
) const{
assert(std::is_sorted(begin(Born_jets), end(Born_jets), rapidity_less{}));
assert(Born_jets.size() >= 2);
const double dy =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
const double R = param_.jet_param.def.R();
const int njets = Born_jets.size();
const double p_J_y_large = (njets-1)*R*R/(2.*dy);
const double p_J_y0 = njets*R/M_PI;
return min(p_J_y_large, p_J_y0, 1.);
}
int PhaseSpacePoint::sample_ng_jets(
int ng, std::vector<fastjet::PseudoJet> const & Born_jets, RNG & ran
){
const double p_J = probability_in_jet(Born_jets);
std::binomial_distribution<> bin_dist(ng, p_J);
const int ng_J = bin_dist(ran);
weight_ *= std::pow(p_J, -ng_J)*std::pow(1 - p_J, ng_J - ng);
return ng_J;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::reshuffle(
std::vector<fastjet::PseudoJet> const & Born_jets,
fastjet::PseudoJet const & q
){
if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets;
const auto jets = resummation_jet_momenta(Born_jets, q);
if(jets.empty()){
weight_ = 0;
return {};
}
// additional Jacobian to ensure Born integration over delta gives 1
weight_ *= resummation_jet_weight(Born_jets, q);
return jets;
}
std::vector<int> PhaseSpacePoint::distribute_jet_partons(
int ng_jets, std::vector<fastjet::PseudoJet> const & jets, RNG & ran
){
size_t first_valid_jet = 0;
size_t num_valid_jets = jets.size();
const double R_eff = 5./3.*param_.jet_param.def.R();
// if there is an unordered jet too far away from the FKL jets
// then extra gluon constituents of the unordered jet would
// violate the FKL rapidity ordering
if((unob_||qqxb_) && jets[0].delta_R(jets[1]) > R_eff){
++first_valid_jet;
--num_valid_jets;
}
else if((unof_||qqxf_) && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){
--num_valid_jets;
}
std::vector<int> np(jets.size(), 1);
for(int i = 0; i < ng_jets; ++i){
++np[first_valid_jet + ran.flat() * num_valid_jets];
}
weight_ *= std::pow(num_valid_jets, ng_jets);
return np;
}
#ifndef NDEBUG
namespace{
bool tagged_FKL_backward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return std::find_if(
begin(jet_partons), end(jet_partons),
[](fastjet::PseudoJet const & p){
return p.user_index() == backward_FKL_uid;
}
) != end(jet_partons);
}
bool tagged_FKL_forward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
// the most forward FKL parton is most likely near the end of jet_partons;
// start search from there
return std::find_if(
jet_partons.rbegin(), jet_partons.rend(),
[](fastjet::PseudoJet const & p){
return p.user_index() == forward_FKL_uid;
}
) != jet_partons.rend();
}
bool tagged_FKL_extremal(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return tagged_FKL_backward(jet_partons) && tagged_FKL_forward(jet_partons);
}
} // namespace anonymous
#endif
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
int ng_jets, size_t qqxbackjet, RNG & ran
){
return split(
jets, distribute_jet_partons(ng_jets, jets, ran), qqxbackjet, ran);
}
bool PhaseSpacePoint::pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const{
if(ext_parton.pt() < param_.min_extparton_pt) return false;
return (ext_parton - jet).pt()/jet.pt() < param_.max_ext_soft_pt_fraction;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np,
size_t qqxbackjet,
RNG & ran
){
assert(! jets.empty());
assert(jets.size() == np.size());
assert(pass_resummation_cuts(jets));
const size_t most_backward_FKL_idx = 0 + unob_ + qqxb_;
const size_t most_forward_FKL_idx = jets.size() - 1 - unof_ - qqxf_;
const auto & jet = param_.jet_param;
const JetSplitter jet_splitter{jet.def, jet.min_pt};
std::vector<fastjet::PseudoJet> jet_partons;
// randomly distribute jet gluons among jets
for(size_t i = 0; i < jets.size(); ++i){
auto split_res = jet_splitter.split(jets[i], np[i], ran);
weight_ *= split_res.weight;
if(weight_ == 0) return {};
assert(
std::all_of(
begin(split_res.constituents), end(split_res.constituents),
is_jet_parton
)
);
const auto first_new_parton = jet_partons.insert(
end(jet_partons),
begin(split_res.constituents), end(split_res.constituents)
);
// mark uno and extremal FKL emissions here so we can check
// their position once all emissions are generated
// also mark qqxmid partons, and apply appropriate pt cut.
auto extremal = end(jet_partons);
if (i == most_backward_FKL_idx){ //FKL backward emission
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(backward_FKL_uid);
}
else if(((unob_ || qqxb_) && i == 0)){
// unordered/qqxb
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index((unob_)?unob_uid:qqxb_uid);
}
else if (i == most_forward_FKL_idx){
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(forward_FKL_uid);
}
else if(((unof_ || qqxf_) && i == jets.size() - 1)){
// unordered/qqxf
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index((unof_)?unof_uid:qqxf_uid);
}
else if((qqxmid_ && i == qqxbackjet)){
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(qqxmid1_uid);
}
else if((qqxmid_ && i == qqxbackjet+1)){
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(qqxmid2_uid);
}
if(
extremal != end(jet_partons)
&& !pass_extremal_cuts(*extremal, jets[i])
){
weight_ = 0;
return {};
}
}
assert(tagged_FKL_extremal(jet_partons));
std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
if(
!extremal_ok(jet_partons)
|| !split_preserved_jets(jets, jet_partons)
){
weight_ = 0.;
return {};
}
return jet_partons;
}
bool PhaseSpacePoint::extremal_ok(
std::vector<fastjet::PseudoJet> const & partons
) const{
assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
if(unob_ && partons.front().user_index() != unob_uid) return false;
if(unof_ && partons.back().user_index() != unof_uid) return false;
if(qqxb_ && partons.front().user_index() != qqxb_uid) return false;
if(qqxf_ && partons.back().user_index() != qqxf_uid) return false;
return
most_backward_FKL(partons).user_index() == backward_FKL_uid
&& most_forward_FKL(partons).user_index() == forward_FKL_uid;
}
bool PhaseSpacePoint::split_preserved_jets(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<fastjet::PseudoJet> const & jet_partons
) const{
assert(std::is_sorted(begin(jets), end(jets), rapidity_less{}));
const auto split_jets = cluster_jets(jet_partons);
// this can happen if two overlapping jets
// are both split into more than one parton
if(split_jets.size() != jets.size()) return false;
for(size_t i = 0; i < split_jets.size(); ++i){
// this can happen if there are two overlapping jets
// and a parton is assigned to the "wrong" jet
if(!nearby_ep(jets[i].rapidity(), split_jets[i].rapidity(), 1e-2)){
return false;
}
}
return true;
}
template<class Particle>
Particle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> const & partons
) const{
return partons[0 + unob_ + qqxb_];
}
template<class Particle>
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
const size_t idx = partons.size() - 1 - unof_ - qqxf_;
assert(idx < partons.size());
return partons[idx];
}
template<class Particle>
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
return partons[0 + unob_ + qqxb_];
}
template<class Particle>
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
const size_t idx = partons.size() - 1 - unof_ - qqxf_;
assert(idx < partons.size());
return partons[idx];
}
bool PhaseSpacePoint::contains_idx(
fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
) const {
auto const & constituents = jet.constituents();
const int idx = parton.user_index();
const bool injet = std::find_if(
begin(constituents), end(constituents),
[idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
) != end(constituents);
const double minpartonjetpt = 1. - param_.max_ext_soft_pt_fraction;
return ((parton.pt()>minpartonjetpt*jet.pt())&&injet);
}
bool PhaseSpacePoint::jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
if(jets.size() != Born_jets.size()) return false;
int in_jet = 0;
for(size_t i = 0; i < jets.size(); ++i){
assert(jets[i].has_constituents());
for(auto && parton: jets[i].constituents()){
if(is_nonjet_parton(parton)) return false;
}
in_jet += jets[i].constituents().size();
}
const int expect_in_jet = std::count_if(
partons.cbegin(), partons.cend(), is_jet_parton
);
if(in_jet != expect_in_jet) return false;
// note that PseudoJet::contains does not work here
if(! (
contains_idx(most_backward_FKL(jets), most_backward_FKL(partons))
&& contains_idx(most_forward_FKL(jets), most_forward_FKL(partons))
)) return false;
if(unob_ && !contains_idx(jets.front(), partons.front())) return false;
if(qqxb_ && !contains_idx(jets.front(), partons.front())) return false;
if(unof_ && !contains_idx(jets.back(), partons.back())) return false;
if(qqxf_ && !contains_idx(jets.back(), partons.back())) return false;
#ifndef NDEBUG
for(size_t i = 0; i < jets.size(); ++i){
assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2));
}
#endif
return true;
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<Particle, 2> const & Born_incoming
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
for(size_t i = 0; i < incoming_.size(); ++i){
incoming_[i].type = Born_incoming[i].type;
}
assert(momentum_conserved());
}
double PhaseSpacePoint::phase_space_normalisation(
int num_Born_jets, int num_out_partons
) const{
return std::pow(16*std::pow(M_PI,3), num_Born_jets - num_out_partons);
}
bool PhaseSpacePoint::momentum_conserved() const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
const double norm = diff.E();
for(auto const & out: outgoing()) diff -= out.p;
return nearby(diff, fastjet::PseudoJet{}, norm);
}
} //namespace HEJ
diff --git a/src/RivetAnalysis.cc b/src/RivetAnalysis.cc
index a9e9d5e..c183288 100644
--- a/src/RivetAnalysis.cc
+++ b/src/RivetAnalysis.cc
@@ -1,186 +1,186 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/RivetAnalysis.hh"
#include "HEJ/ConfigFlags.hh"
#ifdef HEJ_BUILD_WITH_RIVET
#include <ostream>
#include <cstddef>
#include "yaml-cpp/yaml.h"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/Config/RivetConfig.hh"
#ifdef RIVET_ENABLE_HEPMC_3
#include "HepMC3/GenEvent.h"
#include "HEJ/HepMC3Interface.hh"
#else // rivet with HepMC 2
#include "HepMC/GenEvent.h"
#include "HEJ/HepMC2Interface.hh"
#endif
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#endif
namespace HEJ{
std::unique_ptr<Analysis> RivetAnalysis::create(
YAML::Node const & config, LHEF::HEPRUP const & heprup
){
return std::make_unique<RivetAnalysis>(config, heprup);
}
}
#ifdef HEJ_BUILD_WITH_RIVET
namespace HEJ {
#ifdef RIVET_ENABLE_HEPMC_3
using HepMCInterface=HepMC3Interface;
#else
using HepMCInterface=HepMC2Interface;
#endif
struct RivetAnalysis::rivet_info {
rivet_info(std::unique_ptr<Rivet::AnalysisHandler> && h,
- std::string && s, HEJ::HepMCInterface && m):
+ std::string && s, HepMCInterface && m):
handler{std::move(h)}, name{std::move(s)}, hepmc{std::move(m)}
{}
// AnalysisHandler doesn't allow a copy constructor -> use ptr
std::unique_ptr<Rivet::AnalysisHandler> handler;
std::string name;
- HEJ::HepMCInterface hepmc;
+ HepMCInterface hepmc;
};
RivetAnalysis::RivetAnalysis(YAML::Node const & config,
LHEF::HEPRUP const & heprup
):
heprup_{heprup},
first_event_(true)
{
// assert that only supported options are provided
if(!config.IsMap()) throw invalid_type{"Rivet config has to be a map!"};
for(auto const & entry: config){
const auto name = entry.first.as<std::string>();
if(name != "output" && name != "rivet")
throw unknown_option{"Unknown option \'"+name+"\' provided to rivet."};
}
// get output name
if(!config["output"])
throw std::invalid_argument{
"No output was provided to rivet. "
"Either give an output or deactivate rivet."
};
if(!config["output"].IsScalar())
throw invalid_type{
"Output for rivet must be a scalar."
};
output_name_ = config["output"].as<std::string>();
// read in analyses
const auto & name_node = config["rivet"];
switch(name_node.Type()){
case YAML::NodeType::Scalar:
analyses_names_.push_back(name_node.as<std::string>());
break;
case YAML::NodeType::Sequence:
for(YAML::const_iterator it = name_node.begin(); it != name_node.end(); ++it){
analyses_names_.push_back(it->as<std::string>());
}
break;
default:
throw std::invalid_argument{
"No analysis was provided to rivet. "
"Either give an analysis or deactivate rivet."
};
}
}
// it is a bit ugly that we can't do this directly in `initialise`
void RivetAnalysis::init(Event const & event){
#ifdef HEJ_USE_RIVET2
rivet_runs_.reserve(event.variations().size()+1);
#else
(void) event; // shut up compiler
#endif
rivet_runs_.emplace_back( std::make_unique<Rivet::AnalysisHandler>(),
std::string(""), HepMCInterface(heprup_)
);
rivet_runs_.back().handler->addAnalyses(analyses_names_);
#ifdef HEJ_USE_RIVET2
//! scale variation in rivet 2 through multiple analyses
if( !event.variations().empty() ){
for(auto const & vari : event.variations()){
rivet_runs_.emplace_back( std::make_unique<Rivet::AnalysisHandler>(),
"."+to_simple_string(*vari.description), HepMCInterface(heprup_)
);
rivet_runs_.back().handler->addAnalyses(analyses_names_);
}
}
#endif
}
void RivetAnalysis::fill(Event const & event, Event const &){
if(first_event_){
first_event_=false;
init(event);
}
auto hepmc_event(rivet_runs_.front().hepmc(event));
rivet_runs_.front().handler->analyze(hepmc_event);
#ifdef HEJ_USE_RIVET2
for(std::size_t i = 1; i < rivet_runs_.size(); ++i){
auto & run = rivet_runs_[i];
run.hepmc.set_central(hepmc_event, event, i-1);
run.handler->analyze(hepmc_event);
}
#endif
}
void RivetAnalysis::finalise(){
for(auto & run: rivet_runs_){
run.handler->finalize();
run.handler->writeData(output_name_+run.name+std::string(".yoda"));
}
}
} // namespace HEJ
#else // no rivet => create empty analysis
namespace Rivet {
class AnalysisHandler {};
}
namespace HEJ {
struct RivetAnalysis::rivet_info{};
RivetAnalysis::RivetAnalysis(YAML::Node const &, LHEF::HEPRUP const &){
throw std::invalid_argument(
"Failed to create RivetAnalysis: "
"HEJ 2 was built without rivet support"
);
}
void RivetAnalysis::init(Event const &){}
void RivetAnalysis::fill(Event const &, Event const &){}
void RivetAnalysis::finalise(){}
} // namespace HEJ
#endif
namespace HEJ {
RivetAnalysis::~RivetAnalysis() = default;
}
diff --git a/src/ScaleFunction.cc b/src/ScaleFunction.cc
index 04eb410..32d5047 100644
--- a/src/ScaleFunction.cc
+++ b/src/ScaleFunction.cc
@@ -1,171 +1,171 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/ScaleFunction.hh"
#include <cassert>
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
-namespace HEJ{
+namespace HEJ {
double H_T(Event const & ev){
double result = 0.;
for(size_t i = 0; i < ev.outgoing().size(); ++i){
auto const decay_products = ev.decays().find(i);
if(decay_products == end(ev.decays())){
result += ev.outgoing()[i].perp();
}
else{
for(auto const & particle: decay_products->second){
result += particle.perp();
}
}
}
return result;
}
double max_jet_pt(Event const & ev) {
return sorted_by_pt(ev.jets()).front().pt();
}
double jet_invariant_mass(Event const & ev) {
fastjet::PseudoJet sum;
for(const auto & jet: ev.jets()) sum+=jet;
return sum.m();
}
double m_j1j2(Event const & ev) {
const auto jets = sorted_by_pt(ev.jets());
assert(jets.size() >= 2);
return (jets[0] + jets[1]).m();
}
ScaleFunction operator*(double factor, ScaleFunction base_scale) {
base_scale.name_.insert(0, std::to_string(factor) + '*');
auto new_fun =
- [factor,fun{std::move(base_scale.fun_)}](HEJ::Event const & ev) {
+ [factor,fun{std::move(base_scale.fun_)}](Event const & ev) {
return factor*fun(ev);
};
base_scale.fun_ = std::move(new_fun);
return base_scale;
}
ScaleFunction operator*(ScaleFunction factor, ScaleFunction base_scale) {
base_scale.name_.insert(0, factor.name_ + '*');
auto new_fun =
[fun1{std::move(factor.fun_)}, fun2{std::move(base_scale.fun_)}, name{base_scale.name_}]
- (HEJ::Event const & ev) {
+ (Event const & ev) {
return fun1(ev)*fun2(ev);
};
base_scale.fun_ = std::move(new_fun);
return base_scale;
}
ScaleFunction operator/(ScaleFunction base_scale, double denom) {
base_scale.name_.append('/' + std::to_string(denom));
auto new_fun =
- [denom,fun{std::move(base_scale.fun_)}](HEJ::Event const & ev) {
+ [denom,fun{std::move(base_scale.fun_)}](Event const & ev) {
return fun(ev)/denom;
};
base_scale.fun_ = std::move(new_fun);
return base_scale;
}
ScaleFunction operator/(ScaleFunction base_scale, ScaleFunction denom) {
base_scale.name_.append('/' + denom.name_);
auto new_fun =
[fun2{std::move(denom.fun_)}, fun1{std::move(base_scale.fun_)}]
- (HEJ::Event const & ev) {
+ (Event const & ev) {
return fun1(ev)/fun2(ev);
};
base_scale.fun_ = std::move(new_fun);
return base_scale;
}
// TODO: significant logic duplication with operator()
void ScaleGenerator::gen_descriptions() {
if(scales_.empty()) {
throw std::logic_error{"Need at least one scale"};
}
descriptions_.emplace_back(
std::make_shared<ParameterDescription>(scales_.front().name(), 1., 1.)
);
for(auto & base_scale: scales_){
const auto base_name = base_scale.name();
descriptions_.emplace_back(
std::make_shared<ParameterDescription>(base_name, 1., 1.)
);
//multiplicative scale variation
for(double mur_factor: scale_factors_){
for(double muf_factor: scale_factors_){
if(muf_factor == 1. && mur_factor == 1.) continue;
if(
mur_factor/muf_factor < 1/max_scale_ratio_
|| mur_factor/muf_factor > max_scale_ratio_
) continue;
descriptions_.emplace_back(
std::make_shared<ParameterDescription>(
base_name, mur_factor, muf_factor
)
);
}
}
}
}
Event ScaleGenerator::operator()(Event ev) const {
if(! ev.variations().empty()) {
throw std::invalid_argument{"Input event already has scale variation"};
}
assert(!scales_.empty());
assert(!descriptions_.empty());
size_t descr_idx = 0;
const double mu_central = (scales_.front())(ev);
ev.central().mur = mu_central;
ev.central().muf = mu_central;
assert(descr_idx < descriptions_.size());
assert(descriptions_[descr_idx]);
ev.central().description = descriptions_[descr_idx++];
// check if we are doing scale variation at all
if(scales_.size() == 1 && scale_factors_.empty()) return ev;
for(auto & base_scale: scales_){
const double mu_base = base_scale(ev);
assert(descr_idx < descriptions_.size());
assert(descriptions_[descr_idx]);
ev.variations().emplace_back(
EventParameters{
mu_base, mu_base, ev.central().weight, descriptions_[descr_idx++]
}
);
//multiplicative scale variation
for(double mur_factor: scale_factors_){
const double mur = mu_base*mur_factor;
for(double muf_factor: scale_factors_){
if(muf_factor == 1. && mur_factor == 1.) continue;
const double muf = mu_base*muf_factor;
if(
mur/muf < 1/max_scale_ratio_
|| mur/muf > max_scale_ratio_
) continue;
assert(descr_idx < descriptions_.size());
assert(descriptions_[descr_idx]);
ev.variations().emplace_back(
EventParameters{
mur, muf, ev.central().weight, descriptions_[descr_idx++]
}
);
}
}
};
return ev;
}
-}
+} // namespace HEJ
diff --git a/src/Wjets.cc b/src/Wjets.cc
index d941efb..6abe7dd 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,796 +1,796 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
#include "HEJ/Wjets.hh"
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <iostream>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/jets.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/exceptions.hh"
// generated headers
#include "HEJ/currents/jV_j.hh"
#include "HEJ/currents/jVuno_j.hh"
#include "HEJ/currents/jWqqbar_j.hh"
#include "HEJ/currents/jW_qqbar_j.hh"
#include "HEJ/currents/j_Wqqbar_j.hh"
#include "HEJ/currents/jW_jqqbar.hh"
#include "HEJ/currents/jV_juno.hh"
namespace HEJ {
namespace currents {
namespace {
// short hands
using std::abs;
using std::conj;
using std::pow;
// --- Helper Functions ---
// FKL W Helper Functions
double WProp(const HLV & plbar, const HLV & pl,
ParticleProperties const & wprop
){
COM propW = COM(0.,-1.)/( (pl+plbar).m2() - wprop.mass*wprop.mass
+ COM(0.,1.)*wprop.mass*wprop.width);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
/**
* @brief Contraction of W + unordered jet current with FKL current
*
* See eq:wunocurrent in the developer manual for the definition
* of the W + unordered jet current
*/
template<Helicity h1, Helicity h2, Helicity pol>
double jM2Wuno(
HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
HLV const & pa, HLV const & p2, HLV const & pb,
ParticleProperties const & wprop
){
using helicity::minus;
const COM u1 = U1<h1, minus, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
const COM u2 = U2<h1, minus, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
const COM l = L<h1, minus, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
const COM x = u1 - l;
const COM y = u2 + l;
double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real();
amp *= WProp(plbar, pl, wprop);
const HLV q1g = pa-pl-plbar-p1-pg;
const HLV q2 = p2-pb;
const double t1 = q1g.m2();
const double t2 = q2.m2();
amp /= (t1*t2);
return amp;
}
/**
* @brief Contraction of W + extremal qqbar current with FKL current
*
* See eq:crossed in the developer manual for the definition
* of the W + extremal qqbar current.
*
*/
template<Helicity h1, Helicity h2, Helicity hg>
double jM2_W_extremal_qqbar(
HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl,
HLV const & pqbar, HLV const & p3, HLV const & pb,
ParticleProperties const & wprop
){
const COM u1Xcontr = U1X<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM u2Xcontr = U2X<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM lXcontr = LX<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM x = u1Xcontr - lXcontr;
const COM y = u2Xcontr + lXcontr;
double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real();
amp *= WProp(plbar, pl, wprop);
const HLV q1 = pg-pl-plbar-pq-pqbar;
const HLV q2 = p3-pb;
const double t1 = q1.m2();
const double t2 = q2.m2();
amp /= (t1*t2);
return amp;
}
//! 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 const & p1out, HLV plbar, HLV pl,
HLV const & p1in, HLV const & p2out, HLV const & p2in,
bool aqlineb, bool /* aqlinef */,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out);
const double WPropfact = WProp(plbar, pl, wprop);
// we assume that the W is emitted off a quark line
// if this is not the case, we have to apply CP conjugation,
// which is equivalent to swapping lepton and antilepton momenta
if(aqlineb) std::swap(pl, plbar);
- const COM ampm = HEJ::jV_j<minus, minus, minus>(p1in,p1out,p2in,p2out,pl,plbar);
- const COM ampp = HEJ::jV_j<minus, minus, plus>(p1in,p1out,p2in,p2out,pl,plbar);
+ const COM ampm = jV_j<minus, minus, minus>(p1in,p1out,p2in,p2out,pl,plbar);
+ const COM ampp = jV_j<minus, minus, plus>(p1in,p1out,p2in,p2out,pl,plbar);
const double Msqr = std::norm(ampm) + std::norm(ampp);
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return C_F*C_F*WPropfact*Msqr/(q1.m2()*q2.m2()*(N_C*N_C - 1)*4);
}
} // Anonymous Namespace
double ME_W_qQ(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop);
}
double ME_W_qQbar(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true, wprop);
}
double ME_W_qbarQ(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop);
}
double ME_W_qbarQbar(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true, wprop);
}
double ME_W_qg(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop)
*K_g(p2out, p2in)/C_F;
}
double ME_W_qbarg(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop)
*K_g(p2out, p2in)/C_F;
}
namespace {
// helicity amplitude squares for contraction of W current with unordered
// current
template<Helicity h2, Helicity hg>
double ampsq_jW_juno(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & p2,
HLV const & pg,
HLV const & pl, HLV const & plbar
){
using helicity::minus;
const COM u1 = U1_jV<minus, minus, h2, hg>(pa,p1,pb,p2,pg,pl,plbar);
const COM u2 = U2_jV<minus, minus, h2, hg>(pa,p1,pb,p2,pg,pl,plbar);
const COM l = L_jV<minus, minus, h2, hg>(pa,p1,pb,p2,pg,pl,plbar);
return 2.*C_F*std::real((l-u1)*std::conj(l+u2))
+ 2.*C_F*C_F/3.*std::norm(u1+u2)
;
}
/**
* @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.
*
* @TODO use const & for pl plbar
*/
double jW_juno(
HLV const & p1out, HLV plbar, HLV pl, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pg, bool aqlineb,
bool /* aqlinef */,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
// we assume that the W is emitted off a quark line
// if this is not the case, we have to apply CP conjugation,
// which is equivalent to swapping lepton and antilepton momenta
if(aqlineb) std::swap(pl, plbar);
// helicity assignments assume aqlinef == false
// in the end, this is irrelevant since we sum over all helicities
const double ampsq =
+ ampsq_jW_juno<minus, minus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
+ ampsq_jW_juno<minus, plus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
+ ampsq_jW_juno<plus, minus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
+ ampsq_jW_juno<plus, plus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out-pg);
return WProp(plbar, pl, wprop)*ampsq/((16)*(q2.m2()*q1.m2()));
}
} // Anonymous Namespace
double ME_W_unob_qQ(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl, ParticleProperties
const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false, wprop);
}
double ME_W_unob_qQbar(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true, wprop);
}
double ME_W_unob_qbarQ(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false, wprop);
}
double ME_W_unob_qbarQbar(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true, wprop);
}
namespace {
// helicity sum helper for jWuno_j(...)
template<Helicity h1>
double jWuno_j_helsum(
HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
HLV const & pa, HLV const & p2, HLV const & pb,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const double ME2h1pp = jM2Wuno<h1, plus, plus>(
pg, p1, plbar, pl, pa, p2, pb, wprop
);
const double ME2h1pm = jM2Wuno<h1, plus, minus>(
pg, p1, plbar, pl, pa, p2, pb, wprop
);
const double ME2h1mp = jM2Wuno<h1, minus, plus>(
pg, p1, plbar, pl, pa, p2, pb, wprop
);
const double ME2h1mm = jM2Wuno<h1, minus, minus>(
pg, p1, plbar, pl, pa, p2, pb, wprop
);
return ME2h1pp + ME2h1pm + ME2h1mp + ME2h1mm;
}
/**
* @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 const & pg, HLV const & p1out, HLV const & plbar, HLV const & pl,
HLV const & p1in, HLV const & p2out, HLV const & p2in, bool aqlineb,
ParticleProperties const & wprop
){
const auto helsum = aqlineb?
jWuno_j_helsum<helicity::plus>:
jWuno_j_helsum<helicity::minus>;
//Helicity sum and average over initial states
return helsum(pg, p1out,plbar,pl,p1in,p2out,p2in, wprop)/(4.*C_A*C_A);
}
} // Anonymous Namespace
double ME_Wuno_qQ(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qQbar(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qbarQ(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qbarQbar(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qg(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop)
*K_g(p2out, p2in)/C_F;
}
double ME_Wuno_qbarg(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop)
*K_g(p2out, p2in)/C_F;
}
// helicity sum helper for jWqqx_j(...)
template<Helicity h1>
double jWqqx_j_helsum(
HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl,
HLV const & pqbar, HLV const & p3, HLV const & pb,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const double ME2h1pp = jM2_W_extremal_qqbar<h1, plus, plus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double ME2h1pm = jM2_W_extremal_qqbar<h1, plus, minus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double ME2h1mp = jM2_W_extremal_qqbar<h1, minus, plus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double ME2h1mm = jM2_W_extremal_qqbar<h1, minus, minus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
return ME2h1pp + ME2h1pm + ME2h1mp + ME2h1mm;
}
/**
* @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 const & pgin, HLV const & pqout, HLV const & plbar, HLV const & pl,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in, bool aqlinef,
ParticleProperties const & wprop
){
const auto helsum = aqlinef?
jWqqx_j_helsum<helicity::plus>:
jWqqx_j_helsum<helicity::minus>;
//Helicity sum and average over initial states.
double ME2 = helsum(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, wprop)/
(4.*C_A*C_A);
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
double ME_WExqqx_qbarqQ(
HLV const & pgin, HLV const & pqout, HLV const & plbar, HLV const & pl,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop);
}
double ME_WExqqx_qqbarQ(
HLV const & pgin, HLV const & pqbarout, HLV const & plbar, HLV const & pl,
HLV const & pqout, HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop);
}
double ME_WExqqx_qbarqg(
HLV const & pgin, HLV const & pqout, HLV const & plbar, HLV const & pl,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop)
*K_g(p2out,p2in)/C_F;
}
double ME_WExqqx_qqbarg(
HLV const & pgin, HLV const & pqbarout, HLV const & plbar, HLV const & pl,
HLV const & pqout, HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop)
*K_g(p2out,p2in)/C_F;
}
namespace {
template<Helicity h1, Helicity hg>
double jW_jqqbar(
HLV const & pb,
HLV const & p2,
HLV const & p3,
HLV const & pa,
HLV const & p1,
HLV const & pl,
HLV const & plbar
){
static constexpr COM cm1m1 = 8./3.;
static constexpr COM cm2m2 = 8./3.;
static constexpr COM cm3m3 = 6.;
static constexpr COM cm1m2 =-1./3.;
static constexpr COM cm1m3 = COM{0.,-3.};
static constexpr COM cm2m3 = COM{0.,3.};
const COM m1 = jW_qggm1<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
const COM m2 = jW_qggm2<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
const COM m3 = jW_qggm3<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
return real(
cm1m1*pow(abs(m1),2) + cm2m2*pow(abs(m2),2)
+ cm3m3*pow(abs(m3),2) + 2.*real(cm1m2*m1*conj(m2))
)
+ 2.*real(cm1m3*m1*conj(m3))
+ 2.*real(cm2m3*m2*conj(m3)) ;
}
} // Anonymous Namespace
// contraction of W current with extremal qqbar on the other side
double ME_W_Exqqx_QQq(
HLV const & pa, HLV const & pb,
HLV const & p1, HLV const & p2,
HLV const & p3, HLV const & plbar, HLV const & pl, bool aqlinepa,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const double WPropfact = WProp(plbar, pl, wprop);
const double prefactor = 2.*WPropfact/24./4./(
(pa-p1-pl-plbar).m2()*(p2+p3-pb).m2()
);
if(aqlinepa) {
return prefactor*(
jW_jqqbar<plus, minus>(pb,p2,p3,pa,p1,pl,plbar)
+ jW_jqqbar<plus, plus>(pb,p2,p3,pa,p1,pl,plbar)
);
}
return prefactor*(
jW_jqqbar<minus, minus>(pb,p2,p3,pa,p1,pl,plbar)
+ jW_jqqbar<minus, plus>(pb,p2,p3,pa,p1,pl,plbar)
);
}
namespace {
// helper function for matrix element for W + Jets with central qqbar
template<Helicity h1a, Helicity h4b>
double amp_WCenqqx_qq(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & p4,
HLV const & pq, HLV const & pqbar,
HLV const & pl, HLV const & plbar,
HLV const & q11, HLV const & q12
){
const COM sym = M_sym_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
const COM cross = M_cross_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
const COM uncross = M_uncross_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
// Colour factors
static constexpr COM cmsms = 3.;
static constexpr COM cmumu = 4./3.;
static constexpr COM cmcmc = 4./3.;
static constexpr COM cmsmu = COM{0., 3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
static constexpr COM cmumc = -1./6.;
return real(
cmsms*pow(abs(sym),2)
+cmumu*pow(abs(uncross),2)
+cmcmc*pow(abs(cross),2)
)
+2.*real(cmsmu*sym*conj(uncross))
+2.*real(cmsmc*sym*conj(cross))
+2.*real(cmumc*uncross*conj(cross))
;
}
} // Anonymous Namespace
// matrix element for W + Jets with W emission off central qqbar
double ME_WCenqqx_qq(
HLV const & pa, HLV const & pb, HLV const & pl, HLV const & plbar,
std::vector<HLV> const & partons, bool /* aqlinepa */, bool /* aqlinepb */, bool
qqxmarker, int nabove, ParticleProperties const & wprop
){
using helicity::plus;
using helicity::minus;
CLHEP::HepLorentzVector q1 = pa;
for(int i = 0; i <= nabove; ++i){
q1 -= partons[i];
}
const auto qq = split_into_lightlike(q1);
// since q1.m2() < 0 the following assertion is always true
// see documentation for split_into_lightlike
assert(qq.second.e() < 0);
HLV pq, pqbar;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];
}
const HLV p1 = partons.front();
const HLV p4 = partons.back();
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is
// also the choice in qqbar helicity.
// the first helicity label is for aqlinepa == true,
// the second one for aqlinepb == true
// In principle the corresponding helicity should be flipped
// if either of them is false, but we sum up all helicities,
// so this has no effect in the end.
const double amp_mm = amp_WCenqqx_qq<minus, minus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_mp = amp_WCenqqx_qq<minus, plus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pm = amp_WCenqqx_qq<plus, minus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pp = amp_WCenqqx_qq<plus, plus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double t1 = q1.m2();
const double t3 = (q1-pl-plbar-pq-pqbar).m2();
const double amp = WProp(plbar, pl, wprop)*(
amp_mm+amp_mp+amp_pm+amp_pp
)/(9.*4.*t1*t1*t3*t3);
return amp;
}
// helper function for matrix element for W + Jets with central qqbar
// W emitted off extremal parton
// @TODO: code duplication with amp_WCenqqx_qq
template<Helicity h1a, Helicity hqqbar>
double amp_W_Cenqqx_qq(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & pn,
HLV const & pq, HLV const & pqbar,
HLV const & pl, HLV const & plbar,
HLV const & q11, HLV const & q12
){
const COM crossed = M_cross<h1a, hqqbar>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
);
const COM uncrossed = M_qbar<h1a, hqqbar>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
);
const COM sym = M_sym<h1a, hqqbar>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
);
//Colour factors:
static constexpr COM cmsms = 3.;
static constexpr COM cmumu = 4./3.;
static constexpr COM cmcmc = 4./3.;
static constexpr COM cmsmu = COM{0.,3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
static constexpr COM cmumc = -1./6.;
return real(
cmsms*pow(abs(sym),2)
+cmumu*pow(abs(uncrossed),2)
+cmcmc*pow(abs(crossed),2)
)
+2.*real(cmsmu*sym*conj(uncrossed))
+2.*real(cmsmc*sym*conj(crossed))
+2.*real(cmumc*uncrossed*conj(crossed))
;
}
// matrix element for W + Jets with W emission *not* off central qqbar
double ME_W_Cenqqx_qq(
HLV pa, HLV pb, HLV pl,HLV plbar,
std::vector<HLV> partons, bool aqlinepa, bool aqlinepb,
bool qqxmarker, int nabove, int nbelow, bool forwards,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
if (!forwards){ //If Emission from Leg a instead, flip process.
std::swap(pa, pb);
std::reverse(partons.begin(),partons.end());
std::swap(aqlinepa, aqlinepb);
qqxmarker = !qqxmarker;
std::swap(nabove,nbelow);
}
HLV q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
const auto qq = split_into_lightlike(q1);
HLV pq, pqbar;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
// we assume that the W is emitted off a quark line
// if this is not the case, we have to apply CP conjugation,
// which is equivalent to swapping lepton and antilepton momenta
if(aqlinepb) std::swap(pl, plbar);
const HLV p1 = partons.front();
const HLV pn = partons.back();
// helicity labels are for aqlinepa == aqlineb == false,
// In principle the helicities should be flipped
// if either of them is true, but we sum up all helicities,
// so this has no effect in the end.
const double amp_mm = amp_W_Cenqqx_qq<minus, minus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_mp = amp_W_Cenqqx_qq<minus, plus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pm = amp_W_Cenqqx_qq<plus, minus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pp = amp_W_Cenqqx_qq<plus, plus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double t1 = q1.m2();
const double t3 = (q1 - pq - pqbar).m2();
const double amp= WProp(plbar, pl, wprop)*(
amp_mm+amp_mp+amp_pm+amp_pp
)/(9.*4.*t1*t1*t3*t3);
return amp;
}
} // namespace currents
} // namespace HEJ
diff --git a/src/make_RNG.cc b/src/make_RNG.cc
index d96c67f..1315c01 100644
--- a/src/make_RNG.cc
+++ b/src/make_RNG.cc
@@ -1,35 +1,35 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/make_RNG.hh"
#include <algorithm>
#include <locale>
#include "HEJ/exceptions.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/Ranlux64.hh"
namespace HEJ {
- std::unique_ptr<HEJ::RNG> make_RNG(
+ std::unique_ptr<RNG> make_RNG(
std::string const & name,
optional<std::string> const & seed
) {
std::string lname;
std::transform(
begin(name), end(name), std::back_inserter(lname),
[](char c) { return std::tolower(c, std::locale()); }
);
if(lname == "mixmax") {
if(seed) return std::make_unique<Mixmax>(std::stol(*seed));
return std::make_unique<Mixmax>();
}
if(lname == "ranlux64") {
if(seed) return std::make_unique<Ranlux64>(*seed);
return std::make_unique<Ranlux64>();
}
throw std::invalid_argument{"Unknown random number generator: " + name};
}
}
diff --git a/src/resummation_jet.cc b/src/resummation_jet.cc
index 59c772c..f2603f6 100644
--- a/src/resummation_jet.cc
+++ b/src/resummation_jet.cc
@@ -1,114 +1,114 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/resummation_jet.hh"
#include <cassert>
#include <cmath>
#include <cstddef>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include "fastjet/PseudoJet.hh"
#include "HEJ/utility.hh"
-namespace HEJ{
+namespace HEJ {
std::vector<fastjet::PseudoJet> resummation_jet_momenta(
std::vector<fastjet::PseudoJet> const & p_born,
fastjet::PseudoJet const & qperp
) {
// for "new" reshuffling p^B = p + qperp*|p^B|/P^B
double Pperp_born = 0.;
for(auto const & p: p_born) Pperp_born += p.perp();
std::vector<fastjet::PseudoJet> p_res;
p_res.reserve(p_born.size());
for(auto & pB: p_born) {
const double px = pB.px() - qperp.px()*pB.perp()/Pperp_born;
const double py = pB.py() - qperp.py()*pB.perp()/Pperp_born;
const double pperp = std::sqrt(px*px + py*py);
// keep the rapidities fixed
const double pz = pperp*std::sinh(pB.rapidity());
const double E = pperp*std::cosh(pB.rapidity());
p_res.emplace_back(px, py, pz, E);
assert(
- HEJ::nearby_ep(
+ nearby_ep(
p_res.back().rapidity(),
pB.rapidity(),
1e-5
)
);
}
return p_res;
}
namespace {
enum coordinates : std::size_t {
x1, x2
};
namespace ublas = boost::numeric::ublas;
template<class Matrix>
double det(ublas::matrix_expression<Matrix> const& m) {
ublas::permutation_matrix<std::size_t> pivots{m().size1()};
Matrix mLu{m()};
const auto is_singular = lu_factorize(mLu, pivots);
if(is_singular) return 0.;
double det = 1.0;
for (std::size_t i = 0; i < pivots.size(); ++i){
if (pivots(i) != i) det = -det;
det *= mLu(i,i);
}
return det;
}
using ublas::matrix;
}
double resummation_jet_weight(
std::vector<fastjet::PseudoJet> const & p_born,
fastjet::PseudoJet const & qperp
) {
using std::size_t;
static constexpr int num_coordinates = 2;
auto Jacobian = matrix<double>{
num_coordinates*p_born.size(),
num_coordinates*p_born.size()
};
double P_perp = 0.;
for(auto const & J: p_born) P_perp += J.perp();
for(size_t l = 0; l < p_born.size(); ++l){
const double Jl_perp = p_born[l].perp();
for(size_t lp = 0; lp < p_born.size(); ++lp){
const int delta_l = (l == lp);
const auto & Jlp = p_born[lp];
const double Jlp_perp = Jlp.perp();
for(size_t x = x1; x <= x2; ++x){
const double qxy = (x==x1)?qperp.px():qperp.py();
for(size_t xp = x1; xp <= x2; ++xp){
const double Jxy = (xp==x1)?Jlp.px():Jlp.py();
const int delta_x = x == xp;
Jacobian(2*l + x, 2*lp + xp) =
+ delta_l*delta_x
- qxy*Jxy/(P_perp*Jlp_perp) * (delta_l - Jl_perp/P_perp);
}
}
}
}
return det(Jacobian);
}
-}
+} // namespace HEJ

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 3:57 PM (1 d, 15 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805058
Default Alt Text
(288 KB)

Event Timeline