Page MenuHomeHEPForge
Files F9501757

No OneTemporary
Actions

View Options

diff --git a/src/resummation_jet.cc b/src/resummation_jet.cc
index f51ba8d..9f94c1d 100644
--- a/src/resummation_jet.cc
+++ b/src/resummation_jet.cc
@@ -1,112 +1,113 @@
/**
* \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();
+ const int delta_l = (l == lp);
+ const auto & 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 = 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
- );
+ - qxy*Jxy/(P_perp*Jlp_perp) * (delta_l - Jl_perp/P_perp);
}
}
}
}
return det(Jacobian);
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 3:02 PM (5 h, 49 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486806
Default Alt Text
(3 KB)

Event Timeline