Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7879650
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
33 KB
Subscribers
None
View Options
diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
index a2c05e4..f7071c3 100644
--- a/FixedOrderGen/src/PhaseSpacePoint.cc
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -1,976 +1,976 @@
/**
* \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/PDF.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;
}
size_t count_gluons(PhaseSpacePoint::ConstPartonIterator first,
PhaseSpacePoint::ConstPartonIterator last
){
return std::count_if(first, last, [](HEJ::Particle const & p)
{return p.type == HEJ::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
*/
Subleading possible_qqx(
PhaseSpacePoint::ConstPartonIterator first,
const PhaseSpacePoint::ConstReversePartonIterator& last
){
using namespace subleading;
assert( std::distance( first,last.base() )>2 );
Subleading channels = ALL;
channels.reset(eqqx);
channels.reset(cqqx);
auto const ngluon = count_gluons(first,last.base());
if(ngluon < 2) return channels;
if(first->type==HEJ::pid::gluon || last->type==HEJ::pid::gluon){
channels.set(eqqx);
}
if(std::distance(first,last.base())>=4){
channels.set(cqqx);
}
return channels;
}
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_AWZ_proccess(Process const & proc){
return proc.boson && HEJ::is_AWZ_boson(*proc.boson);
}
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)) );
}
Subleading ensure_AWZ(
double & subl_chance, bool & allow_strange,
HEJ::ParticleID const boson,
PhaseSpacePoint::ConstPartonIterator first_parton,
PhaseSpacePoint::ConstPartonIterator last_parton
){
auto channels = subleading::ALL;
if(std::none_of(first_parton, last_parton,
[&boson](HEJ::Particle const & p){
return can_couple_to_W(p, boson);})) {
// enforce qqx if A/W/Z can't couple somewhere else
// this is ensured to work through filter_partons in reconstruct_incoming
channels.reset(subleading::uno);
assert(channels.any());
subl_chance = 1.;
// strange not allowed for W
if(std::abs(boson)== HEJ::pid::Wp) allow_strange = false;
}
return channels;
}
} // namespace
void PhaseSpacePoint::turn_to_subl(
Subleading const channels,
bool const can_be_uno_backward, bool const can_be_uno_forward,
bool const allow_strange,
HEJ::RNG & ran
){
double const nchannels = channels.count();
double const step = 1./nchannels;
weight_*=nchannels;
unsigned selected = subleading::first;
double rnd = nchannels>1?ran.flat():0.;
// @TODO optimise this sampling
for(; selected<=subleading::last; ++selected){
assert(rnd>=0);
if(channels[selected]){
if(rnd<step) break;
rnd-=step;
}
}
switch(selected){
case subleading::uno:
return turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
case subleading::cqqx:
return turn_to_cqqx(allow_strange, ran);
case subleading::eqqx:
return turn_to_eqqx(allow_strange, ran);
default:
throw std::logic_error{"unreachable"};
}
}
void PhaseSpacePoint::maybe_turn_to_subl(
double chance, Subleading channels, Process const & proc, HEJ::RNG & ran
){
using namespace HEJ;
if(proc.njets <= 2) return;
assert(outgoing_.size() >= 2);
// decide what kind of subleading process is allowed
bool const can_be_uno_backward = uno_possible(cbegin_partons(),
outgoing_.cbegin());
bool const can_be_uno_forward = uno_possible(crbegin_partons(),
outgoing_.crbegin());
if(channels[subleading::uno]){
channels.set(subleading::uno, can_be_uno_backward || can_be_uno_forward);
}
channels &= possible_qqx(cbegin_partons(), crbegin_partons());
bool allow_strange = true;
if(is_AWZ_proccess(proc)) {
channels &= ensure_AWZ(chance, allow_strange, *proc.boson,
cbegin_partons(), cend_partons());
}
std::size_t const nchannels = channels.count();
// no subleading
if(nchannels==0) return;
if(ran.flat() >= chance){
weight_ /= 1 - chance;
return;
}
weight_ /= chance;
// select channel
return turn_to_subl( channels, can_be_uno_backward, can_be_uno_forward,
allow_strange, ran);
}
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) return;
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_qqx_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_cqqx(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 qqx");
auto flavour = select_qqx_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_eqqx(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 qqx");
auto flavour = select_qqx_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 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
if(subl_chance == 0.)
subl_channels.reset();
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
);
// pre fill flavour with gluons
for( auto it = std::make_move_iterator(partons.begin());
it != std::make_move_iterator(partons.end());
++it
){
outgoing_.emplace_back(HEJ::Particle{HEJ::pid::gluon, *it, {}});
}
if(status_ != Status::good) return;
if(proc.boson){ // decay boson
auto const & 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);
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)
);
}
}
// 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_chance, subl_channels, 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;
maybe_turn_to_subl(subl_chance, subl_channels, proc, ran);
if(proc.boson) couple_boson(*proc.boson, 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
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) == HEJ::pid::Wp){
weight_/=M_PI*M_PI*16.; //Corrects B-W factors, see git issue 132
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 {
/// partons are ordered: even = anti, 0 = gluon
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
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;
}
PhaseSpacePoint::part_mask init_allowed(HEJ::ParticleID const id){
if(std::abs(id) == HEJ::pid::proton)
return ~0;
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 ~1; // not a 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, Subleading const subl_channels,
std::array<part_mask,2> allowed_partons,
HEJ::RNG & ran
){
assert(proc.boson);
auto couple_parton = allowed_quarks(*proc.boson);
+ // eqqx possible if one incoming is a gluon
+ if(proc.njets >= 3 && subl_channels[subleading::eqqx]){
+ couple_parton.set(pid_to_index(HEJ::ParticleID::gluon));
+ }
if( // coupling possible through input
allowed_partons[0] == (couple_parton&allowed_partons[0])
|| allowed_partons[1] == (couple_parton&allowed_partons[1])
// cqqx possible
|| (proc.njets >= 4 && subl_channels[subleading::cqqx])
){
return allowed_partons;
}
- // eqqx only possible if one incoming is a gluon
- if(proc.njets >= 3 && subl_channels[subleading::eqqx]){
- couple_parton.set(pid_to_index(HEJ::ParticleID::gluon));
- }
// 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_eqqx(
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 qqx."};
}
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 first_quarks = first_beam;
first_quarks.reset(gluon_idx);
auto & second_beam = allowed_partons[1];
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 qqx => no incoming gluon
* b) 2j => no incoming gluon
* c) >3j => can couple OR is gluon => 2 gluons become qqx later
* 3. pure eqqx 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,
double const subl_chance, Subleading const subl_channels,
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
if(proc.boson && is_AWZ_proccess(proc)){
allowed_partons = incoming_AWZ(proc, subl_channels, allowed_partons, ran);
}
// special case: pure subleading
if(subl_chance!=1.){
return allowed_partons;
}
auto other_channels = subl_channels;
// pure eqqx
other_channels.reset(subleading::eqqx);
if(other_channels.none()){
return incoming_eqqx(allowed_partons, ran);
}
other_channels = subl_channels;
// pure uno
other_channels.reset(subleading::uno);
if(other_channels.none()){
return incoming_uno(allowed_partons, ran);
}
return allowed_partons;
}
void PhaseSpacePoint::reconstruct_incoming(
Process const & proc,
double const subl_chance, Subleading const subl_channels,
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, subl_chance, subl_channels, 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];
}
}
assert(momentum_conserved(1e-7));
}
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"};
}
void PhaseSpacePoint::couple_boson(
HEJ::ParticleID const boson, HEJ::RNG & ran
){
if(std::abs(boson) != HEJ::pid::Wp) return; // only matters for W
/// @TODO this could be use to sanity check gamma and Z
// find all possible quarks
std::vector<PartonIterator> allowed_parts;
for(auto part_it=begin_partons(); part_it!=end_partons(); ++part_it){
// Wp -> up OR anti-down, Wm -> anti-up OR down, no bottom
if ( can_couple_to_W(*part_it, boson) )
allowed_parts.push_back(part_it);
}
if(allowed_parts.empty()){
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<HEJ::ParticleID>( allowed_parts[idx]->type - W_charge );
}
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;
}
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);
}
} // namespace HEJFOG
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:40 PM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806115
Default Alt Text
(33 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment