Page MenuHomeHEPForge
Files F19251167

JetSplitter.cc
No OneTemporary
Actions

Size
6 KB
Referenced Files
None
Subscribers
None

JetSplitter.cc
View Options

/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/JetSplitter.hh"
#include <array>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <numeric>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/PseudoJet.hh"
#include "HEJ/Constants.hh"
#include "HEJ/RNG.hh"
#include "HEJ/exceptions.hh"
namespace HEJ {
namespace {
template<class Iterator>
bool same_pt_and_rapidity(
Iterator begin, Iterator end,
fastjet::PseudoJet const & jet
){
constexpr double ep = 1e-2;
const fastjet::PseudoJet reconstructed_jet = std::accumulate(
begin, end, fastjet::PseudoJet{}
);
return
(std::abs(reconstructed_jet.pt() - jet.pt()) < ep)
&& (std::abs(reconstructed_jet.rapidity() - jet.rapidity()) < ep)
;
}
bool all_in_one_jet(
std::vector<fastjet::PseudoJet> const & partons,
fastjet::JetDefinition const & jet_def, double min_jet_pt
){
fastjet::ClusterSequence ev(partons, jet_def);
const std::vector<fastjet::PseudoJet> testjet = ev.inclusive_jets(min_jet_pt);
return testjet.size() == 1u
&& testjet[0].constituents().size() == partons.size();
}
/** @brief sample y-phi distance to jet pt axis for a jet splitting into two
* partons
*
* @param wt Multiplied by the weight of the sampling point
* @param ran Random number generator
* @returns The distance in units of the jet radius
*/
double sample_distance_2p(double & wt, RNG & ran) {
static constexpr double x_small = 0.1;
static constexpr double p_small = 0.4;
const double pR = ran.flat();
if(pR < p_small){
wt *= x_small/p_small;
return x_small/p_small*pR;
}
wt *= (1-x_small)/(1-p_small);
return (1-x_small)/(1-p_small)*(pR-p_small) + x_small;
}
} // namespace
JetSplitter::JetSplitter(
fastjet::JetDefinition jet_def, double min_pt
):
min_jet_pt_{min_pt},
jet_def_{std::move(jet_def)}
{}
JetSplitter::SplitResult JetSplitter::split(
fastjet::PseudoJet const & j2split, int ncons, RNG & ran
) const{
if(ncons <= 0) {
throw std::invalid_argument{
"number of requested jet constituents less than 1"
};
}
double swt = 1.;
std::vector<fastjet::PseudoJet> jcons;
if(ncons == 1){
jcons.emplace_back(j2split);
jcons.back().set_user_index(0);
return {jcons, swt};
}
if(ncons == 2){
return Split2(j2split, ran);
}
const double R_max = R_FACTOR*jet_def_.R();
assert(R_max < M_PI);
double pt_remaining = j2split.pt();
const double phi_jet = j2split.phi();
const double y_jet = j2split.rapidity();
for(int i = 0; i < ncons - 1; ++i){
/**
* Generate rapidity and azimuthal angle with a distance
* R = sqrt(delta_y^2 + delta_phi^2) < R_max
* from the jet centre
*/
const double R = R_max*ran.flat();
const double theta = 2*M_PI*ran.flat();
const double delta_phi = R*std::cos(theta);
const double delta_y = R*std::sin(theta);
/**
* Generate pt such that the total contribution of all partons
* along the jet pt axis does not exceed the jet pt
*/
const double pt_max = pt_remaining/std::cos(delta_phi);
assert(pt_max > 0);
if(pt_max < CMINPT) return {}; // no pt remaining for this parton
const double pt = (pt_max - CMINPT)*ran.flat() + CMINPT;
pt_remaining -= pt*std::cos(delta_phi);
jcons.emplace_back(
pt*std::cos(phi_jet + delta_phi), pt*std::sin(phi_jet + delta_phi),
pt*std::sinh(y_jet + delta_y), pt*std::cosh(y_jet + delta_y)
);
jcons.back().set_user_index(i);
swt *= 2*M_PI*R*R_max*pt*(pt_max - CMINPT);
}
const fastjet::PseudoJet p_total = std::accumulate(
jcons.cbegin(), jcons.cend(), fastjet::PseudoJet{}
);
// Calculate the pt of the last parton
const double last_px = j2split.px() - p_total.px();
const double last_py = j2split.py() - p_total.py();
const double last_pt = std::sqrt(last_px*last_px + last_py*last_py);
if(last_pt < CMINPT) return {};
// Calculate the rapidity of the last parton using the requirement that the
// new jet must have the same rapidity as the LO jet.
const double exp_2y_jet = (j2split.e() + j2split.pz())/(j2split.e() - j2split.pz());
const double bb = (p_total.e()+p_total.pz()) - exp_2y_jet*(p_total.e()-p_total.pz());
const double lasty = std::log(
( -bb+std::sqrt(bb*bb+4.*exp_2y_jet*last_pt*last_pt) ) / ( 2.*last_pt ));
jcons.emplace_back(
last_px, last_py, last_pt*std::sinh(lasty), last_pt*std::cosh(lasty)
);
jcons.back().set_user_index(ncons-1);
assert(same_pt_and_rapidity(jcons.cbegin(), jcons.cend(), j2split));
// Test that the last parton is not too far away from the jet centre.
if (jcons.back().delta_R(j2split) > R_max) return {};
if(! all_in_one_jet(jcons, jet_def_, min_jet_pt_)) return {};
return {jcons, swt};
}
JetSplitter::SplitResult JetSplitter::Split2(
fastjet::PseudoJet const & j2split, RNG & ran
) const {
static constexpr std::size_t ncons = 2;
std::vector<fastjet::PseudoJet> jcons(ncons);
std::array<double, ncons> R{};
std::array<double, ncons> phi{};
std::array<double, ncons> y{};
std::array<double, ncons> pt{};
double wt = 1;
const double theta = 2*M_PI*ran.flat(); // angle in y-phi plane
// empiric observation: we are always within the jet radius
R[0] = sample_distance_2p(wt, ran)*jet_def_.R();
R[1] = -sample_distance_2p(wt, ran)*jet_def_.R();
for(std::size_t i = 0; i <= 1; ++i){
phi[i] = j2split.phi() + R[i]*std::cos(theta);
y[i] = j2split.rapidity() + R[i]*std::sin(theta);
}
for(std::size_t i = 0; i <= 1; ++i){
pt[i] = (j2split.py() - std::tan(phi[1-i])*j2split.px())/
(std::sin(phi[i]) - std::tan(phi[1-i])*std::cos(phi[i]));
if(pt[i] < CMINPT) return {};
jcons[i].reset_PtYPhiM(pt[i], y[i], phi[i]);
jcons[i].set_user_index(i);
}
assert(same_pt_and_rapidity(jcons.cbegin(), jcons.cend(), j2split));
if(! all_in_one_jet(jcons, jet_def_, min_jet_pt_)) return {};
wt *= 2*M_PI*pt[0]*R[0]*jet_def_.R()*jet_def_.R();
// from transformation of delta(R[1] - ...) to delta(pt[0] - ...)
const double dphi0 = phi[0] - j2split.phi();
const double ptJ = j2split.pt();
const double jacobian = std::cos(theta)*pt[1]*pt[1]/(ptJ*std::sin(dphi0));
return {jcons, jacobian*wt};
}
} // namespace HEJ

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 5:48 AM (1 d, 10 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6480832
Default Alt Text
JetSplitter.cc (6 KB)

Event Timeline

HEPForge.org · .