2
#ifndef RIVET_FastJets_HH
3
#define RIVET_FastJets_HH
5
#include "Rivet/Projection.hh"
6
#include "Rivet/Projections/JetAlg.hh"
7
#include "Rivet/Projections/FinalState.hh"
8
#include "Rivet/Particle.hh"
9
#include "Rivet/Jet.hh"
11
#include "fastjet/JetDefinition.hh"
12
#include "fastjet/AreaDefinition.hh"
13
#include "fastjet/ClusterSequence.hh"
14
#include "fastjet/ClusterSequenceArea.hh"
15
#include "fastjet/PseudoJet.hh"
17
#include "fastjet/SISConePlugin.hh"
18
#include "fastjet/ATLASConePlugin.hh"
19
#include "fastjet/CMSIterativeConePlugin.hh"
20
#include "fastjet/CDFJetCluPlugin.hh"
21
#include "fastjet/CDFMidPointPlugin.hh"
22
#include "fastjet/D0RunIIConePlugin.hh"
23
#include "fastjet/TrackJetPlugin.hh"
24
#include "fastjet/JadePlugin.hh"
25
//#include "fastjet/PxConePlugin.hh"
30
/// Make a 3-momentum vector from a FastJet pseudo-jet
31
inline Vector3 momentum3(const fastjet::PseudoJet& pj) {
32
return Vector3(pj.px(), pj.py(), pj.pz());
35
/// Make a 4-momentum vector from a FastJet pseudo-jet
36
inline FourMomentum momentum(const fastjet::PseudoJet& pj) {
37
return FourMomentum(pj.E(), pj.px(), pj.py(), pj.pz());
41
/// Typedef for a collection of PseudoJets.
42
typedef vector<fastjet::PseudoJet> PseudoJets;
46
/////////////////////////
50
/// Project out jets found using the FastJet package jet algorithms.
51
class FastJets : public JetAlg {
54
/// Wrapper enum for selected Fastjet jet algorithms.
55
enum JetAlgName { KT, CAM, SISCONE, ANTIKT,
58
CDFJETCLU, CDFMIDPOINT, D0ILCONE,
59
JADE, DURHAM, TRACKJET };
62
/// @name Constructors etc.
65
/// "Wrapped" argument constructor using Rivet enums for most common
66
/// jet alg choices (including some plugins). For the built-in algs,
67
/// E-scheme recombination is used. For full control of
68
/// FastJet built-in jet algs, use the native arg constructor.
69
FastJets(const FinalState& fsp, JetAlgName alg,
70
double rparameter, double seed_threshold=1.0)
71
: JetAlg(fsp), _adef(0) { _init1(alg, rparameter, seed_threshold); }
73
/// Native argument constructor, using FastJet alg/scheme enums.
74
FastJets(const FinalState& fsp, fastjet::JetAlgorithm type,
75
fastjet::RecombinationScheme recom, double rparameter)
76
: JetAlg(fsp), _adef(0) { _init2(type, recom, rparameter); }
78
/// Explicitly pass in an externally-constructed plugin (must be heap-allocated, Rivet will delete)
79
FastJets(const FinalState& fsp, fastjet::JetDefinition::Plugin* plugin)
80
: JetAlg(fsp), _adef(0) { _init3(plugin); }
81
/// Explicitly pass in an externally-constructed plugin (must be heap-allocated, Rivet will delete)
82
FastJets(const FinalState& fsp, fastjet::JetDefinition::Plugin& plugin)
83
: JetAlg(fsp), _adef(0) { _init3(&plugin); }
86
/// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
87
FastJets(JetAlgName alg, double rparameter, double seed_threshold=1.0)
88
: _adef(0) { _init1(alg, rparameter, seed_threshold); }
89
/// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
90
FastJets(fastjet::JetAlgorithm type, fastjet::RecombinationScheme recom, double rparameter)
91
: _adef(0) { _init2(type, recom, rparameter); }
92
/// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
93
FastJets(fastjet::JetDefinition::Plugin* plugin)
94
: _adef(0) { _init3(plugin); }
95
/// Same thing as above, but without an FS (for when we want to pass the particles directly to the calc method)
96
FastJets(fastjet::JetDefinition::Plugin& plugin)
97
: _adef(0) { _init3(&plugin); }
100
/// Clone on the heap.
101
virtual const Projection* clone() const {
102
return new FastJets(*this);
110
/// Reset the projection. Jet def, etc. are unchanged.
113
/// @brief Use provided jet area definition
114
/// @warning adef is NOT copied, the user must ensure that it remains valid!
115
/// Provide an adef null pointer to re-disable jet area calculation
116
void useJetArea(fastjet::AreaDefinition* adef) {
120
/// Number of jets above the \f$ p_\perp \f$ cut.
121
size_t numJets(double ptmin = 0.0) const;
124
size_t size() const {
128
/// Get the jets (unordered).
129
Jets _jets(double ptmin = 0.0) const;
131
/// Get the pseudo jets (unordered).
132
PseudoJets pseudoJets(double ptmin = 0.0) const;
134
/// Get the pseudo jets, ordered by \f$ p_T \f$.
135
PseudoJets pseudoJetsByPt(double ptmin = 0.0) const {
136
return sorted_by_pt(pseudoJets(ptmin));
139
/// Get the pseudo jets, ordered by \f$ E \f$.
140
PseudoJets pseudoJetsByE(double ptmin = 0.0) const {
141
return sorted_by_E(pseudoJets(ptmin));
144
/// Get the pseudo jets, ordered by rapidity.
145
PseudoJets pseudoJetsByRapidity(double ptmin = 0.0) const {
146
return sorted_by_rapidity(pseudoJets(ptmin));
149
/// Return the cluster sequence (FastJet-specific).
150
const fastjet::ClusterSequence* clusterSeq() const {
154
/// Return the cluster sequence (FastJet-specific).
155
const fastjet::ClusterSequenceArea* clusterSeqArea() const {
156
/// @todo Throw error if no area def? Or just blindly call dynamic_cast?
157
if (_adef == 0) return (fastjet::ClusterSequenceArea*) 0;
158
return dynamic_cast<fastjet::ClusterSequenceArea*>(_cseq.get());
161
/// Return the jet definition (FastJet-specific).
162
const fastjet::JetDefinition& jetDef() const {
166
/// Return the area definition (FastJet-specific). May be null.
167
const fastjet::AreaDefinition* areaDef() const {
171
/// Get the subjet splitting variables for the given jet.
172
vector<double> ySubJet(const fastjet::PseudoJet& jet) const;
174
/// @brief Split a jet a la PRL100,242001(2008).
175
/// Based on code from G.Salam, A.Davison.
176
fastjet::PseudoJet splitJet(fastjet::PseudoJet jet, double& last_R) const;
178
/// @brief Filter a jet a la PRL100,242001(2008).
179
/// Based on code from G.Salam, A.Davison.
180
fastjet::PseudoJet filterJet(fastjet::PseudoJet jet, double& stingy_R, const double def_R) const;
184
Jets _pseudojetsToJets(const PseudoJets& pjets) const;
186
/// Shared utility functions to implement constructor behaviour
187
void _init1(JetAlgName alg, double rparameter, double seed_threshold);
188
void _init2(fastjet::JetAlgorithm type,
189
fastjet::RecombinationScheme recom, double rparameter);
190
void _init3(fastjet::JetDefinition::Plugin* plugin);
194
/// Perform the projection on the Event.
195
void project(const Event& e);
197
/// Compare projections.
198
int compare(const Projection& p) const;
202
/// Do the calculation locally (no caching).
203
void calc(const ParticleVector& ps);
209
fastjet::JetDefinition _jdef;
211
/// Pointer to user-handled area definition
212
fastjet::AreaDefinition* _adef;
215
shared_ptr<fastjet::ClusterSequence> _cseq;
217
/// FastJet external plugin
218
shared_ptr<fastjet::JetDefinition::Plugin> _plugin;
220
/// Map of vectors of y scales. This is mutable so we can use caching/lazy evaluation.
221
mutable map<int, vector<double> > _yscales;
223
/// set of particles sorted by their PT2
224
//set<Particle, ParticleBase::byPTAscending> _particles;
225
map<int, Particle> _particles;