Page MenuHomeHEPForge
Files F19251802

resummation_jet.cc
No OneTemporary
Actions

Size
2 KB
Referenced Files
None
Subscribers
None

resummation_jet.cc
View Options

/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/resummation_jet.hh"
#include <assert.h>
#include <math.h>
#include <stdio.h>
#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 & 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 & 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 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(
HEJ::nearby_ep(
p_res.back().rapidity(),
pB.rapidity(),
1e-5
)
);
}
return p_res;
}
namespace{
enum coordinates : size_t {
x1, x2
};
namespace ublas = boost::numeric::ublas;
template<class Matrix>
double det(ublas::matrix_expression<Matrix> const& m) {
ublas::permutation_matrix<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;
}
double resummation_jet_weight(
std::vector<fastjet::PseudoJet> const & p_born,
fastjet::PseudoJet const & qperp
) {
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 = l == lp;
const double Jlp_perp = p_born[lp].perp();
for(size_t x = x1; x <= x2; ++x){
for(size_t xp = x1; xp <= x2; ++xp){
const int delta_x = x == xp;
Jacobian(2*l + x, 2*lp + xp) =
+ delta_l*delta_x
- qperp[x]*p_born[lp][xp]/(P_perp*Jlp_perp)*(
+ delta_l - Jl_perp/P_perp
);
}
}
}
}
return det(Jacobian);
}
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 6:10 AM (1 d, 19 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6469446
Default Alt Text
resummation_jet.cc (2 KB)

Event Timeline

HEPForge.org · .