Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7879835
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
35 KB
Subscribers
None
View Options
diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc
index 32bb082..a52dd9e 100644
--- a/src/PhaseSpacePoint.cc
+++ b/src/PhaseSpacePoint.cc
@@ -1,984 +1,986 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/PhaseSpacePoint.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdlib>
#include <functional>
#include <iterator>
#include <limits>
#include <numeric>
#include <random>
#include <tuple>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "HEJ/Constants.hh"
#include "HEJ/Event.hh"
#include "HEJ/JetSplitter.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/RNG.hh"
#include "HEJ/event_types.hh"
#include "HEJ/kinematics.hh"
#include "HEJ/resummation_jet.hh"
#include "HEJ/utility.hh"
namespace HEJ {
namespace {
constexpr int MAX_JET_USER_IDX = PhaseSpacePoint::NG_MAX;
bool is_nonjet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() > MAX_JET_USER_IDX;
}
bool is_jet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() <= MAX_JET_USER_IDX;
}
namespace user_idx {
//! user indices for partons with extremal rapidity
enum ID: int {
qqbar_mid1 = -9,
qqbar_mid2 = -8,
qqbarb = -7,
qqbarf = -6,
unob = -5,
unof = -4,
backward_fkl = -3,
forward_fkl = -2,
};
} // namespace user_idx
using UID = user_idx::ID;
double phase_space_normalisation(
int num_Born_jets, int num_out_partons
){
return std::pow(16.*std::pow(M_PI,3), num_Born_jets - num_out_partons);
}
} // namespace
Event::EventData to_EventData(PhaseSpacePoint psp){
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),
rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double nan = std::numeric_limits<double>::quiet_NaN();
result.decays = std::move(psp).decays_; // NOLINT(bugprone-use-after-move)
result.parameters.central = {nan, nan, psp.weight()}; // NOLINT(bugprone-use-after-move)
return result;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
return sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
}
bool PhaseSpacePoint::pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const{
return cluster_jets(jets).size() == jets.size();
}
namespace {
// find iterators to central qqbar emission
auto get_central_qqbar(Event const & ev) {
// find born quarks (ignore extremal partons)
auto const firstquark = std::find_if(
std::next(ev.begin_partons()), std::prev(ev.end_partons(), 2),
[](Particle const & s){ return (is_anyquark(s)); }
);
// assert that it is a q-q_bar pair.
assert(std::distance(firstquark, ev.end_partons()) != 2);
assert(
( is_quark(*firstquark) && is_antiquark(*std::next(firstquark)) )
|| ( is_antiquark(*firstquark) && is_quark(*std::next(firstquark)) )
);
return std::make_pair(firstquark, std::next(firstquark));
}
//! returns index of most backward q-qbar jet
template<class Iterator>
int get_back_quark_jet(Event const & ev, Iterator firstquark){
// find jets at FO corresponding to the quarks
// technically this isn't necessary for LO
std::vector<int> const born_indices{ ev.particle_jet_indices() };
const auto firstquark_idx = std::distance(ev.begin_partons(), firstquark);
int const firstjet_idx = born_indices[firstquark_idx];
assert(firstjet_idx>0);
assert( born_indices[firstquark_idx+1] == firstjet_idx+1 );
return firstjet_idx;
}
//! returns index of most backward q-qbar jet
int getBackQuarkJet(Event const & ev){
const auto firstquark = get_central_qqbar(ev).first;
return get_back_quark_jet(ev, firstquark);
}
} // namespace
double PhaseSpacePoint::estimate_emission_rapidity_range(
Event const & ev
) const {
assert(std::is_sorted(begin(ev.jets()), end(ev.jets()), rapidity_less{}));
const double ymin = is_backward_g_to_h(ev)?
ev.outgoing().front().rapidity():
most_backward_FKL(ev.jets()).rapidity();
const double ymax = is_forward_g_to_h(ev)?
ev.outgoing().back().rapidity():
most_forward_FKL(ev.jets()).rapidity();
double delta_y = ymax - ymin;
// neglect tiny probability for emission between central qqbar pair
if(ev.type() == event_type::central_qqbar) {
const int qjet = getBackQuarkJet(ev);
delta_y -= ev.jets()[qjet+1].rapidity() - ev.jets()[qjet].rapidity();
}
assert(delta_y >= 0);
return delta_y;
}
double PhaseSpacePoint::estimate_ng_mean(Event const & ev) const {
// Formula derived from fit in arXiv:1805.04446 (see Fig. 2)
constexpr double GLUONS_PER_RAPIDITY = 0.975052;
return GLUONS_PER_RAPIDITY*estimate_emission_rapidity_range(ev);
}
int PhaseSpacePoint::sample_ng(Event const & event, RNG & ran){
const double ng_mean = estimate_ng_mean(event);
std::poisson_distribution<int> dist(ng_mean);
const int ng = dist(ran);
assert(ng >= 0);
assert(ng < NG_MAX);
weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng);
return ng;
}
void PhaseSpacePoint::boost_AWZH_boson_from(
std::vector<fastjet::PseudoJet> const & boosted_bosons, Event const & event
){
auto const & from = event.outgoing();
auto find_AWZH = [](Particle const & p){ return is_AWZH_boson(p); };
size_t boosted_idx = 0;
for(
auto AWZH_boson = std::find_if(begin(from), end(from), find_AWZH);
AWZH_boson != end(from);
AWZH_boson = std::find_if(++AWZH_boson, end(from), find_AWZH),
++boosted_idx
) {
auto insertion_point = std::lower_bound(
begin(outgoing_), end(outgoing_), *AWZH_boson, rapidity_less{}
);
- auto original_boson = *AWZH_boson;
- auto const & boosted_boson = boosted_bosons[boosted_idx];
- auto new_boson = original_boson;
- new_boson.p = boosted_boson;
-
-
+ auto original_boson = *AWZH_boson;
+ auto boosted_boson = boosted_bosons[boosted_idx];
+ auto new_boson = original_boson;
+ new_boson.p = boosted_boson;
+
+ // insert boson at new_idx
+ const int new_idx = std::distance(begin(outgoing_), insertion_point);
+ assert(new_idx >= 0); // insert invalidates distance
+ outgoing_.insert(insertion_point, new_boson);
+
+ // copy decay products
+ const int idx = std::distance(begin(from), AWZH_boson);
+ assert(idx >= 0);
+ const auto decay_it = event.decays().find(idx);
+ if(decay_it != end(event.decays())){
+ assert(outgoing_[new_idx].type == AWZH_boson->type);
+
+ auto decay_particles = decay_it -> second;
+ fastjet::PseudoJet sum;
+ for(auto & particle: decay_particles){
+ auto & p = particle.p;
+ // boost _to_ rest frame of input boson
+ p.unboost(original_boson.p);
+ // then boost _from_ rest frame of shuffled boson
+ p.boost(boosted_boson);
+ if(p.E() < std::abs(p.pz())){
+ throw std::underflow_error("Reshuffled decay with E<|pz|");
+ }
+ sum += p;
- // copy & boost decay products
- const int idx = std::distance(begin(from), AWZH_boson);
- assert(idx >= 0);
- const auto decay_it = event.decays().find(idx);
- if(decay_it == end(event.decays()))
- return;
- // insert boson at new_idx
- const int new_idx = std::distance(begin(outgoing_), insertion_point);
- assert(new_idx >= 0);
+ }
+ if(!nearby(boosted_boson, sum, boosted_boson.E())){
+ throw std::underflow_error("Boson and decays momenta do not match after reshuffling");
+ }
- assert(outgoing_[new_idx].type == AWZH_boson->type);
- auto decay_particles = decay_it -> second;
- fastjet::PseudoJet sum;
- for(auto & particle: decay_particles){
- auto & p = particle.p;
- // boost _to_ rest frame of input boson
- p.unboost(original_boson.p);
- // then boost _from_ rest frame of shuffled boson
- p.boost(boosted_boson);
- if(p.E() < std::abs(p.pz())){
- throw std::underflow_error("Reshuffled decay with E<|pz|");
- }
- sum += p;
- }
- if(!nearby(boosted_boson, sum, boosted_boson.E())){
- throw std::underflow_error("Boson and decays momenta do not match after reshuffling");
- }
- decays_.emplace(new_idx, decay_particles);
- }
- }
+ decays_.emplace(new_idx, decay_particles);
+ }
+ }
+ }
namespace {
template<class ConstIterator, class Iterator>
void label_extremal_qqbar(
ConstIterator born_begin, ConstIterator born_end,
Iterator first_out
){
// find born quarks
const auto firstquark = std::find_if(
born_begin, born_end-1,
[](Particle const & s){ return (is_anyquark(s)); }
);
assert(firstquark != born_end-1);
const auto secondquark = std::find_if(
firstquark+1, born_end,
[](Particle const & s){ return (is_anyquark(s)); }
);
assert(secondquark != born_end);
assert( ( is_quark(*firstquark) && is_antiquark(*secondquark) )
|| ( is_antiquark(*firstquark) && is_quark(*secondquark) ));
assert(first_out->type == ParticleID::gluon);
assert((first_out+1)->type == ParticleID::gluon);
// copy type from born
first_out->type = firstquark->type;
(first_out+1)->type = secondquark->type;
}
} // namespace
void PhaseSpacePoint::label_qqbar(Event const & event){
assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
assert(filter_partons(outgoing_).size() == outgoing_.size());
if(qqbarb_){
label_extremal_qqbar(event.outgoing().cbegin(), event.outgoing().cend(),
outgoing_.begin() );
return;
}
if(qqbarf_){ // same as qqbarb with reversed order
label_extremal_qqbar( event.outgoing().crbegin(), event.outgoing().crend(),
outgoing_.rbegin() );
return;
}
// central qqbar
const auto firstquark = get_central_qqbar(event).first;
// find jets at FO corresponding to the quarks
// technically this isn't necessary for LO
const auto firstjet_idx = get_back_quark_jet(event, firstquark);
// find corresponding jets after resummation
fastjet::ClusterSequence cs{to_PseudoJet(outgoing_), param_.jet_param.def};
auto const jets = fastjet::sorted_by_rapidity(
cs.inclusive_jets( param_.jet_param.min_pt ));
std::vector<int> const resum_indices{ cs.particle_jet_indices({jets}) };
// assert that jets didn't move
assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx )->rapidity(),
jets[ firstjet_idx ].rapidity(), 1e-2) );
assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx+1 )->rapidity(),
jets[ firstjet_idx+1 ].rapidity(), 1e-2) );
// find last partons in first (central) jet
size_t idx_out = 0;
for(size_t i=resum_indices.size()-2; i>0; --i)
if(resum_indices[i] == firstjet_idx){
idx_out = i;
break;
}
assert(idx_out != 0);
// check that there is sufficient pt in jets from the quarks
const double minpartonjetpt = 1. - param_.soft_pt_regulator;
if (outgoing_[idx_out].p.pt()<minpartonjetpt*( event.jets().cbegin()+firstjet_idx )->pt()){
weight_=0.;
status_ = StatusCode::wrong_jets;
return;
}
if (outgoing_[idx_out+1].p.pt()<minpartonjetpt*( event.jets().cbegin()+firstjet_idx+1 )->pt()){
weight_=0.;
status_ = StatusCode::wrong_jets;
return;
}
// check that no additional emission between jets
// such configurations are possible if we have an gluon gets generated
// inside the rapidities of the qqbar chain, but clusted to a
// differnet/outside jet. Changing this is non trivial
if(resum_indices[idx_out+1] != resum_indices[idx_out]+1){
weight_=0.;
status_ = StatusCode::gluon_in_qqbar;
return;
}
outgoing_[idx_out].type = firstquark->type;
outgoing_[idx_out+1].type = std::next(firstquark)->type;
}
void PhaseSpacePoint::label_quarks(Event const & ev){
const auto WZEmit = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & s){ return (std::abs(s.type) == pid::Wp || s.type == pid::Z_photon_mix); }
);
if (WZEmit != end(ev.outgoing())){
if(!qqbarb_) {
const size_t backward_FKL_idx = unob_?1:0;
const auto backward_FKL = std::next(ev.begin_partons(), backward_FKL_idx);
outgoing_[backward_FKL_idx].type = backward_FKL->type;
}
if(!qqbarf_) {
const size_t forward_FKL_idx = unof_?1:0;
const auto forward_FKL = std::prev(ev.end_partons(), 1+forward_FKL_idx);
outgoing_.rbegin()[unof_].type = forward_FKL->type; // NOLINT
}
} else {
if(!is_backward_g_to_h(ev)) {
most_backward_FKL(outgoing_).type = ev.incoming().front().type;
}
if(!is_forward_g_to_h(ev)) {
most_forward_FKL(outgoing_).type = ev.incoming().back().type;
}
}
if(qqbar_mid_||qqbarb_||qqbarf_){
label_qqbar(ev);
}
}
PhaseSpacePoint::PhaseSpacePoint(
Event const & ev, PhaseSpacePointConfig conf, RNG & ran
):
unob_{ev.type() == event_type::unob},
unof_{ev.type() == event_type::unof},
qqbarb_{ev.type() == event_type::qqbar_exb},
qqbarf_{ev.type() == event_type::qqbar_exf},
qqbar_mid_{ev.type() == event_type::qqbar_mid},
param_{std::move(conf)},
status_{unspecified}
{
// legacy code: override new variable with old
if(param_.max_ext_soft_pt_fraction){
param_.soft_pt_regulator = *param_.max_ext_soft_pt_fraction;
param_.max_ext_soft_pt_fraction = {};
}
weight_ = 1;
auto const & Born_jets = ev.jets();
const int ng = sample_ng(ev, ran);
weight_ /= std::tgamma(ng + 1);
const int ng_jets = sample_ng_jets(ev, ng, ran);
std::vector<fastjet::PseudoJet> out_partons = gen_non_jet(
ng - ng_jets, CMINPT, param_.jet_param.min_pt, ran
);
const auto qperp = std::accumulate(
begin(out_partons), end(out_partons),
fastjet::PseudoJet{}
);
std::vector<fastjet::PseudoJet> jets;
std::vector<fastjet::PseudoJet> bosons;
std::tie(jets, bosons) = reshuffle(ev, qperp);
if(weight_ == 0.) {
status_ = failed_reshuffle;
return;
}
if(! pass_resummation_cuts(jets)){
status_ = failed_resummation_cuts;
weight_ = 0.;
return;
}
// split jets in multiple partons
std::vector<fastjet::PseudoJet> jet_partons = split(
ev, jets, ng_jets, ran
);
if(weight_ == 0.) {
status_ = StatusCode::failed_split;
return;
}
const double ymin = is_backward_g_to_h(ev)?
ev.outgoing().front().rapidity():
most_backward_FKL(jet_partons).rapidity()
;
const double ymax = is_forward_g_to_h(ev)?
ev.outgoing().back().rapidity():
most_forward_FKL(jet_partons).rapidity()
;
if(qqbar_mid_){
const int qqbar_backjet = getBackQuarkJet(ev);
rescale_qqbar_rapidities(
out_partons, jets,
ymin, ymax,
qqbar_backjet
);
}
else{
rescale_rapidities(out_partons, ymin, ymax);
}
if(! cluster_jets(out_partons).empty()){
weight_ = 0.;
status_ = StatusCode::empty_jets;
return;
}
std::sort(begin(out_partons), end(out_partons), rapidity_less{});
assert(
std::is_sorted(begin(jet_partons), end(jet_partons), rapidity_less{})
);
const auto first_jet_parton = out_partons.insert(
end(out_partons), begin(jet_partons), end(jet_partons)
);
std::inplace_merge(
begin(out_partons), first_jet_parton, end(out_partons), rapidity_less{}
);
if(! jets_ok(ev, out_partons)){
weight_ = 0.;
status_ = StatusCode::wrong_jets;
return;
}
weight_ *= phase_space_normalisation(Born_jets.size(), out_partons.size());
outgoing_.reserve(out_partons.size() + 1); // one slot for possible A, W, Z, H
for( auto it = std::make_move_iterator(out_partons.begin());
it != std::make_move_iterator(out_partons.end());
++it
){
outgoing_.emplace_back( Particle{pid::gluon, *it, {}});
}
assert(!outgoing_.empty());
label_quarks(ev);
if(weight_ == 0.) {
//! @TODO optimise s.t. this is not possible
// status is handled internally
return;
}
// reattach boson & decays
if(!bosons.empty()){
try {
boost_AWZH_boson_from(bosons, ev);
} catch (std::underflow_error const & e){
weight_ = 0.;
status_ = StatusCode::failed_reshuffle;
return;
}
}
reconstruct_incoming(ev.incoming());
status_ = StatusCode::good;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_non_jet(
int const ng_non_jet, double const ptmin, double const ptmax, RNG & ran
){
// heuristic parameters for pt sampling
const double ptpar = 1.3 + ng_non_jet/5.;
const double temp1 = std::atan((ptmax - ptmin)/ptpar);
std::vector<fastjet::PseudoJet> partons(ng_non_jet);
for(int i = 0; i < ng_non_jet; ++i){
const double r1 = ran.flat();
const double pt = ptmin + ptpar*std::tan(r1*temp1);
const double temp2 = std::cos(r1*temp1);
const double phi = 2*M_PI*ran.flat();
weight_ *= 2.0*M_PI*pt*ptpar*temp1/(temp2*temp2);
// we don't know the allowed rapidity span yet,
// set a random value to be rescaled later on
const double y = ran.flat();
partons[i].reset_PtYPhiM(pt, y, phi);
// Set user index higher than any jet-parton index
// in order to assert that these are not inside jets
partons[i].set_user_index(i + 1 + NG_MAX);
assert(ptmin-1e-5 <= partons[i].pt() && partons[i].pt() <= ptmax+1e-5);
}
assert(std::all_of(partons.cbegin(), partons.cend(), is_nonjet_parton));
return sorted_by_rapidity(partons);
}
void PhaseSpacePoint::rescale_qqbar_rapidities(
std::vector<fastjet::PseudoJet> & out_partons,
std::vector<fastjet::PseudoJet> const & jets,
const double ymin1, const double ymax2,
const int qqbar_backjet
){
const double ymax1 = jets[qqbar_backjet].rapidity();
const double ymin2 = jets[qqbar_backjet+1].rapidity();
constexpr double ep = 1e-7;
const double tot_y = ymax1 - ymin1 + ymax2 - ymin2;
std::vector<std::reference_wrapper<fastjet::PseudoJet>> refpart(
out_partons.begin(), out_partons.end());
double ratio = (ymax1 - ymin1)/tot_y;
const auto gap{ std::find_if(refpart.begin(), refpart.end(),
[ratio](fastjet::PseudoJet const & p){
return (p.rapidity()>=ratio);} ) };
double ymin = ymin1;
double ymax = ymax1;
double dy = ymax - ymin - 2*ep;
double offset = 0.;
for(auto it_part=refpart.begin(); it_part<refpart.end(); ++it_part){
if(it_part == gap){
ymin = ymin2;
ymax = ymax2;
dy = ymax - ymin - 2*ep;
offset = ratio;
ratio = 1-ratio;
}
fastjet::PseudoJet & part = *it_part;
assert(offset <= part.rapidity() && part.rapidity() < ratio+offset);
const double y = ymin + ep + dy*((part.rapidity()-offset)/ratio);
part.reset_momentum_PtYPhiM(part.pt(), y, part.phi());
weight_ *= tot_y-4.*ep;
assert(ymin <= part.rapidity() && part.rapidity() <= ymax);
}
assert(is_sorted(begin(out_partons), end(out_partons), rapidity_less{}));
}
void PhaseSpacePoint::rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
){
constexpr double ep = 1e-7;
for(auto & parton: partons){
assert(0 <= parton.rapidity() && parton.rapidity() <= 1);
const double dy = ymax - ymin - 2*ep;
const double y = ymin + ep + dy*parton.rapidity();
parton.reset_momentum_PtYPhiM(parton.pt(), y, parton.phi());
weight_ *= dy;
assert(ymin <= parton.rapidity() && parton.rapidity() <= ymax);
}
}
namespace {
template<typename T, typename... Rest>
auto min(T const & a, T const & b, Rest&&... r) {
using std::min;
return min(a, min(b, std::forward<Rest>(r)...));
}
}
double PhaseSpacePoint::probability_in_jet(Event const & ev) const{
const double dy = estimate_emission_rapidity_range(ev);
const double R = param_.jet_param.def.R();
// jets into which we predominantly emit
const int njets = ev.jets().size() - unof_ - unob_ - qqbarb_ - qqbarf_; //NOLINT
assert(njets >= 1);
const size_t nextremal_jets = std::min(njets, 2);
const double p_J_y_large = (njets - nextremal_jets/2.)*R*R/(2.*dy);
const double p_J_y0 = njets*R/M_PI;
return min(p_J_y_large, p_J_y0, 1.);
}
int PhaseSpacePoint::sample_ng_jets(Event const & event, int ng, RNG & ran){
const double p_J = probability_in_jet(event);
std::binomial_distribution<> bin_dist(ng, p_J);
const int ng_J = bin_dist(ran);
weight_ *= std::pow(p_J, -ng_J)*std::pow(1 - p_J, ng_J - ng);
return ng_J;
}
std::pair< std::vector<fastjet::PseudoJet>,
std::vector<fastjet::PseudoJet> >
PhaseSpacePoint::reshuffle(
Event const & ev,
fastjet::PseudoJet const & q
){
// Create a copy of the outgoing momenta not containing decay products
std::vector<fastjet::PseudoJet const *> born_momenta;
born_momenta.reserve(ev.jets().size());
std::transform(ev.jets().cbegin(), ev.jets().cend(),
back_inserter(born_momenta),
[](fastjet::PseudoJet const & t) { return &t; });
auto bosons = filter_AWZH_bosons(ev.outgoing());
std::vector<fastjet::PseudoJet const *> p_boson_momenta;
std::transform(bosons.cbegin(), bosons.cend(),
back_inserter(p_boson_momenta),
[](Particle const & t) { return &(t.p); });
std::vector<fastjet::PseudoJet> boson_momenta;
std::transform(bosons.cbegin(), bosons.cend(),
back_inserter(boson_momenta),
[](Particle const & t) { return t.p; });
// reshuffle all momenta
if(q == fastjet::PseudoJet{0, 0, 0, 0}) return {ev.jets(), boson_momenta};
// add boson to reshuffling
if(!bosons.empty()) {
born_momenta.insert( born_momenta.end(), p_boson_momenta.begin(), p_boson_momenta.end() );
}
auto shuffle_momenta = resummation_jet_momenta(born_momenta, q);
if(shuffle_momenta.empty()){
weight_ = 0;
return {};
}
// additional Jacobian to ensure Born integration over delta gives 1
weight_ *= resummation_jet_weight(born_momenta, q);
// take out boson again
if(!bosons.empty()) {
std::vector<fastjet::PseudoJet> shuffle_bosons;
for(size_t i = 0; i < bosons.size(); ++i) {
shuffle_bosons.push_back(shuffle_momenta.back());
shuffle_momenta.pop_back();
}
std::reverse(shuffle_bosons.begin(), shuffle_bosons.end());
return {shuffle_momenta, shuffle_bosons};
}
return {shuffle_momenta, {}};
}
std::vector<int> PhaseSpacePoint::distribute_jet_partons(
int ng_jets, std::vector<fastjet::PseudoJet> const & jets, RNG & ran
){
size_t first_valid_jet = 0;
size_t num_valid_jets = jets.size();
const double R_eff = 5./3.*param_.jet_param.def.R();
// if there is an unordered jet too far away from the FKL jets
// then extra gluon constituents of the unordered jet would
// violate the FKL rapidity ordering
if((unob_||qqbarb_) && jets[0].delta_R(jets[1]) > R_eff){
++first_valid_jet;
--num_valid_jets;
}
else if((unof_||qqbarf_) && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){
--num_valid_jets;
}
std::vector<int> np(jets.size(), 1);
for(int i = 0; i < ng_jets; ++i){
++np[first_valid_jet + ran.flat() * num_valid_jets];
}
weight_ *= std::pow(num_valid_jets, ng_jets);
return np;
}
#ifndef NDEBUG
namespace {
bool tagged_FKL_backward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return std::find_if(
begin(jet_partons), end(jet_partons),
[](fastjet::PseudoJet const & p){
return p.user_index() == UID::backward_fkl;
}
) != end(jet_partons);
}
bool tagged_FKL_forward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
// the most forward FKL parton is most likely near the end of jet_partons;
// start search from there
return std::find_if(
jet_partons.rbegin(), jet_partons.rend(),
[](fastjet::PseudoJet const & p){
return p.user_index() == UID::forward_fkl;
}
) != jet_partons.rend();
}
} // namespace
#endif
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
Event const & Born_event,
std::vector<fastjet::PseudoJet> const & jets,
int ng_jets
, RNG & ran
){
return split(
Born_event,
jets,
distribute_jet_partons(ng_jets, jets, ran),
ran
);
}
bool PhaseSpacePoint::pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const{
if(ext_parton.pt() < param_.min_extparton_pt) return false;
return (ext_parton - jet).pt()/jet.pt() < param_.soft_pt_regulator;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
Event const & Born_event,
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np,
RNG & ran
){
assert(! jets.empty());
assert(jets.size() == np.size());
assert(pass_resummation_cuts(jets));
constexpr auto no_such_jet_idx = std::numeric_limits<std::size_t>::max();
const size_t most_backward_FKL_idx = is_backward_g_to_h(Born_event)?
no_such_jet_idx: // we have backward Higgs instead of FKL jet
(0 + unob_ + qqbarb_); // NOLINT
const size_t most_forward_FKL_idx = is_forward_g_to_h(Born_event)?
no_such_jet_idx: // we have forward Higgs instead of FKL jet
(jets.size() - 1 - unof_ - qqbarf_); // NOLINT
const size_t qqbar_jet_idx = qqbar_mid_?
getBackQuarkJet(Born_event):
no_such_jet_idx;
auto const & jet = param_.jet_param;
const JetSplitter jet_splitter{jet.def, jet.min_pt};
std::vector<fastjet::PseudoJet> jet_partons;
// randomly distribute jet gluons among jets
for(size_t i = 0; i < jets.size(); ++i){
auto split_res = jet_splitter.split(jets[i], np[i], ran);
weight_ *= split_res.weight;
if(weight_ == 0) return {};
assert(
std::all_of(
begin(split_res.constituents), end(split_res.constituents),
is_jet_parton
)
);
const auto first_new_parton = jet_partons.insert(
end(jet_partons),
begin(split_res.constituents), end(split_res.constituents)
);
// mark uno and extremal FKL emissions here so we can check
// their position once all emissions are generated
// also mark qqbar_mid partons, and apply appropriate pt cut.
auto extremal = end(jet_partons);
if (i == most_backward_FKL_idx){ //FKL backward emission
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(UID::backward_fkl);
}
else if(((unob_ || qqbarb_) && i == 0)){
// unordered/qqbarb
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index((unob_)?UID::unob:UID::qqbarb);
}
else if (i == most_forward_FKL_idx){
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(UID::forward_fkl);
}
else if(((unof_ || qqbarf_) && i == jets.size() - 1)){
// unordered/qqbarf
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index((unof_)?UID::unof:UID::qqbarf);
}
else if((qqbar_mid_ && i == qqbar_jet_idx)){
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(UID::qqbar_mid1);
}
else if((qqbar_mid_ && i == qqbar_jet_idx+1)){
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(UID::qqbar_mid2);
}
if(
extremal != end(jet_partons)
&& !pass_extremal_cuts(*extremal, jets[i])
){
weight_ = 0;
return {};
}
}
assert(is_backward_g_to_h(Born_event) || tagged_FKL_backward(jet_partons));
assert(is_forward_g_to_h(Born_event) || tagged_FKL_forward(jet_partons));
std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
if(
!extremal_ok(Born_event, jet_partons)
|| !split_preserved_jets(jets, jet_partons)
){
weight_ = 0.;
return {};
}
return jet_partons;
}
bool PhaseSpacePoint::extremal_ok(
Event const & Born_event,
std::vector<fastjet::PseudoJet> const & partons
) const{
assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
if(unob_ && partons.front().user_index() != UID::unob) return false;
if(unof_ && partons.back().user_index() != UID::unof) return false;
if(qqbarb_ && partons.front().user_index() != UID::qqbarb) return false;
if(qqbarf_ && partons.back().user_index() != UID::qqbarf) return false;
if(is_backward_g_to_h(Born_event)) {
if(partons.front().rapidity() < Born_event.outgoing().front().rapidity()){
return false;
}
} else if(most_backward_FKL(partons).user_index() != UID::backward_fkl) {
return false;
}
if(is_forward_g_to_h(Born_event)) {
return partons.back().rapidity() <= Born_event.outgoing().back().rapidity();
}
return most_forward_FKL(partons).user_index() == UID::forward_fkl;
}
bool PhaseSpacePoint::split_preserved_jets(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<fastjet::PseudoJet> const & jet_partons
) const{
assert(std::is_sorted(begin(jets), end(jets), rapidity_less{}));
const auto split_jets = cluster_jets(jet_partons);
// this can happen if two overlapping jets
// are both split into more than one parton
if(split_jets.size() != jets.size()) return false;
for(size_t i = 0; i < split_jets.size(); ++i){
// this can happen if there are two overlapping jets
// and a parton is assigned to the "wrong" jet
if(!nearby_ep(jets[i].rapidity(), split_jets[i].rapidity(), 1e-2)){
return false;
}
}
return true;
}
template<class Particle>
Particle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> const & partons
) const{
return partons[0 + unob_ + qqbarb_];
}
template<class Particle>
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
const size_t idx = partons.size() - 1 - unof_ - qqbarf_;
assert(idx < partons.size());
return partons[idx];
}
template<class Particle>
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
return partons[0 + unob_ + qqbarb_];
}
template<class Particle>
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
const size_t idx = partons.size() - 1 - unof_ - qqbarf_;
assert(idx < partons.size());
return partons[idx];
}
bool PhaseSpacePoint::contains_idx(
fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
) const {
auto const & constituents = jet.constituents();
const int idx = parton.user_index();
const bool injet = std::find_if(
begin(constituents), end(constituents),
[idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
) != end(constituents);
const double minpartonjetpt = 1. - param_.soft_pt_regulator;
return ((parton.pt()>minpartonjetpt*jet.pt())&&injet);
}
bool PhaseSpacePoint::jets_ok(
Event const & Born_event,
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
if(jets.size() != Born_event.jets().size()) return false;
int in_jet = 0;
for(auto const & jet : jets){
assert(jet.has_constituents());
for(auto && parton: jet.constituents()){
if(is_nonjet_parton(parton)) return false;
}
in_jet += jet.constituents().size();
}
const int expect_in_jet = std::count_if(
partons.cbegin(), partons.cend(), is_jet_parton
);
if(in_jet != expect_in_jet) return false;
// note that PseudoJet::contains does not work here
if(
!is_backward_g_to_h(Born_event) &&
!contains_idx(most_backward_FKL(jets), most_backward_FKL(partons))
) return false;
if(
!is_forward_g_to_h(Born_event)
&& !contains_idx(most_forward_FKL(jets), most_forward_FKL(partons))
) return false;
if(unob_ && !contains_idx(jets.front(), partons.front())) return false;
if(qqbarb_ && !contains_idx(jets.front(), partons.front())) return false;
if(unof_ && !contains_idx(jets.back(), partons.back())) return false;
if(qqbarf_ && !contains_idx(jets.back(), partons.back())) return false;
#ifndef NDEBUG
for(size_t i = 0; i < jets.size(); ++i){
assert(nearby_ep(jets[i].rapidity(), Born_event.jets()[i].rapidity(), 1e-2));
}
#endif
return true;
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<Particle, 2> const & Born_incoming
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
for(size_t i = 0; i < incoming_.size(); ++i){
incoming_[i].type = Born_incoming[i].type;
}
assert(momentum_conserved());
}
bool PhaseSpacePoint::momentum_conserved() const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
const double norm = diff.E();
for(auto const & out: outgoing()) diff -= out.p;
return nearby(diff, fastjet::PseudoJet{}, norm);
}
} //namespace HEJ
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 9:14 PM (1 d, 30 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806218
Default Alt Text
(35 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment