Page MenuHomeHEPForge
Files F8723730

No OneTemporary
Actions

View Options

diff --git a/src/resummation_jet_momenta.cc b/src/resummation_jet_momenta.cc
index 8db9bf8..c2ff1c5 100644
--- a/src/resummation_jet_momenta.cc
+++ b/src/resummation_jet_momenta.cc
@@ -1,174 +1,173 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <vector>
#include <array>
#include <algorithm>
-// #include "RHEJ/gsl_wrapper.hh"
#include "RHEJ/resummation_jet_momenta.hh"
#include "RHEJ/debug.hh"
#include "RHEJ/gsl_wrapper.hh"
-using namespace gsl;
-namespace RHEJ{
+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;
}
- namespace{
-
// 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
){
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(nearby_ep(pdiff.px(), 0, 1e-5));
assert(nearby_ep(pdiff.py(), 0, 1e-5));
}
#endif
}
- }
+}
+
+namespace RHEJ{
std::vector<fastjet::PseudoJet> resummation_jet_momenta(
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(
p_res.back().rapidity(),
p_born[jet].rapidity(),
1e-5
)
);
}
assert_pt_conservation(p_born, q, p_res);
return p_res;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:14 PM (21 h, 54 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242432
Default Alt Text
(5 KB)

Event Timeline