diff --git a/include/RHEJ/resummation_jet_momenta.hh b/include/RHEJ/resummation_jet_momenta.hh
index 03213f1..12cda67 100644
--- a/include/RHEJ/resummation_jet_momenta.hh
+++ b/include/RHEJ/resummation_jet_momenta.hh
@@ -1,21 +1,27 @@
/** \file
* \brief Function to calculate the momenta of resummation jets
*/
#pragma once
#include "RHEJ/utility.hh"
namespace RHEJ{
+ enum class ReshufflingAlgorithm {
+ linear_in_born,
+ linear_in_resummed
+ };
+
/**
* \brief Calculate the resummation jet momenta
- * @param p_born Born Jet Momenta
+ * @param p_born born Jet Momenta
* @param qperp Sum of non-jet Parton Transverse Momenta
* @returns Resummation Jet Momenta
*/
std::vector<fastjet::PseudoJet> resummation_jet_momenta(
std::vector<fastjet::PseudoJet> const & p_born,
- fastjet::PseudoJet const & qperp
+ fastjet::PseudoJet const & qperp,
+ ReshufflingAlgorithm algo = ReshufflingAlgorithm::linear_in_born
);
}
diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc
index 3ac02ab..f60b75f 100644
--- a/src/PhaseSpacePoint.cc
+++ b/src/PhaseSpacePoint.cc
@@ -1,537 +1,559 @@
#include "RHEJ/PhaseSpacePoint.hh"
#include <random>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
#include "RHEJ/Constants.hh"
#include "RHEJ/resummation_jet_momenta.hh"
#include "RHEJ/Jacobian.hh"
#include "RHEJ/uno.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/kinematics.hh"
namespace RHEJ{
namespace {
constexpr int max_jet_user_idx = PhaseSpacePoint::ng_max;
bool is_nonjet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() > max_jet_user_idx;
}
bool is_jet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() <= max_jet_user_idx;
}
// user indices for partons with extremal rapidity
constexpr int unob_idx = -5;
constexpr int unof_idx = -4;
constexpr int backward_FKL_idx = -3;
constexpr int forward_FKL_idx = -2;
}
namespace {
double estimate_ng_mean(std::vector<fastjet::PseudoJet> const & Born_jets){
const double delta_y =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
assert(delta_y > 0);
// Formula derived from fit in reversed HEJ intro paper
return 0.975052*delta_y;
}
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
return cs.inclusive_jets(param_.jet_param.min_pt);
}
bool PhaseSpacePoint::pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const{
return cluster_jets(jets).size() == jets.size();
}
int PhaseSpacePoint::sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets){
const double ng_mean = estimate_ng_mean(Born_jets);
std::poisson_distribution<int> dist(ng_mean);
const int ng = dist(ran_.get());
assert(ng >= 0);
assert(ng < ng_max);
weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng);
return ng;
}
void PhaseSpacePoint::copy_AWZH_boson_from(Event const & event){
auto const & from = event.outgoing();
const auto AWZH_boson = std::find_if(
begin(from), end(from),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson == end(from)) return;
auto insertion_point = std::lower_bound(
begin(outgoing_), end(outgoing_), *AWZH_boson, rapidity_less{}
);
outgoing_.insert(insertion_point, *AWZH_boson);
// copy decay products
const int idx = std::distance(begin(from), AWZH_boson);
const auto decay_it = event.decays().find(idx);
if(decay_it != end(event.decays())){
const int new_idx = std::distance(begin(outgoing_), insertion_point);
assert(outgoing_[new_idx].type == AWZH_boson->type);
decays_.emplace(new_idx, decay_it->second);
}
assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
}
PhaseSpacePoint::PhaseSpacePoint(
Event const & ev, PhaseSpacePointConfig conf, RHEJ::RNG & ran
):
unob_{has_unob_gluon(ev.incoming(), ev.outgoing())},
unof_{!unob_ && has_unof_gluon(ev.incoming(), ev.outgoing())},
param_{std::move(conf)},
ran_{ran}
{
weight_ = 1;
const auto Born_jets = sorted_by_rapidity(ev.jets());
const int ng = sample_ng(Born_jets);
weight_ /= std::tgamma(ng + 1);
const int ng_jets = sample_ng_jets(ng, Born_jets);
std::vector<fastjet::PseudoJet> out_partons = gen_non_jet(
ng - ng_jets, CMINPT, param_.jet_param.min_pt
);
{
const auto qperp = std::accumulate(
begin(out_partons), end(out_partons),
fastjet::PseudoJet{}
);
const auto jets = reshuffle(Born_jets, qperp);
if(weight_ == 0.) return;
if(! pass_resummation_cuts(jets)){
weight_ = 0.;
return;
}
std::vector<fastjet::PseudoJet> jet_partons = split(jets, ng_jets);
if(weight_ == 0.) return;
rescale_rapidities(
out_partons,
most_backward_FKL(jet_partons).rapidity(),
most_forward_FKL(jet_partons).rapidity()
);
if(! cluster_jets(out_partons).empty()){
weight_ = 0.;
return;
}
std::sort(begin(out_partons), end(out_partons), rapidity_less{});
assert(
std::is_sorted(begin(jet_partons), end(jet_partons), rapidity_less{})
);
const auto first_jet_parton = out_partons.insert(
end(out_partons), begin(jet_partons), end(jet_partons)
);
std::inplace_merge(
begin(out_partons), first_jet_parton, end(out_partons), rapidity_less{}
);
}
if(! jets_ok(Born_jets, out_partons)){
weight_ = 0.;
return;
}
weight_ *= phase_space_normalisation(Born_jets.size(), out_partons.size());
outgoing_.reserve(out_partons.size() + 1); // one slot for possible A, W, Z, H
for(auto & p: out_partons){
outgoing_.emplace_back(Particle{pid::gluon, std::move(p)});
}
most_backward_FKL(outgoing_).type = ev.incoming().front().type;
most_forward_FKL(outgoing_).type = ev.incoming().back().type;
copy_AWZH_boson_from(ev);
assert(!outgoing_.empty());
reconstruct_incoming(ev.incoming());
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_non_jet(
int count, double ptmin, double ptmax
){
// heuristic parameters for pt sampling
const double ptpar = 1.3 + count/5.;
const double temp1 = atan((ptmax - ptmin)/ptpar);
std::vector<fastjet::PseudoJet> partons(count);
for(size_t i = 0; i < (size_t) count; ++i){
const double r1 = ran_.get().flat();
const double pt = ptmin + ptpar*tan(r1*temp1);
const double temp2 = cos(r1*temp1);
const double phi = 2*M_PI*ran_.get().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_.get().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 partons;
}
void PhaseSpacePoint::rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
){
constexpr double ep = 1e-7;
for(auto & parton: partons){
assert(0 <= parton.rapidity() && parton.rapidity() <= 1);
const double dy = ymax - ymin - 2*ep;
const double y = ymin + ep + dy*parton.rapidity();
parton.reset_momentum_PtYPhiM(parton.pt(), y, parton.phi());
weight_ *= dy;
assert(ymin <= parton.rapidity() && parton.rapidity() <= ymax);
}
}
namespace {
template<typename T, typename... Rest>
auto min(T const & a, T const & b, Rest&&... r) {
using std::min;
return min(a, min(b, std::forward<Rest>(r)...));
}
}
double PhaseSpacePoint::probability_in_jet(
std::vector<fastjet::PseudoJet> const & Born_jets
) const{
assert(std::is_sorted(begin(Born_jets), end(Born_jets), rapidity_less{}));
assert(Born_jets.size() >= 2);
const double dy =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
const double R = param_.jet_param.def.R();
const int njets = Born_jets.size();
const double p_J_y_large = (njets-1)*R*R/(2.*dy);
const double p_J_y0 = njets*R/M_PI;
return min(p_J_y_large, p_J_y0, 1.);
}
int PhaseSpacePoint::sample_ng_jets(
int ng, std::vector<fastjet::PseudoJet> const & Born_jets
){
const double p_J = probability_in_jet(Born_jets);
std::binomial_distribution<> bin_dist(ng, p_J);
const int ng_J = bin_dist(ran_.get());
weight_ *= std::pow(p_J, -ng_J)*std::pow(1 - p_J, ng_J - ng);
return ng_J;
}
+ namespace {
+ double reshuffling_weight(
+ std::vector<fastjet::PseudoJet> const & Born_jets,
+ std::vector<fastjet::PseudoJet> const & jets,
+ fastjet::PseudoJet const & q,
+ ReshufflingAlgorithm algo
+ ) {
+ switch(algo) {
+ case ReshufflingAlgorithm::linear_in_born:
+ // transform delta functions to integration over resummation momenta
+ return 1/Jacobian(jets, q);
+ case ReshufflingAlgorithm::linear_in_resummed:
+ // additional Jacobian to ensure Born integration over delta gives 1
+ return Jacobian(Born_jets, -1.*q);
+ default:;
+ }
+ throw std::logic_error{"unreachable"};
+ }
+ }
+
std::vector<fastjet::PseudoJet>
PhaseSpacePoint::reshuffle(
std::vector<fastjet::PseudoJet> const & Born_jets,
fastjet::PseudoJet const & q
){
+ static constexpr auto algo = ReshufflingAlgorithm::linear_in_born;
+
if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets;
- std::vector<fastjet::PseudoJet> jets = resummation_jet_momenta(Born_jets, q);
+ const auto jets = resummation_jet_momenta(Born_jets, q, algo);
if(jets.empty()){
weight_ = 0;
return {};
}
- // transform delta functions to integration over resummation momenta
- weight_ /= Jacobian(jets, q);
+
+ weight_ *= reshuffling_weight(Born_jets, jets, q, algo);
return jets;
}
std::vector<int> PhaseSpacePoint::distribute_jet_partons(
int ng_jets, std::vector<fastjet::PseudoJet> const & jets
){
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_ && jets[0].delta_R(jets[1]) > R_eff){
++first_valid_jet;
--num_valid_jets;
}
else if(unof_ && 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_.get().flat() * num_valid_jets];
}
weight_ *= std::pow(num_valid_jets, ng_jets);
return np;
}
#ifndef NDEBUG
namespace{
bool tagged_FKL_backward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return std::find_if(
begin(jet_partons), end(jet_partons),
[](fastjet::PseudoJet const & p){
return p.user_index() == backward_FKL_idx;
}
) != end(jet_partons);
}
bool tagged_FKL_forward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
// the most forward FKL parton is most likely near the end of jet_partons;
// start search from there
return std::find_if(
jet_partons.rbegin(), jet_partons.rend(),
[](fastjet::PseudoJet const & p){
return p.user_index() == forward_FKL_idx;
}
) != jet_partons.rend();
}
bool tagged_FKL_extremal(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return tagged_FKL_backward(jet_partons) && tagged_FKL_forward(jet_partons);
}
} // namespace anonymous
#endif
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
int ng_jets
){
return split(jets, distribute_jet_partons(ng_jets, jets));
}
bool PhaseSpacePoint::pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const{
if(ext_parton.pt() < param_.min_extparton_pt) return false;
return (ext_parton - jet).pt()/jet.pt() < param_.max_ext_soft_pt_fraction;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np
){
assert(! jets.empty());
assert(jets.size() == np.size());
assert(pass_resummation_cuts(jets));
const size_t most_backward_FKL_idx = 0 + unob_;
const size_t most_forward_FKL_idx = jets.size() - 1 - unof_;
const auto & jet = param_.jet_param;
const JetSplitter jet_splitter{jet.def, jet.min_pt, ran_};
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]);
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
auto extremal = end(jet_partons);
if((unob_ && i == 0) || i == most_backward_FKL_idx){
// unordered or FKL backward emission
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(
(i == most_backward_FKL_idx)?backward_FKL_idx:unob_idx
);
}
else if((unof_ && i == jets.size() - 1) || i == most_forward_FKL_idx){
// unordered or FKL forward emission
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
extremal->set_user_index(
(i == most_forward_FKL_idx)?forward_FKL_idx:unof_idx
);
}
if(
extremal != end(jet_partons)
&& !pass_extremal_cuts(*extremal, jets[i])
){
weight_ = 0;
return {};
}
}
assert(tagged_FKL_extremal(jet_partons));
std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
if(
!extremal_ok(jet_partons)
|| !split_preserved_jets(jets, jet_partons)
){
weight_ = 0.;
return {};
}
return jet_partons;
}
bool PhaseSpacePoint::extremal_ok(
std::vector<fastjet::PseudoJet> const & partons
) const{
assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
if(unob_ && partons.front().user_index() != unob_idx) return false;
if(unof_ && partons.back().user_index() != unof_idx) return false;
return
most_backward_FKL(partons).user_index() == backward_FKL_idx
&& most_forward_FKL(partons).user_index() == forward_FKL_idx;
}
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 = sorted_by_rapidity(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_];
}
template<class Particle>
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
const size_t idx = partons.size() - 1 - unof_;
assert(idx < partons.size());
return partons[idx];
}
template<class Particle>
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
return partons[0 + unob_];
}
template<class Particle>
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
const size_t idx = partons.size() - 1 - unof_;
assert(idx < partons.size());
return partons[idx];
}
namespace {
bool contains_idx(
fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
){
auto const & constituents = jet.constituents();
const int idx = parton.user_index();
return std::find_if(
begin(constituents), end(constituents),
[idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
) != end(constituents);
}
}
/**
* final jet test:
* - number of jets must match Born kinematics
* - no partons designated as nonjet may end up inside jets
* - all other outgoing partons *must* end up inside jets
* - the extremal (in rapidity) partons must be inside the extremal jets
* - rapidities must be the same (by construction)
*/
bool PhaseSpacePoint::jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
if(jets.size() != Born_jets.size()) return false;
int in_jet = 0;
for(size_t i = 0; i < jets.size(); ++i){
assert(jets[i].has_constituents());
for(auto && parton: jets[i].constituents()){
if(is_nonjet_parton(parton)) return false;
}
in_jet += jets[i].constituents().size();
}
const int expect_in_jet = std::count_if(
partons.cbegin(), partons.cend(), is_jet_parton
);
if(in_jet != expect_in_jet) return false;
// note that PseudoJet::contains does not work here
if(! (
contains_idx(most_backward_FKL(jets), most_backward_FKL(partons))
&& contains_idx(most_forward_FKL(jets), most_forward_FKL(partons))
)) return false;
if(unob_ && !contains_idx(jets.front(), partons.front())) return false;
if(unof_ && !contains_idx(jets.back(), partons.back())) return false;
for(size_t i = 0; i < jets.size(); ++i){
assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2));
}
return true;
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<Particle, 2> const & Born_incoming
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
for(size_t i = 0; i < incoming_.size(); ++i){
incoming_[i].type = Born_incoming[i].type;
}
assert(momentum_conserved());
}
double PhaseSpacePoint::phase_space_normalisation(
int num_Born_jets, int num_out_partons
) const{
return pow(16*pow(M_PI,3), num_Born_jets - num_out_partons);
}
bool PhaseSpacePoint::momentum_conserved() const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
const double norm = diff.E();
for(auto const & out: outgoing()) diff -= out.p;
return nearby(diff, fastjet::PseudoJet{}, norm);
}
} //namespace RHEJ
diff --git a/src/resummation_jet_momenta.cc b/src/resummation_jet_momenta.cc
index fc02847..61206a0 100644
--- a/src/resummation_jet_momenta.cc
+++ b/src/resummation_jet_momenta.cc
@@ -1,173 +1,218 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <vector>
#include <array>
#include <algorithm>
#include "RHEJ/resummation_jet_momenta.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/gsl_wrapper.hh"
namespace {
using namespace gsl;
enum Coordinate{
x = 0,
y = 1,
};
struct BornParameters{
std::vector< std::array<double, 2> > pt_born;
std::array<double, 2> q;
};
double pt_abs(double px, double py){
return sqrt(px*px + py*py);
}
//! calculate momentum difference for jets
/**
* After reshuffling, we have the condition
* Delta = pt_born - pt_res + q* |pt_res|/P_perp = 0
* for each jet, where P_perp is the sum of |pt_res| over all jets
* This function calculates Delta and stores the result in diff
*
* Since the gsl_vectors are one-dimensional, indices are a bit funny;
* diff->data[2*i + X]
* corresponds to the X component of the momentum vector for jet i
*/
int calc_jet_diff(const gsl_vector * pt_resum, void *data, gsl_vector * diff){
assert(diff->size == pt_resum->size);
assert(diff->size % 2 == 0);
auto param = static_cast<BornParameters const *>(data);
auto const & pt_born = param->pt_born;
auto pt_res = pt_resum->data;
auto const & q = param->q;
const size_t n_jets = pt_resum->size/2;
double P_perp = 0.;
for(size_t jet = 0; jet < n_jets; ++jet){
P_perp += pt_abs(pt_res[2*jet + x], pt_res[2*jet + y]);
}
for(size_t jet = 0; jet < n_jets; ++jet){
const double pt_res_abs = pt_abs(pt_res[2*jet + x], pt_res[2*jet + y]);
for(Coordinate X: {x,y}){
diff->data[2*jet + X] =
pt_born[jet][X] - pt_res[2*jet + X] - q[X]*pt_res_abs/P_perp;
}
}
return GSL_SUCCESS;
}
// computes resummation jet pt from Born jet pt
// if this fails an empty vector is returned
std::vector< std::array<double ,2> > reshuffle(BornParameters const & p){
constexpr int max_iterations = 1000;
const size_t n_jets = p.pt_born.size();
gsl_multiroot_function f = {
&calc_jet_diff, 2*n_jets,
const_cast<BornParameters *>(&p)
};
Vector pt_resum{2*n_jets};
// initial values for solver - resummation pt = Born pt
for (size_t jet = 0; jet < n_jets; ++jet) {
pt_resum[2*jet+x] = p.pt_born[jet][x];
pt_resum[2*jet+y] = p.pt_born[jet][y];
}
MultirootFsolver solver{
gsl_multiroot_fsolver_hybrids,
&f, std::move(pt_resum)
};
// solve equations
int status = GSL_CONTINUE;
for(
int iterations = 1;
status == GSL_CONTINUE && iterations < max_iterations;
++iterations
){
status = solver.iterate();
if(status == GSL_EBADFUNC || status == GSL_ENOPROG) return {};
status = solver.test_residual(1e-7);
}
if(status != GSL_SUCCESS) return {};
std::vector< std::array<double ,2> > res_pt(n_jets);
for(size_t jet = 0; jet < n_jets; ++jet) {
res_pt[jet][x] = gsl_vector_get(solver.x(), 2*jet + x);
res_pt[jet][y] = gsl_vector_get(solver.x(), 2*jet + y);
}
return res_pt;
}
// check that pt_i^B == pt_i + qt_i for each jet
void assert_pt_conservation(
std::vector<fastjet::PseudoJet> const & jetvects,
fastjet::PseudoJet const & qperp,
std::vector<fastjet::PseudoJet> const & shuffledmomenta
){
RHEJ::ignore(jetvects, qperp, shuffledmomenta);
#ifndef NDEBUG
assert(jetvects.size() == shuffledmomenta.size());
double Pperp = 0;
for(auto const & p: shuffledmomenta) Pperp += p.perp();
for(size_t i = 0; i < jetvects.size(); ++i){
const auto qperp_i = qperp*shuffledmomenta[i].perp()/Pperp;
const auto pdiff = jetvects[i] - shuffledmomenta[i] - qperp_i;
assert(RHEJ::nearby_ep(pdiff.px(), 0, 1e-5));
assert(RHEJ::nearby_ep(pdiff.py(), 0, 1e-5));
}
#endif
}
-}
-
-namespace RHEJ{
- std::vector<fastjet::PseudoJet> resummation_jet_momenta(
+ // for "old" reshuffling p^B = p + q*|p|/P
+ std::vector<fastjet::PseudoJet> jet_momenta_linear_in_born(
std::vector<fastjet::PseudoJet> const & p_born,
fastjet::PseudoJet const & q
) {
std::vector<fastjet::PseudoJet> p_res;
p_res.reserve(p_born.size());
BornParameters r;
r.q[x] = q.px();
r.q[y] = q.py();
r.pt_born.resize(p_born.size());
for (size_t jet = 0; jet < p_born.size(); ++jet) {
r.pt_born[jet][x] = p_born[jet].px();
r.pt_born[jet][y] = p_born[jet].py();
}
const auto pt_reshuffled = reshuffle(r);
if(pt_reshuffled.empty()) return {};
// Construct the new 4-momenta
for (size_t jet = 0; jet < p_born.size(); ++jet) {
const double px = pt_reshuffled[jet][x];
const double py = pt_reshuffled[jet][y];
const double pperp = sqrt(px*px + py*py);
// keep the rapidities fixed
const double pz = pperp*sinh(p_born[jet].rapidity());
const double E = pperp*cosh(p_born[jet].rapidity());
p_res.emplace_back(px, py, pz, E);
assert(
- nearby_ep(
+ RHEJ::nearby_ep(
p_res.back().rapidity(),
p_born[jet].rapidity(),
1e-5
)
);
}
assert_pt_conservation(p_born, q, p_res);
return p_res;
}
+
+ // for "new" reshuffling p^B = p + q*|p^B|/P^B
+ std::vector<fastjet::PseudoJet> jet_momenta_linear_in_resummed(
+ std::vector<fastjet::PseudoJet> const & p_born,
+ fastjet::PseudoJet const & q
+ ) {
+ double Pperp_born = 0.;
+ for(auto const & p: p_born) Pperp_born += p.perp();
+
+ std::vector<fastjet::PseudoJet> p_res;
+ p_res.reserve(p_born.size());
+ for(auto & pB: p_born) {
+ const double px = pB.px() - q.px()*pB.perp()/Pperp_born;
+ const double py = pB.py() - q.py()*pB.perp()/Pperp_born;
+ const double pperp = sqrt(px*px + py*py);
+ // keep the rapidities fixed
+ const double pz = pperp*sinh(pB.rapidity());
+ const double E = pperp*cosh(pB.rapidity());
+ p_res.emplace_back(px, py, pz, E);
+ assert(
+ RHEJ::nearby_ep(
+ p_res.back().rapidity(),
+ pB.rapidity(),
+ 1e-5
+ )
+ );
+ }
+ return p_res;
+ }
+}
+
+namespace RHEJ{
+
+ std::vector<fastjet::PseudoJet> resummation_jet_momenta(
+ std::vector<fastjet::PseudoJet> const & p_born,
+ fastjet::PseudoJet const & q,
+ ReshufflingAlgorithm algo
+ ) {
+ switch(algo) {
+ case ReshufflingAlgorithm::linear_in_born:
+ return jet_momenta_linear_in_born(p_born, q);
+ case ReshufflingAlgorithm::linear_in_resummed:
+ return jet_momenta_linear_in_resummed(p_born, q);
+ default:;
+ }
+ throw std::logic_error{"unreachable"};
+ }
}