Page MenuHomeHEPForge
Files F7876992

No OneTemporary
Actions

View Options

diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
index 8b9cdde..a924155 100644
--- a/FixedOrderGen/src/PhaseSpacePoint.cc
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -1,1057 +1,1057 @@
/**
* \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 <cstdint>
#include <cstdlib>
#include <iostream>
#include <iterator>
#include <limits>
#include <tuple>
#include <type_traits>
#include <utility>
#include "Subleading.hh"
#include "Utility.hh"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/PDF.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Particle.hh"
#include "HEJ/RNG.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/kinematics.hh"
#include "HEJ/utility.hh"
#include "JetParameters.hh"
#include "Process.hh"
namespace HEJFOG {
HEJ::Event::EventData to_EventData(PhaseSpacePoint psp){
//! @TODO Same function already in HEJ
HEJ::Event::EventData result;
result.incoming = std::move(psp).incoming_; // NOLINT(bugprone-use-after-move)
result.outgoing = std::move(psp).outgoing_; // NOLINT(bugprone-use-after-move)
// 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{}
)
);
assert(result.outgoing.size() >= 2);
result.decays = std::move(psp).decays_; // NOLINT(bugprone-use-after-move)
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double nan = std::numeric_limits<double>::quiet_NaN();
result.parameters.central = {nan, nan, psp.weight()}; // NOLINT(bugprone-use-after-move)
return result;
}
PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::begin_partons() const {
return cbegin_partons();
}
PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::cbegin_partons() const {
return {HEJ::is_parton, cbegin(outgoing()), cend(outgoing())};
}
PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::end_partons() const {
return cend_partons();
}
PhaseSpacePoint::ConstPartonIterator PhaseSpacePoint::cend_partons() const {
return {HEJ::is_parton, cend(outgoing()), cend(outgoing())};
}
PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::rbegin_partons() const {
return crbegin_partons();
}
PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::crbegin_partons() const {
return std::reverse_iterator<ConstPartonIterator>( cend_partons() );
}
PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::rend_partons() const {
return crend_partons();
}
PhaseSpacePoint::ConstReversePartonIterator PhaseSpacePoint::crend_partons() const {
return std::reverse_iterator<ConstPartonIterator>( cbegin_partons() );
}
PhaseSpacePoint::PartonIterator PhaseSpacePoint::begin_partons() {
return {HEJ::is_parton, begin(outgoing_), end(outgoing_)};
}
PhaseSpacePoint::PartonIterator PhaseSpacePoint::end_partons() {
return {HEJ::is_parton, end(outgoing_), end(outgoing_)};
}
PhaseSpacePoint::ReversePartonIterator PhaseSpacePoint::rbegin_partons() {
return std::reverse_iterator<PartonIterator>( end_partons() );
}
PhaseSpacePoint::ReversePartonIterator PhaseSpacePoint::rend_partons() {
return std::reverse_iterator<PartonIterator>( begin_partons() );
}
namespace {
template<class PartonIt>
bool can_swap_to_uno(PartonIt first_parton){
assert(is_parton(*first_parton));
assert(std::next(first_parton)->type == HEJ::pid::gluon);
return HEJ::is_anyquark(*first_parton);
}
template<class Iterator>
bool allow_strange(
Iterator first_parton, Iterator last_parton,
Process const & proc
) {
const auto boson = proc.boson;
return !(
boson
&& std::abs(*boson)== HEJ::pid::Wp
&& std::none_of(
first_parton, last_parton,
[boson](HEJ::Particle const & p){
return vector_boson_can_couple_to(*boson, p.type);
})
);
}
}
void PhaseSpacePoint::turn_to_subl(
subleading::Channels const channel,
Process const & proc,
HEJ::RNG & ran
) {
using namespace HEJ;
assert(proc.njets > 2);
switch(channel) {
case subleading::unordered: {
bool const can_be_uno_backward = can_swap_to_uno(outgoing().cbegin());
bool const can_be_uno_forward = can_swap_to_uno(outgoing().crbegin());
return turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
}
case subleading::central_qqbar: {
return turn_to_cqqbar(proc, ran);
}
case subleading::extremal_qqbar: {
return turn_to_eqqbar(proc, ran);
}
default:;
}
throw std::logic_error{"enum value not covered"};
}
void PhaseSpacePoint::turn_to_uno(
const bool can_be_uno_backward, const bool can_be_uno_forward,
HEJ::RNG & ran
){
if(can_be_uno_backward && can_be_uno_forward){
weight_ *= 2.;
if(ran.flat() < 0.5){
return std::swap(begin_partons()->type, std::next(begin_partons())->type);
}
return std::swap(rbegin_partons()->type, std::next(rbegin_partons())->type);
}
if(can_be_uno_backward){
return std::swap(begin_partons()->type, std::next(begin_partons())->type);
}
assert(can_be_uno_forward);
std::swap(rbegin_partons()->type, std::next(rbegin_partons())->type);
}
//! select flavour of quark
HEJ::ParticleID PhaseSpacePoint::select_qqbar_flavour(
Process const & proc, HEJ::RNG & ran
){
bool const strange_allowed = allow_strange(
outgoing().begin(),
outgoing().end(),
proc
);
const double r1 = 2.*ran.flat()-1.;
const double max_flavour = strange_allowed?HEJ::N_F:HEJ::N_F-1;
weight_ *= max_flavour*2;
double const flavour = HEJ::pid::down + std::floor(std::abs(r1)*max_flavour);
return static_cast<HEJ::ParticleID>(flavour*(r1<0.?-1:1));
}
void PhaseSpacePoint::turn_to_cqqbar(Process const & proc, HEJ::RNG & ran){
// we assume all FKL partons to be gluons
auto first = std::next(begin_partons());
auto last = std::next(rbegin_partons());
auto const ng = std::distance(first, last.base());
if(ng < 2)
throw std::logic_error("not enough gluons to create qqbar");
auto flavour = select_qqbar_flavour(proc, ran);
// select gluon for switch
if(ng!=2){
const double steps = 1./(ng-1.);
weight_ /= steps;
for(auto rnd = ran.flat(); rnd>steps; ++first){
rnd-=steps;
}
}
first->type = flavour;
std::next(first)->type = anti(flavour);
}
void PhaseSpacePoint::turn_to_eqqbar(Process const & proc, HEJ::RNG & ran){
/// find first and last gluon in FKL chain
auto first = begin_partons();
const bool can_forward = !is_anyquark(*first);
auto last = rbegin_partons();
const bool can_backward = !is_anyquark(*last);
if(std::distance(first, last.base()) < 2)
throw std::logic_error("not enough gluons to create qqbar");
auto flavour = select_qqbar_flavour(proc, ran);
// select gluon for switch
if(can_forward && !can_backward){
first->type = flavour;
std::next(first)->type = anti(flavour);
return;
}
if(!can_forward && can_backward){
last->type = flavour;
std::next(last)->type = anti(flavour);
return;
}
assert(can_forward && can_backward);
weight_*=2.;
if(ran.flat()>0.5){
first->type = flavour;
std::next(first)->type = anti(flavour);
return;
}
last->type = flavour;
std::next(last)->type = anti(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};
}
Decay PhaseSpacePoint::select_decay_channel(
std::vector<Decay> const & decays,
HEJ::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"};
}
namespace {
//! generate decay products of a boson
std::vector<HEJ::Particle> decay_boson(
HEJ::Particle const & parent,
std::vector<HEJ::ParticleID> const & decays,
HEJ::RNG & ran
){
if(decays.size() != 2){
throw HEJ::not_implemented{
"only decays into two particles are implemented"
};
}
std::vector<HEJ::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.; // Jacobian Factors for W in line 418
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;
}
} // namespace
std::vector<HEJ::Particle> PhaseSpacePoint::decay_channel(
HEJ::Particle const & parent,
std::vector<Decay> const & decays,
HEJ::RNG & ran
){
const auto channel = select_decay_channel(decays, ran);
return decay_boson(parent, channel.products, ran);
}
namespace {
//! adds a particle to target (in correct rapidity ordering)
//! @returns positon of insertion
auto insert_particle(std::vector<HEJ::Particle> & target,
HEJ::Particle && particle
){
const auto pos = std::upper_bound(
begin(target),end(target),particle,HEJ::rapidity_less{}
);
target.insert(pos, std::move(particle));
return pos;
}
} // namespace
PhaseSpacePoint::PhaseSpacePoint(
Process const & proc,
JetParameters const & jet_param,
HEJ::PDF & pdf, double E_beam,
double const subl_chance,
Subleading const subl_channels,
ParticlesDecayMap const & particle_decays,
HEJ::EWConstants const & ew_parameters,
HEJ::RNG & ran
){
assert(proc.njets >= 2);
status_ = Status::good;
weight_ = 1;
// ensure that all setting are consistent
// a more thorough check is in EventGenerator
assert(subl_chance > 0. || subl_channels.none());
const std::size_t 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
);
if(status_ != Status::good) return;
// pre fill flavour with gluons
for(auto && parton: std::move(partons)){
outgoing_.emplace_back(HEJ::Particle{HEJ::pid::gluon, std::move(parton), {}});
}
const std::optional<subleading::Channels> channel = select_channel(
subl_channels,
subl_chance,
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, channel, pdf, E_beam, jet_param.min_pt, ran);
if(status_ != Status::good) return;
// set outgoing states
begin_partons()->type = incoming_[0].type;
rbegin_partons()->type = incoming_[1].type;
if(channel) turn_to_subl(*channel, proc, ran);
if(proc.boson) {
add_boson(proc, particle_decays, ew_parameters, ran);
// total momentum has changed, update & check incoming momenta
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
const auto [xa, xb] = xa_xb(E_beam);
if (xa>1. || xb>1.){
weight_=0;
status_ = Status::too_much_energy;
return;
}
}
assert(weight_ != 0 || status() == Status::good);
assert(momentum_conserved(1e-7));
}
void PhaseSpacePoint::add_boson(
Process const & proc,
ParticlesDecayMap const & particle_decays,
HEJ::EWConstants const & ew_parameters,
HEJ::RNG & ran
) {
auto const & boson_prop = (*proc.boson == HEJ::pid::Z_photon_mix
? ew_parameters.prop(HEJ::pid::Z)
: 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);
auto const & 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_channel(outgoing_[boson_idx], boson_decay->second, ran)
);
}
}
// pt generation, see eq:pt_sampling in developer manual
double PhaseSpacePoint::gen_hard_pt(
const int np , const double ptmin, const double ptmax, const double /* y */,
HEJ::RNG & ran
){
// heuristic parameter for pt sampling, see eq:pt_par in developer manual
const double ptpar = ptmin + np/5.;
const double arctan = std::atan((ptmax - ptmin)/ptpar);
const double xpt = ran.flat();
const double pt = ptmin + ptpar*std::tan(xpt*arctan);
const double cosine = std::cos(xpt*arctan);
weight_ *= pt*ptpar*arctan/(cosine*cosine);
return pt;
}
// pt generation below minimum resummaton jet pt,
// see eq:FO_pt_sampling in developer manual
double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, HEJ::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,
HEJ::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_ = 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,
HEJ::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_ = Status::not_enough_jets;
return {};
}
std::sort(begin(partons), end(partons), HEJ::rapidity_less{});
return partons;
}
HEJ::Particle PhaseSpacePoint::gen_boson(
HEJ::ParticleID bosonid, double mass, double width,
HEJ::RNG & ran
){
// Usual phase space measure for single particle,
// see eq:single_particle_phase_space in developer manual
weight_ /= 16.*std::pow(M_PI, 3);
const double y = (bosonid == HEJ::pid::Higgs)?
gen_Higgs_rapidity(ran):
gen_V_boson_rapidity(bosonid, ran);
// off-shell sampling, see eq:breit-wigner-sampling in developer manual
const double r1 = ran.flat();
const double s_boson = mass*(
mass + width*std::tan(M_PI/2.*r1 + (r1-1.)*std::atan(mass/width))
);
if(bosonid != HEJ::pid::Higgs){
// two-particle phase space,
// compare eq:two_particle_phase_space and eq:one_particle_phase_space
// in developer manual
weight_/=M_PI*M_PI*16.;
// Jacobian, see eq:breit-wigner-sampling_jacobian in developer manual
weight_*= mass*width*( M_PI+2.*std::atan(mass/width) )
/ ( 1. + std::cos( M_PI*r1 + 2.*(r1-1.)*std::atan(mass/width) ) );
}
return {
bosonid,
gen_last_momentum(outgoing_, s_boson, y),
{}
};
}
namespace {
bool is_uno_forward(
std::array<HEJ::Particle, 2> const & incoming,
std::vector<HEJ::Particle> const & partons
) {
if(partons.size() < 3) return false;
assert(
std::is_sorted(
begin(incoming), end(incoming),
HEJ::rapidity_less{}
)
);
assert(
std::is_sorted(
begin(partons), end(partons),
HEJ::rapidity_less{}
)
);
assert(is_parton(*partons.rbegin()));
assert(is_parton(*std::next(partons.rbegin())));
return
HEJ::is_anyquark(incoming.back())
&& partons.rbegin()->type == HEJ::pid::gluon
&& HEJ::is_anyquark(*std::next(partons.rbegin()));
}
bool is_uno_backward(
std::array<HEJ::Particle, 2> const & incoming,
std::vector<HEJ::Particle> const & partons
) {
if(partons.size() < 3) return false;
assert(
std::is_sorted(
begin(incoming), end(incoming),
HEJ::rapidity_less{}
)
);
assert(
std::is_sorted(
begin(partons), end(partons),
HEJ::rapidity_less{}
)
);
assert(is_parton(*partons.begin()));
assert(is_parton(*std::next(partons.begin())));
return
HEJ::is_anyquark(incoming.front())
&& partons.begin()->type == HEJ::pid::gluon
&& HEJ::is_anyquark(*std::next(partons.begin()));
}
}
double PhaseSpacePoint::gen_Higgs_rapidity(HEJ::RNG & ran) {
constexpr double STDDEV_Y_H = 1.6;
constexpr double Y_H_UNO_OFFSET = 1.6;
if(is_uno_forward(incoming(), outgoing())) {
return outgoing_[outgoing().size() - 2].rapidity()
- Y_H_UNO_OFFSET
- random_normal_trunc(
STDDEV_Y_H,
ran,
-Y_H_UNO_OFFSET,
std::numeric_limits<double>::max()
);
} else if(is_uno_backward(incoming(), outgoing())) {
return outgoing_[1].rapidity()
+ Y_H_UNO_OFFSET
+ random_normal_trunc(
STDDEV_Y_H,
ran,
-Y_H_UNO_OFFSET,
std::numeric_limits<double>::max()
);
}
return random_normal(STDDEV_Y_H, ran);
}
namespace {
template<class Iterator, class Predicate>
Iterator find_if_nth(
Iterator begin, Iterator end, Predicate p, std::size_t n
) {
return std::find_if(
begin, end,
[p,n](auto const & i) mutable {
if(!p(i)) return false;
if(n > 0) {
--n;
return false;
} else {
return true;
}
}
);
}
}
double PhaseSpacePoint::gen_V_boson_rapidity(
HEJ::ParticleID const V_boson,
HEJ::RNG & ran
){
constexpr double STDDEV_Y_V = 1.2;
const auto can_couple = [V_boson](HEJ::Particle const & p) {
return vector_boson_can_couple_to(V_boson, p.type);
};
const std::size_t n_pot_emitters = std::count_if(
outgoing().begin(), outgoing().end(),
can_couple
);
assert(n_pot_emitters > 0);
const std::size_t emitter_idx = ran.flat() * n_pot_emitters;
auto emitter = find_if_nth(
outgoing_.begin(), outgoing_.end(),
can_couple,
emitter_idx
);
assert(emitter != outgoing_.end());
// For a W we get a flavour change in the emitter,
// so we have to _sum_ over all possible emitters
// and have to adjust the weight
// since we only pick a single term in that sum here.
// For a Z/photon we do not get distinct configurations,
// so there is no sum and no weight adjustment.
if(std::abs(V_boson) == HEJ::pid::Wp) {
emitter->type = static_cast<HEJ::ParticleID>(
emitter->type - HEJ::charge(V_boson).numerator()
);
weight_ *= n_pot_emitters;
}
const double dy = random_normal(STDDEV_Y_V, ran);
return emitter->rapidity() + dy;
}
namespace {
/// partons are ordered: even = anti, 0 = gluon
constexpr HEJ::ParticleID index_to_pid(size_t i){
if(!i) return HEJ::pid::gluon;
return static_cast<HEJ::ParticleID>( i%2 ? (i+1)/2 : -i/2 );
}
/// partons are ordered: even = anti, 0 = gluon
constexpr size_t pid_to_index(HEJ::ParticleID id){
if(id==HEJ::pid::gluon) return 0;
return id>0 ? id*2-1 : std::abs(id)*2;
}
constexpr PhaseSpacePoint::part_mask GLUON = 1u << pid_to_index(HEJ::pid::gluon);
constexpr PhaseSpacePoint::part_mask ALL = ~0;
PhaseSpacePoint::part_mask init_allowed(HEJ::ParticleID const id){
if(std::abs(id) == HEJ::pid::proton)
return ALL;
PhaseSpacePoint::part_mask out{0};
if(HEJ::is_parton(id))
out[pid_to_index(id)] = true;
return out;
}
/// decides which "index" (see index_to_pid) are allowed for process
PhaseSpacePoint::part_mask allowed_quarks(HEJ::ParticleID const boson){
if(std::abs(boson) != HEJ::pid::Wp){
return ~GLUON;
}
// special case W:
// Wp: anti-down or up-type quark, no b/t
// Wm: down or anti-up-type quark, no b/t
return boson>0?0b00011001100 // NOLINT(readability-magic-numbers)
:0b00100110010; // NOLINT(readability-magic-numbers)
}
} // namespace
std::array<PhaseSpacePoint::part_mask,2> PhaseSpacePoint::incoming_AWZ(
Process const & proc,
std::array<part_mask,2> allowed_partons,
HEJ::RNG & ran
){
assert(proc.boson);
auto couple_parton = allowed_quarks(*proc.boson);
if( // coupling possible through input
allowed_partons[0] == (couple_parton&allowed_partons[0])
|| allowed_partons[1] == (couple_parton&allowed_partons[1])
){
return allowed_partons;
}
// only first can couple
if( (allowed_partons[0]&couple_parton).any()
&&(allowed_partons[1]&couple_parton).none()
){
return {allowed_partons[0]&couple_parton, allowed_partons[1]};
}
// only second can couple
if( (allowed_partons[0]&couple_parton).none()
&& (allowed_partons[1]&couple_parton).any()
){
return {allowed_partons[0], allowed_partons[1]&couple_parton};
}
// both can couple
if( (allowed_partons[0]&couple_parton).any()
&& (allowed_partons[1]&couple_parton).any()
){
double rnd = ran.flat();
weight_*=3.;
if(rnd<1./3.){
return {
allowed_partons[0] & couple_parton,
allowed_partons[1] & ~couple_parton
};
}
if(rnd<2./3.){
return {
allowed_partons[0] & ~couple_parton,
allowed_partons[1] & couple_parton
};
}
return {
allowed_partons[0] & couple_parton,
allowed_partons[1] & couple_parton
};
}
throw std::invalid_argument{"Incoming state not allowed."};
}
std::array<PhaseSpacePoint::part_mask,2> PhaseSpacePoint::incoming_eqqbar(
std::array<part_mask,2> allowed_partons, HEJ::RNG & ran
){
auto const gluon_idx = pid_to_index(HEJ::pid::gluon);
auto & first_beam = allowed_partons[0];
auto & second_beam = allowed_partons[1];
if(first_beam[gluon_idx] && !second_beam[gluon_idx]){
first_beam.reset();
first_beam.set(gluon_idx);
return allowed_partons;
}
if(!first_beam[gluon_idx] && second_beam[gluon_idx]) {
second_beam.reset();
second_beam.set(gluon_idx);
return allowed_partons;
}
if(first_beam[gluon_idx] && second_beam[gluon_idx]) {
// both beams can be gluons
// if one can't be a quark everything is good
auto first_quarks = first_beam;
first_quarks.reset(gluon_idx);
auto second_quarks = second_beam;
second_quarks.reset(gluon_idx);
if(first_quarks.none() || second_quarks.none()){
return allowed_partons;
}
// else choose one to be a gluon
double rnd = ran.flat();
weight_*=3.;
if(rnd<1./3.){
allowed_partons[0].reset();
allowed_partons[0].set(gluon_idx);
allowed_partons[1].reset(gluon_idx);
} else if(rnd<2./3.){
allowed_partons[1].reset();
allowed_partons[1].set(gluon_idx);
allowed_partons[0].reset(gluon_idx);
} else {
allowed_partons[0].reset();
allowed_partons[0].set(gluon_idx);
allowed_partons[1].reset();
allowed_partons[1].set(gluon_idx);
}
return allowed_partons;
}
throw std::invalid_argument{
"Incoming state not allowed for pure extremal qqbar."};
}
std::array<PhaseSpacePoint::part_mask,2> PhaseSpacePoint::incoming_uno(
std::array<part_mask,2> allowed_partons, HEJ::RNG & ran
){
auto const gluon_idx = pid_to_index(HEJ::pid::gluon);
auto & first_beam = allowed_partons[0];
auto & second_beam = allowed_partons[1];
auto first_quarks = first_beam;
first_quarks.reset(gluon_idx);
auto second_quarks = second_beam;
second_quarks.reset(gluon_idx);
if(first_quarks.any() && second_quarks.none()){
first_beam.reset(gluon_idx);
return allowed_partons;
}
if(first_quarks.none() && second_quarks.any()) {
second_beam.reset(gluon_idx);
return allowed_partons;
}
if(first_quarks.any() && second_quarks.any()) {
// both beams can be quarks
// if one can't be gluon everything is good
if(!first_beam[gluon_idx] || !second_beam[gluon_idx]){
return allowed_partons;
}
// else choose one to be a quark
double rnd = ran.flat();
weight_*=3.;
if(rnd<1./3.){
allowed_partons[0].reset(gluon_idx);
allowed_partons[1].reset();
allowed_partons[1].set(gluon_idx);
} else if(rnd<2./3.){
allowed_partons[1].reset(gluon_idx);
allowed_partons[0].reset();
allowed_partons[0].set(gluon_idx);
} else {
allowed_partons[0].reset(gluon_idx);
allowed_partons[1].reset(gluon_idx);
}
return allowed_partons;
}
throw std::invalid_argument{
"Incoming state not allowed for pure unordered."};
}
/**
* @brief Returns list of all allowed initial states partons
*
* 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 qqbar => no incoming gluon
* b) 2j => no incoming gluon
* c) >3j => can couple OR is gluon => 2 gluons become qqbar later
* 3. pure eqqbar requires at least one gluon
* 4. pure uno requires at least one quark
*/
std::array<PhaseSpacePoint::part_mask,2> PhaseSpacePoint::allowed_incoming(
Process const & proc,
std::optional<subleading::Channels> channel,
HEJ::RNG & ran
){
// all possible incoming states
std::array<part_mask,2> allowed_partons{
init_allowed(proc.incoming[0]),
init_allowed(proc.incoming[1])
};
// special case: A/W/Z without qqbar
if(
is_AWZ_process(proc)
&& (!channel || *channel == subleading::unordered)
) {
return incoming_AWZ(proc, allowed_partons, ran);
}
if(!channel) return allowed_partons;
switch(*channel) {
case subleading::central_qqbar:
return allowed_partons;
case subleading::extremal_qqbar:
return incoming_eqqbar(allowed_partons, ran);
case subleading::unordered:
return incoming_uno(allowed_partons, ran);
default:;
};
throw std::logic_error{"enum value not covered"};
}
std::pair<double, double> PhaseSpacePoint::xa_xb(
double const E_beam
) const {
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;
return std::make_pair(xa, xb);
}
void PhaseSpacePoint::reconstruct_incoming(
Process const & proc,
std::optional<subleading::Channels> const channel,
HEJ::PDF & pdf, double E_beam,
double uf,
HEJ::RNG & ran
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
const auto [xa, xb] = xa_xb(E_beam);
// abort if phase space point is outside of collider energy reach
if(xa > 1. || xb > 1.) {
weight_=0;
status_ = Status::too_much_energy;
return;
}
auto const & ids = proc.incoming;
std::array<part_mask,2> allowed_partons
= allowed_incoming(proc, channel, ran);
for(size_t i = 0; i < 2; ++i){
if(ids[i] == HEJ::pid::proton || ids[i] == HEJ::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];
}
}
}
HEJ::ParticleID PhaseSpacePoint::generate_incoming_id(
size_t const beam_idx, double const x, double const uf,
HEJ::PDF & pdf, part_mask allowed_partons, HEJ::RNG & ran
){
std::array<double,NUM_PARTONS> pdf_wt{};
pdf_wt[0] = allowed_partons[0]?
std::fabs(pdf.pdfpt(beam_idx,x,uf,HEJ::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;
+ HEJ::C_F/HEJ::C_A*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"};
}
double PhaseSpacePoint::random_normal( double stddev, HEJ::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;
}
double PhaseSpacePoint::random_normal_trunc(
const double stddev,
HEJ::RNG & ran,
const double min,
const double max
) {
assert(min <= max);
double result;
std::uint64_t trials = 0;
do {
++trials;
const double r1 = ran.flat();
const double r2 = ran.flat();
const double lninvr1 = -std::log(r1);
result = stddev*std::sqrt(2.*lninvr1)*std::cos(2.*M_PI*r2);
} while(result < min || result > max);
weight_ *= exp(result*result/(2*stddev*stddev))*std::sqrt(2.*M_PI)*stddev / trials;
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 HEJ::nearby_ep(diff, fastjet::PseudoJet{}, ep);
}
std::optional<subleading::Channels> PhaseSpacePoint::select_channel(
Subleading const subl_channels,
double const chance,
HEJ::RNG & ran
) {
const double r = ran.flat();
if(r > chance) {
weight_ /= 1 - chance;
return {};
}
// pick a random subleading channel (flat)
const std::size_t num_channels = subl_channels.count();
weight_ *= num_channels / chance;
std::size_t channel_idx = r / chance * num_channels;
for(
unsigned channel = subleading::first;
channel <= subleading::last;
++channel
) {
if(subl_channels.test(channel)) {
if(channel_idx == 0) {
return {static_cast<subleading::Channels>(channel)};
} else {
--channel_idx;
}
}
}
throw std::logic_error{"unreachable"};
}
} // namespace HEJFOG

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 2:49 PM (1 d, 13 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804815
Default Alt Text
(34 KB)

Event Timeline