Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7878609
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
46 KB
Subscribers
None
View Options
diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh
index 1b08985..0a19159 100644
--- a/FixedOrderGen/include/PhaseSpacePoint.hh
+++ b/FixedOrderGen/include/PhaseSpacePoint.hh
@@ -1,323 +1,319 @@
/** \file PhaseSpacePoint.hh
* \brief Contains the PhaseSpacePoint Class
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <bitset>
#include <cstddef>
#include <unordered_map>
#include <optional>
#include <vector>
#include "boost/iterator/filter_iterator.hpp"
#include "HEJ/Event.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Particle.hh"
#include "fastjet/PseudoJet.hh"
#include "Decay.hh"
#include "Status.hh"
#include "Subleading.hh"
namespace HEJ {
class EWConstants;
class PDF;
class RNG;
}
namespace HEJFOG {
struct JetParameters;
struct Process;
//! A point in resummation phase space
class PhaseSpacePoint{
public:
//! Iterator over partons
using ConstPartonIterator = boost::filter_iterator<
bool (*)(HEJ::Particle const &),
std::vector<HEJ::Particle>::const_iterator
>;
//! Reverse Iterator over partons
using ConstReversePartonIterator = std::reverse_iterator<
ConstPartonIterator>;
//! Default PhaseSpacePoint Constructor
PhaseSpacePoint() = delete;
//! PhaseSpacePoint Constructor
/**
* @param proc The process to generate
* @param jet_param Jet defintion & cuts
* @param pdf The pdf set (used for sampling)
* @param E_beam Energie of the beam
* @param subl_chance Chance to turn a potentially unordered
* emission into an actual one
* @param subl_channels Possible subleading channels.
* see HEJFOG::Subleading
* @param particle_properties Properties of producted boson
*
* Initially, only FKL phase space points are generated. subl_chance gives
* the chance of turning one emissions into a subleading configuration,
* i.e. either unordered or central quark/anti-quark pair. Unordered
* emissions require that the most extremal emission in any direction is
* a quark or anti-quark and the next emission is a gluon. Quark/anti-quark
* pairs are only generated for W processes. At most one subleading
* emission will be generated in this way.
*/
PhaseSpacePoint(
Process const & proc,
JetParameters const & jet_param,
HEJ::PDF & pdf, double E_beam,
double subl_chance,
Subleading subl_channels,
ParticlesDecayMap const & particle_decays,
HEJ::EWConstants const & ew_parameters,
HEJ::RNG & ran
);
//! Get Weight Function
/**
* @returns Weight of Event
*/
double weight() const{
return weight_;
}
Status status() const{
return status_;
}
//! Get Incoming Function
/**
* @returns Incoming Particles
*/
std::array<HEJ::Particle, 2> const & incoming() const{
return incoming_;
}
//! Get Outgoing Function
/**
* @returns Outgoing Particles
*/
std::vector<HEJ::Particle> const & outgoing() const{
return outgoing_;
}
std::unordered_map<std::size_t, std::vector<HEJ::Particle>> const & decays() const{
return decays_;
}
//! 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;
private:
static constexpr std::size_t NUM_PARTONS = 11; //!< Number of partons
public:
using part_mask = std::bitset<NUM_PARTONS>; //!< Selection mask for partons
private:
friend HEJ::Event::EventData to_EventData(PhaseSpacePoint psp);
/**
* @internal
* @brief Generate LO parton momentum
*
* @param np Number of partons to generate
* @param is_pure_jets If true ensures momentum conservation in x and y
* @param jet_param Jet properties to fulfil
* @param max_pt max allowed pt for a parton (typically E_CMS)
* @param ran Random Number Generator
*
* @returns Momentum of partons
*
* Ensures that each parton is in its own jet.
* Generation is independent of parton flavour. Output is sorted in rapidity.
*/
std::vector<fastjet::PseudoJet> gen_LO_partons(
int np, bool is_pure_jets,
JetParameters const & jet_param,
double max_pt,
HEJ::RNG & ran
);
void add_boson(
Process const & proc,
ParticlesDecayMap const & particle_decays,
HEJ::EWConstants const & ew_parameters,
HEJ::RNG & ran
);
HEJ::Particle gen_boson(
HEJ::ParticleID bosonid, double mass, double width,
HEJ::RNG & ran
);
double gen_V_boson_rapidity(
HEJ::ParticleID V_boson,
HEJ::RNG & ran
);
double gen_Higgs_rapidity(HEJ::RNG & ran);
template<class ParticleMomenta>
fastjet::PseudoJet gen_last_momentum(
ParticleMomenta const & other_momenta,
double mass_square, double y
) const;
bool jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const;
/**
* @internal
* @brief Generate incoming partons according to the PDF
*
* @param uf Scale used in the PDF
*/
void reconstruct_incoming(
Process const & proc,
std::optional<subleading::Channels> channel,
HEJ::PDF & pdf, double E_beam,
double uf,
HEJ::RNG & ran
);
HEJ::ParticleID generate_incoming_id(
std::size_t beam_idx, double x, double uf, HEJ::PDF & pdf,
part_mask allowed_partons, HEJ::RNG & ran
);
//! @brief Returns list of all allowed initial states partons
std::array<part_mask,2> allowed_incoming(
Process const & proc,
std::optional<subleading::Channels> channel,
HEJ::RNG & ran
);
//! Ensure that incoming state compatible with A/W/Z production
//! Should not be called for final-state qqbar configurations,
//! since A/W/Z emission doesn't constrain the initial state there
std::array<part_mask,2> incoming_AWZ(
Process const & proc,
std::array<part_mask,2> allowed_partons,
HEJ::RNG & ran
);
//! Ensure that incoming state compatible with extremal qqbar
std::array<part_mask,2> incoming_eqqbar(
std::array<part_mask,2> allowed_partons, HEJ::RNG & ran);
//! Ensure that incoming state compatible with unordered
std::array<part_mask,2> incoming_uno(
std::array<part_mask,2> allowed_partons, HEJ::RNG & ran);
bool momentum_conserved(double ep) const;
bool extremal_FKL_ok(
std::vector<fastjet::PseudoJet> const & partons
) const;
double random_normal(double stddev, HEJ::RNG & ran);
double random_normal_trunc(
double stddev,
HEJ::RNG & ran,
double min,
double max
);
/**
* @internal
* @brief Turns a FKL configuration into a subleading one
*
* @param chance Chance to switch to subleading configuration
* @param channels Allowed channels for subleading process
* @param proc Process to decide which subleading
* configurations are allowed
*
* With a chance of "chance" the FKL configuration is either turned into
* a unordered configuration or, for A/W/Z bosons, a configuration with
* a central quark/anti-quark pair.
*/
void turn_to_subl(
subleading::Channels channel,
Process const & proc,
HEJ::RNG & ran
);
- void turn_to_subl(
- Subleading channels,
- bool can_be_uno_backward, bool can_be_uno_forward, bool allow_strange,
- HEJ::RNG & ran);
void turn_to_uno(bool can_be_uno_backward, bool can_be_uno_forward,
HEJ::RNG & ran);
HEJ::ParticleID select_qqbar_flavour(bool allow_strange, HEJ::RNG & ran);
void turn_to_cqqbar(bool allow_strange, HEJ::RNG & ran);
void turn_to_eqqbar(bool allow_strange, HEJ::RNG & ran);
//! decay where we select the decay channel
std::vector<HEJ::Particle> decay_channel(
HEJ::Particle const & parent,
std::vector<Decay> const & decays,
HEJ::RNG & ran
);
Decay select_decay_channel(
std::vector<Decay> const & decays,
HEJ::RNG & ran
);
double gen_hard_pt(
int np, double ptmin, double ptmax, double y,
HEJ::RNG & ran
);
double gen_soft_pt(int np, double max_pt, HEJ::RNG & ran);
double gen_parton_pt(
int count, JetParameters const & jet_param, double max_pt, double y,
HEJ::RNG & ran
);
//! Iterator over partons (non-const)
using PartonIterator = boost::filter_iterator<
bool (*)(HEJ::Particle const &),
std::vector<HEJ::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::optional<subleading::Channels> select_channel(
Subleading subl_channels,
double chance,
HEJ::RNG & ran
);
double weight_;
Status status_;
std::array<HEJ::Particle, 2> incoming_;
std::vector<HEJ::Particle> outgoing_;
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<std::size_t, std::vector<HEJ::Particle>> decays_;
};
//! Extract HEJ::Event::EventData from PhaseSpacePoint
HEJ::Event::EventData to_EventData(PhaseSpacePoint psp);
} // namespace HEJFOG
diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
index 098fa91..ac13f69 100644
--- a/FixedOrderGen/src/PhaseSpacePoint.cc
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -1,1081 +1,1064 @@
/**
* \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 {
bool can_swap_to_uno(
HEJ::Particle const & p1, HEJ::Particle const & p2
){
assert(is_parton(p1) && is_parton(p2));
return p1.type != HEJ::pid::gluon
&& p2.type == HEJ::pid::gluon;
}
template<class PartonIt, class OutIt>
bool uno_possible(
PartonIt first_parton, OutIt first_out
){
using namespace HEJ;
// Special case: Higgs can not be outside of uno
if(first_out->type == pid::Higgs
|| std::next(first_out)->type==pid::Higgs){
return false;
}
// decide what kind of subleading process is allowed
return can_swap_to_uno( *first_parton, *std::next(first_parton) );
}
- bool is_up_type(HEJ::Particle const & part){
- return is_anyquark(part) && ((std::abs(part.type)%2) == 0);
- }
- bool is_down_type(HEJ::Particle const & part){
- return is_anyquark(part) && ((std::abs(part.type)%2) != 0);
- }
- bool can_couple_to_W(
- HEJ::Particle const & part, HEJ::pid::ParticleID const W_id
- ){
- const int W_charge = W_id>0?1:-1;
- return std::abs(part.type)<HEJ::pid::b
- && ( (W_charge*part.type > 0 && is_up_type(part))
- || (W_charge*part.type < 0 && is_down_type(part)) );
- }
- } // namespace
-
- namespace {
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 can_couple_to_W(p, *boson);
+ 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 = uno_possible(cbegin_partons(),
outgoing_.cbegin());
bool const can_be_uno_forward = uno_possible(crbegin_partons(),
outgoing_.crbegin());
return turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
}
case subleading::central_qqbar: {
bool const strange_allowed = allow_strange(
begin_partons(),
end_partons(),
proc
);
return turn_to_cqqbar(strange_allowed, ran);
}
case subleading::extremal_qqbar: {
bool const strange_allowed = allow_strange(
begin_partons(),
end_partons(),
proc
);
return turn_to_eqqbar(strange_allowed, 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(
const bool allow_strange, HEJ::RNG & ran
){
const double r1 = 2.*ran.flat()-1.;
const double max_flavour = allow_strange?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(const bool allow_strange, 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(allow_strange, 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(const bool allow_strange, 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(allow_strange, 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 incoming momenta
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
if(incoming().front().E() + incoming().back().E() > 2 * E_beam) {
weight_ = 0;
status_ = Status::too_much_energy;
return;
}
}
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;
}
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;
if(is_uno_forward(incoming(), outgoing())) {
return random_normal_trunc(
STDDEV_Y_H,
ran,
std::numeric_limits<double>::lowest(),
outgoing_[outgoing().size() - 2].rapidity()
);
} else if(is_uno_backward(incoming(), outgoing())) {
return random_normal_trunc(
STDDEV_Y_H,
ran,
outgoing_[1].rapidity(),
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.8;
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_proccess(proc)
+ 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"};
}
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_);
// 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_ = 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;
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
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:26 PM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3789616
Default Alt Text
(46 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment