Page MenuHomeHEPForge
Files F8724393

JetShape.cc
No OneTemporary
Actions

JetShape.cc
View Options

// -*- C++ -*-
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/JetShape.hh"
namespace Rivet {
// Constructor.
JetShape::JetShape(const JetAlg& jetalg,
double rmin, double rmax, size_t nbins,
double ptmin, double ptmax,
double absrapmin, double absrapmax,
RapScheme rapscheme)
: _rapscheme(rapscheme)
{
setName("JetShape");
_binedges = linspace(nbins, rmin, rmax);
_ptcuts = make_pair(ptmin, ptmax);
_rapcuts = make_pair(absrapmin, absrapmax);
addProjection(jetalg, "Jets");
}
// Constructor.
JetShape::JetShape(const JetAlg& jetalg,
vector<double> binedges,
double ptmin, double ptmax,
double absrapmin, double absrapmax,
RapScheme rapscheme)
: _binedges(binedges), _rapscheme(rapscheme)
{
setName("JetShape");
_ptcuts = make_pair(ptmin, ptmax);
_rapcuts = make_pair(absrapmin, absrapmax);
addProjection(jetalg, "Jets");
}
int JetShape::compare(const Projection& p) const {
const int jcmp = mkNamedPCmp(p, "Jets");
if (jcmp != EQUIVALENT) return jcmp;
const JetShape& other = pcast<JetShape>(p);
const int ptcmp = cmp(ptMin(), other.ptMin()) || cmp(ptMax(), other.ptMax());
if (ptcmp != EQUIVALENT) return ptcmp;
const int rapcmp = cmp(_rapcuts.first, other._rapcuts.first) || cmp(_rapcuts.second, other._rapcuts.second);
if (rapcmp != EQUIVALENT) return rapcmp;
int bincmp = cmp(numBins(), other.numBins());
if (bincmp != EQUIVALENT) return bincmp;
for (size_t i = 0; i < _binedges.size(); ++i) {
bincmp = cmp(_binedges[i], other._binedges[i]);
if (bincmp != EQUIVALENT) return bincmp;
}
return EQUIVALENT;
}
void JetShape::clear() {
_diffjetshapes.clear();
}
void JetShape::calc(const Jets& jets) {
clear();
foreach (const Jet& j, jets) {
// Apply jet cuts
FourMomentum pj = j.momentum();
if (!inRange(pj.pT(), _ptcuts)) continue;
/// @todo Introduce a better (i.e. more safe and general) eta/y selection mechanism: MomentumFilter
if (_rapscheme == PSEUDORAPIDITY && !inRange(fabs(pj.eta()), _rapcuts)) continue;
if (_rapscheme == RAPIDITY && !inRange(fabs(pj.rapidity()), _rapcuts)) continue;
// Fill bins
vector<double> bins(numBins(), 0.0);
foreach (const Particle& p, j.particles()) {
const double dR = deltaR(pj, p.momentum(), _rapscheme);
const int dRindex = index_between(dR, _binedges);
if (dRindex == -1) continue; //< Out of histo range
bins[dRindex] += p.momentum().pT();
}
// Add bin vector for this jet to the diffjetshapes container
_diffjetshapes += bins;
}
// Normalize to total pT
foreach (vector<double>& binsref, _diffjetshapes) {
double integral = 0.0;
for (size_t i = 0; i < numBins(); ++i) {
integral += binsref[i];
}
if (integral > 0) {
for (size_t i = 0; i < numBins(); ++i) {
binsref[i] /= integral;
}
} else {
// It's just-about conceivable that a jet would have no particles in the given Delta(r) range...
MSG_DEBUG("No pT contributions in jet Delta(r) range: weird!");
}
}
}
void JetShape::project(const Event& e) {
const Jets jets = applyProjection<JetAlg>(e, "Jets").jets(_ptcuts.first, _ptcuts.second,
-_rapcuts.second, _rapcuts.second, _rapscheme);
calc(jets);
}
}

File Metadata

Mime Type
text/x-c
Expires
Mon, Jan 20, 11:19 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4239999
Default Alt Text
JetShape.cc (3 KB)

Event Timeline