Page MenuHomeHEPForge
Files F19252187

resummation_jet.cc
No OneTemporary
Actions

Size
3 KB
Referenced Files
None
Subscribers
None

resummation_jet.cc
View Options

/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/resummation_jet.hh"
#include <cassert>
#include <cmath>
#include <cstddef>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include "fastjet/PseudoJet.hh"
#include "HEJ/utility.hh"
namespace HEJ {
std::vector<fastjet::PseudoJet> resummation_jet_momenta(
std::vector<fastjet::PseudoJet const *> const & p_born,
fastjet::PseudoJet const & qperp
) {
// for "new" reshuffling p^B = p + qperp*|p^B|/P^B
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 const & pB: p_born) {
const double px = pB->px() - qperp.px()*pB->perp()/Pperp_born;
const double py = pB->py() - qperp.py()*pB->perp()/Pperp_born;
const double mperp = std::sqrt(px*px + py*py + pB->m2());
// keep the rapidities fixed
const double pz = mperp*std::sinh(pB->rapidity());
const double E = mperp*std::cosh(pB->rapidity());
p_res.emplace_back(px, py, pz, E);
assert(
nearby_ep(
p_res.back().rapidity(),
pB->rapidity(),
1e-5
)
);
}
return p_res;
}
namespace {
enum coordinates : std::size_t {
x1, x2
};
namespace ublas = boost::numeric::ublas;
template<class Matrix>
double det(ublas::matrix_expression<Matrix> const& m) {
ublas::permutation_matrix<std::size_t> pivots{m().size1()};
Matrix mLu{m()};
const auto is_singular = lu_factorize(mLu, pivots);
if(is_singular) return 0.;
double det = 1.0;
for (std::size_t i = 0; i < pivots.size(); ++i){
if (pivots(i) != i) det = -det;
det *= mLu(i,i);
}
return det;
}
using ublas::matrix;
} // namespace
double resummation_jet_weight(
std::vector<fastjet::PseudoJet const *> const & p_born,
fastjet::PseudoJet const & qperp
) {
using std::size_t;
static constexpr int num_coordinates = 2;
auto Jacobian = matrix<double>{
num_coordinates*p_born.size(),
num_coordinates*p_born.size()
};
double P_perp = 0.;
for(auto const & J: p_born) P_perp += J->perp();
for(size_t l = 0; l < p_born.size(); ++l){
const double Jl_perp = p_born[l]->perp();
for(size_t lp = 0; lp < p_born.size(); ++lp){
const int delta_l = static_cast<int>(l == lp);
auto const & Jlp = p_born[lp];
const double Jlp_perp = Jlp->perp();
for(size_t x = x1; x <= x2; ++x){
const double qxy = (x==x1)?qperp.px():qperp.py();
for(size_t xp = x1; xp <= x2; ++xp){
const double Jxy = (xp==x1)?Jlp->px():Jlp->py();
const int delta_x = static_cast<int>(x == xp);
Jacobian(2*l + x, 2*lp + xp) =
+ delta_l*delta_x
- qxy*Jxy/(P_perp*Jlp_perp) * (delta_l - Jl_perp/P_perp);
}
}
}
}
return det(Jacobian);
}
} // namespace HEJ

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 6:15 AM (5 h, 46 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6475914
Default Alt Text
resummation_jet.cc (3 KB)

Event Timeline

HEPForge.org · .