Loading...
Searching...
No Matches
fastjet_timing_plugins.cc
//STARTHEADER
// $Id$
//
// Copyright (c) 2005-2018, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development and are described in hep-ph/0512210. If you use
// FastJet as part of work towards a scientific publication, please
// include a citation to the FastJet paper.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//ENDHEADER
//----------------------------------------------------------------------
/// fastjet_timing.cc: Program to help time and test the fastjet package
///
/// It reads files containing multiple events in the format
/// p1x p1y p1z E1
/// p2x p2y p2z E2
/// ...
/// #END
///
/// An example input file containing 10 events is included as
/// data/Pythia-PtMin1000-LHC-10ev.dat
///
/// Usage:
/// fastjet_timing [-strategy NUMBER] [-repeat nrepeats] [-massive] \
/// [-combine nevents] [-r Rparameter] [-incl ptmin] [...] \
/// < data_file
///
/// where the clustering can be repeated to aid timing and multiple
/// events can be combined to get to larger multiplicities. Some options:
///
/// Options for reading
/// -------------------
///
/// -nev n number of events to run
///
/// -combine n for combining multiple events from the data file in order
/// to get a single high-multipicity event to run.
///
/// -massless read in only the 3-momenta and deduce energies assuming
/// that particles are massless
///
/// -dense adds dense ghost coverage
///
/// -repeat n repeats each event n times
///
/// -nhardest n keep only the n hardest particles in the event
///
/// -file name read from the corresponding file rather than stdin.
/// (The file will be reopened for each new jet alg.; in
/// constrast, if you use stdin, each new alg will take a
/// new event).
///
/// Output Options
/// --------------
///
/// -incl ptmin output of all inclusive jets with pt > ptmin is obtained
/// with the -incl option.
///
/// -repeat-incl ptmin
/// same as -incl ptmin but do it for each repetition
/// of the clustering
///
/// -excld dcut output of all exclusive jets as obtained in a clustering
/// with dcut
///
/// -excly ycut output of all exclusive jets as obtained in a clustering
/// with ycut
///
/// -excln n output of clustering to n exclusive jets
///
/// -ee-print print things as px,py,pz,E
///
/// -get-all-dij print out all dij values
/// -get-all-yij print out all yij values
///
/// -const show jet constituents (works with excl jets)
///
/// -write for writing out detailed clustering sequence (valuable
/// for testing purposes)
///
/// -unique_write writes out the sequence of dij's according to the
/// "unique_history_order" (useful for verifying consistency
/// between different clustering strategies).
///
/// -root file sends output to file that can be read in with the script in
/// root/ so as to show a lego-plot of the event
///
/// -cones show extra info about internal steps for SISCone
///
/// -area calculate areas. Additional options include
/// -area:active
/// -area:passive
/// -area:explicit
/// -area:voronoi Rfact
/// -area:repeat nrepeat
/// -ghost-area area
/// -ghost-maxrap maxrap
/// -area:fj2 place ghosts as in fj2
///
/// -bkgd calculate the background density. Additional options include
/// -bkgd:csab use the old ClusterSequenceAreaBase methods
/// -bkgd:jetmedian use the new JetMedianBackgroundEstimator class
/// -bkgd:fj2 force jetmedian to calculate sigma as in fj2
/// -bkgd:gridmedian use GridMedianBackgroundEstimator with grid up to ghost_maxrap-ktR and grid spacing of 2ktR
///
/// -compare-strategy STRAT
/// compares the output of the default strategy (possibly as specified
/// with -strategy) with that from STRAT. Currently compares the history.
/// Algorithms
/// ----------
/// -all-algs runs all algorithms
///
/// -kt switch to the longitudinally invariant kt algorithm
/// Note: this is the default one.
///
/// -cam switch to the inclusive Cambridge/Aachen algorithm --
/// note that the option -excld dcut provides a clustering
/// up to the dcut which is the minimum squared
/// distance between any pair of jets.
///
/// -antikt switch to the anti-kt clustering algorithm
///
/// -genkt switch to the genkt algorithm
/// you can provide the parameter of the alg as an argument to
/// -genkt (1 by default)
///
/// -eekt switch to the e+e- kt algorithm
///
/// -eegenkt switch to the genkt algorithm
/// you can provide the parameter of the alg as an argument to
/// -ee_genkt (1 by default)
///
/// plugins (don't delete this line)
///
/// -pxcone switch to the PxCone jet algorithm
///
/// -siscone switch to the SISCone jet algorithm (seedless cones)
/// -sisconespheri switch to the Spherical SISCone jet algorithm (seedless cones)
///
/// -midpoint switch to CDF's midpoint code
/// -jetclu switch to CDF's jetclu code
///
/// -d0runipre96cone switch to the D0RunIpre96Cone plugin
/// -d0runicone switch to the D0RunICone plugin
///
/// -d0runiicone switch to D0's run II midpoint cone
///
/// -trackjet switch to the TrackJet plugin
///
/// -atlascone switch to the ATLASCone plugin
///
/// -eecambridge switch to the EECambridge plugin
///
/// -jade switch to the Jade plugin
///
/// -cmsiterativecone switch to the CMSIterativeCone plugin
///
/// -gridjet switch to the GridJet plugin
///
/// end of plugins (don't delete this line)
///
///
/// Options for running algs
/// ------------------------
///
/// -r sets the radius of the jet algorithm (default = 1.0)
///
/// -overlap | -f sets the overlap fraction in cone algs with split-merge
///
/// -seed sets the seed threshold
///
/// -strategy N indicate stratgey from the enum fastjet::Strategy (see
/// fastjet/JetDefinition.hh).
///
#ifndef __FJCORE__
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
#include "fastjet/tools/GridMedianBackgroundEstimator.hh"
#include "fastjet/tools/Subtractor.hh"
#include "fastjet/Selector.hh"
#else
#include "fjcore.hh"
#endif
#include<iostream>
#include<iomanip>
#include<sstream>
#include<fstream>
#include<valarray>
#include<vector>
#include <cstdlib>
//#include<cstddef> // for size_t
#include "CmdLine.hh"
#ifndef __FJCORE__
// get info on how fastjet was configured
#include "fastjet/config.h"
#endif
// include the installed plugins (don't delete this line)
#ifdef FASTJET_ENABLE_PLUGIN_SISCONE
#include "fastjet/SISConePlugin.hh"
#include "fastjet/SISConeSphericalPlugin.hh"
#endif
#ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
#include "fastjet/CDFMidPointPlugin.hh"
#include "fastjet/CDFJetCluPlugin.hh"
#endif
#ifdef FASTJET_ENABLE_PLUGIN_PXCONE
#include "fastjet/PxConePlugin.hh"
#endif
#ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
#include "fastjet/D0RunIIConePlugin.hh"
#endif
#ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
#include "fastjet/TrackJetPlugin.hh"
#endif
#ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
#include "fastjet/ATLASConePlugin.hh"
#endif
#ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
#include "fastjet/EECambridgePlugin.hh"
#endif
#ifdef FASTJET_ENABLE_PLUGIN_JADE
#include "fastjet/JadePlugin.hh"
#endif
#ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
#include "fastjet/CMSIterativeConePlugin.hh"
#endif
#ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
#include "fastjet/D0RunIpre96ConePlugin.hh"
#include "fastjet/D0RunIConePlugin.hh"
#endif
#ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
#include "fastjet/GridJetPlugin.hh"
#endif
// end of installed plugins inclusion (don't delete this line)
using namespace std;
// to avoid excessive typing, use the fastjet namespace
#ifndef __FJCORE__
using namespace fastjet;
#else
using namespace fjcore;
#endif
inline double pow2(const double x) {return x*x;}
// pretty print the jets and their subjets
void print_jets_and_sub (const vector<PseudoJet> & jets, double dcut);
#ifndef __FJCORE__
void print_jets_bkgd(const vector<PseudoJet> &jets,
const vector<PseudoJet> &subtracted_jets,
BackgroundEstimatorBase * bge_ptr,
bool do_subtractor);
#endif // __FJCORE__
// have various kinds of subjet finding, to test consistency among them
//
// this is needed in print_jets_and_sub and declaring it in the
// function scope results in errors with older intel compilers (due to
// the overloaded == operator in PseudoJet which results in the "a
// template argument may not reference a local type" error)
enum SubType {subtype_internal, subtype_newclust_dcut, subtype_newclust_R};
void do_compare_strategy(int iev,
const vector<PseudoJet> & particles,
const JetDefinition & jet_def,
const ClusterSequence & cs,
int compare_strategy);
string rootfile;
CmdLine * cmdline_p;
bool do_areas;
/// sort and pretty print jets, with exact behaviour depending on
/// whether ee_print is true or not
bool ee_print = false;
void print_jets(const vector<PseudoJet> & jets, bool show_const = false);
bool found_unavailable = false;
void is_unavailable(const string & algname) {
cerr << algname << " requested, but not available for this compilation" << endl;
found_unavailable = true;
//exit(0);
}
/// a program to test and time a range of algorithms as implemented or
/// wrapped in fastjet
int main (int argc, char ** argv) {
CmdLine cmdline(argc,argv);
cout << "# " << cmdline.command_line() << endl;
ClusterSequence::print_banner();
cmdline_p = &cmdline;
// allow the use to specify the Strategy either through the
// -clever or the -strategy options (both will take numerical
// values); the latter will override the former.
Strategy strategy = Strategy(cmdline.int_val("-strategy",
cmdline.int_val("-clever", Best)));
int repeat = cmdline.int_val("-repeat",1);
int combine = cmdline.int_val("-combine",1);
bool write = cmdline.present("-write");
bool unique_write = cmdline.present("-unique_write");
bool hydjet = cmdline.present("-hydjet");
double ktR = cmdline.double_val("-r",1.0);
ktR = cmdline.double_val("-R",ktR); // allow -r and -R
double inclkt = cmdline.double_val("-incl",-1.0);
double repeatinclkt = cmdline.double_val("-repeat-incl",-1.0);
int excln = cmdline.int_val ("-excln",-1);
double excld = cmdline.double_val("-excld",-1.0);
double excly = cmdline.double_val("-excly",-1.0);
ee_print = cmdline.present("-ee-print");
bool get_all_dij = cmdline.present("-get-all-dij");
bool get_all_yij = cmdline.present("-get-all-yij");
double subdcut = cmdline.double_val("-subdcut",-1.0);
double rapmax = cmdline.double_val("-rapmax",1.0e305);
if (cmdline.present("-etamax")) {
cerr << "WARNING: -etamax options actually sets maximum rapidity (and overrides -rapmax)\n";
rapmax = cmdline.double_val("-etamax");
}
bool show_constituents = cmdline.present("-const");
bool massless = cmdline.present("-massless");
int nev = cmdline.int_val("-nev",1);
int skip = cmdline.int_val("-skip",0);
bool add_dense_coverage = cmdline.present("-dense");
double ghost_maxrap = cmdline.value("-ghost-maxrap",5.0);
bool all_algs = cmdline.present("-all-algs");
// have the option of comparing the clustering results to those
// obtained with a different clustering strategy; the misuse the
// "plugin_strategy" to indicate that no comparison is needed.
// Does not currently support areas
int compare_strategy = cmdline.value<int>("-compare-strategy", plugin_strategy);
Selector particles_sel = (cmdline.present("-nhardest"))
? SelectorNHardest(cmdline.value<unsigned int>("-nhardest"))
: SelectorIdentity();
do_areas = cmdline.present("-area");
#ifndef __FJCORE__
AreaDefinition area_def;
if (do_areas) {
assert(!write); // it's incompatible
GhostedAreaSpec ghost_spec(ghost_maxrap,
cmdline.value("-area:repeat", 1),
cmdline.value("-ghost-area", 0.01));
if (cmdline.present("-area:fj2")) ghost_spec.set_fj2_placement(true);
if (cmdline.present("-area:explicit")) {
area_def = AreaDefinition(active_area_explicit_ghosts, ghost_spec);
} else if (cmdline.present("-area:passive")) {
area_def = AreaDefinition(passive_area, ghost_spec);
} else if (cmdline.present("-area:voronoi")) {
double Rfact = cmdline.value<double>("-area:voronoi");
area_def = AreaDefinition(voronoi_area,
VoronoiAreaSpec(Rfact));
} else {
cmdline.present("-area:active"); // allow, but do not require, arg
area_def = AreaDefinition(active_area, ghost_spec);
}
}
#else
do_areas=false;
#endif
bool do_bkgd = cmdline.present("-bkgd"); // background estimation
#ifndef __FJCORE__
bool do_bkgd_csab = false, do_bkgd_jetmedian = false, do_bkgd_fj2 = false;
bool do_bkgd_gridmedian = false;
bool do_bkgd_localrange = false;
bool do_subtractor = false;
double bkgd_alt_ktR = -1.0;
BackgroundRescalingYPolynomial * bkgd_rescaling = 0;
Selector bkgd_range;
if (do_bkgd) {
bkgd_range = SelectorAbsRapMax(ghost_maxrap - ktR);
if (cmdline.present("-bkgd:csab")) {do_bkgd_csab = true;}
else if (cmdline.present("-bkgd:jetmedian")) {do_bkgd_jetmedian = true;
do_bkgd_fj2 = cmdline.present("-bkgd:fj2");
do_bkgd_localrange = cmdline.present("-bkgd:localrange");
bkgd_alt_ktR = cmdline.value("-bkgd:alt-ktR", bkgd_alt_ktR);
if (do_bkgd_localrange) bkgd_range = SelectorStrip(1.5);
} else if (cmdline.present("-bkgd:gridmedian")) {
do_bkgd_gridmedian = true;
} else {
throw Error("with the -bkgd option, some particular background must be specified (csab or jetmedian)");
}
if (cmdline.present("-bkgd:rescaling")) {
bkgd_rescaling = new BackgroundRescalingYPolynomial(1.157,0,-0.0266,0,0.000048);
}
assert(do_areas || do_bkgd_gridmedian);
do_subtractor = cmdline.present("-subtractor");
if (do_subtractor) assert(do_areas);
}
#else
do_bkgd = false;
#endif
bool show_cones = cmdline.present("-cones"); // only works for siscone
// for cone algorithms
// allow -f and -overlap
double overlap_threshold = cmdline.double_val("-overlap",0.5);
overlap_threshold = cmdline.double_val("-f",overlap_threshold);
double seed_threshold = cmdline.double_val("-seed",1.0);
#ifdef __FJCORE__
show_cones = false;
#endif
// for ee algorithms, allow to specify ycut
double ycut = cmdline.double_val("-ycut",0.08);
// for printing jets to a file for reading by root
rootfile = cmdline.value<string>("-root","");
// out default scheme is the E_scheme
RecombinationScheme scheme = E_scheme;
// The following option causes the Cambridge algo to be used.
// Note that currently the only output that works sensibly here is
// "-incl 0"
vector<JetDefinition> jet_defs;
if (all_algs || cmdline.present("-cam") || cmdline.present("-CA")) {
jet_defs.push_back( JetDefinition(cambridge_algorithm, ktR, scheme, strategy));
}
if (all_algs || cmdline.present("-antikt")) {
jet_defs.push_back( JetDefinition(antikt_algorithm, ktR, scheme, strategy));
}
if (all_algs || cmdline.present("-genkt")) {
double p;
if (cmdline.present("-genkt")) p = cmdline.value<double>("-genkt");
else p = -0.5;
jet_defs.push_back( JetDefinition(genkt_algorithm, ktR, p, scheme, strategy));
}
if (all_algs || cmdline.present("-eekt")) {
jet_defs.push_back( JetDefinition(ee_kt_algorithm));
}
if (all_algs || cmdline.present("-eegenkt")) {
double p;
if (cmdline.present("-eegenkt")) p = cmdline.value<double>("-eegenkt");
else p = -0.5;
jet_defs.push_back( JetDefinition(ee_genkt_algorithm, ktR, p, scheme, strategy));
// checking if one asks to run a plugin (don't delete this line)
}
if (all_algs || cmdline.present("-midpoint")) {
#ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
typedef CDFMidPointPlugin MPPlug; // for brevity
double cone_area_fraction = 1.0;
int max_pair_size = 2;
int max_iterations = 100;
MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
if (cmdline.present("-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
if (cmdline.present("-sm-pt")) sm_scale = MPPlug::SM_pt; // default
if (cmdline.present("-sm-mt")) sm_scale = MPPlug::SM_mt;
if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et;
jet_defs.push_back( JetDefinition( new CDFMidPointPlugin (
seed_threshold, ktR,
cone_area_fraction, max_pair_size,
max_iterations, overlap_threshold,
sm_scale)));
#else // FASTJET_ENABLE_PLUGIN_CDFCONES
is_unavailable("MidPoint");
#endif // FASTJET_ENABLE_PLUGIN_CDFCONES
}
if (all_algs || cmdline.present("-pxcone")) {
#ifdef FASTJET_ENABLE_PLUGIN_PXCONE
double min_jet_energy = 5.0;
// mode: 1=e+e-, 2=pp
int mode = cmdline.value("-pxcone-mode", 2);
bool E_scheme_jets = false;
jet_defs.push_back( JetDefinition( new PxConePlugin (
ktR, min_jet_energy,
overlap_threshold, E_scheme_jets, mode)));
#else // FASTJET_ENABLE_PLUGIN_PXCONE
is_unavailable("PxCone");
#endif // FASTJET_ENABLE_PLUGIN_PXCONE
}
if (all_algs || cmdline.present("-jetclu")) {
#ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
jet_defs.push_back( JetDefinition( new CDFJetCluPlugin (
ktR, overlap_threshold, seed_threshold)));
#else // FASTJET_ENABLE_PLUGIN_CDFCONES
is_unavailable("JetClu");
#endif // FASTJET_ENABLE_PLUGIN_CDFCONES
}
if (all_algs || cmdline.present("-siscone") || cmdline.present("-sisconespheri")) {
#ifdef FASTJET_ENABLE_PLUGIN_SISCONE
typedef SISConePlugin SISPlug; // for brevity
int npass = cmdline.value("-npass",0);
if (all_algs || cmdline.present("-siscone")) {
double sisptmin = cmdline.value("-sisptmin",0.0);
bool cache = cmdline.present("-cache");
SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass,sisptmin,cache);
if (cmdline.present("-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
if (cmdline.present("-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
if (cmdline.present("-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
if (cmdline.present("-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
// cause it to use the jet-definition's own recombiner
plugin->set_use_jet_def_recombiner(true);
jet_defs.push_back( JetDefinition(plugin));
}
if (all_algs || cmdline.present("-sisconespheri")) {
double sisEmin = cmdline.value("-sisEmin",0.0);
SISConeSphericalPlugin * plugin =
new SISConeSphericalPlugin(ktR, overlap_threshold,npass,sisEmin);
if (cmdline.present("-ghost-sep")) {
plugin->set_ghost_separation_scale(cmdline.value<double>("-ghost-sep"));
}
jet_defs.push_back( JetDefinition(plugin));
}
#else // FASTJET_ENABLE_PLUGIN_SISCONE
is_unavailable("SISCone");
#endif // FASTJET_ENABLE_PLUGIN_SISCONE
}
if (all_algs || cmdline.present("-d0runiicone")) {
#ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
double min_jet_Et = 6.0; // was 8 GeV in earlier work
jet_defs.push_back( JetDefinition(new D0RunIIConePlugin(ktR,min_jet_Et)));
#else // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
is_unavailable("D0RunIICone");
#endif // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
}
if (all_algs || cmdline.present("-trackjet")) {
#ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
jet_defs.push_back( JetDefinition(new TrackJetPlugin(ktR)));
#else // FASTJET_ENABLE_PLUGIN_TRACKJET
is_unavailable("TrackJet");
#endif // FASTJET_ENABLE_PLUGIN_TRACKJET
}
if (all_algs || cmdline.present("-atlascone")) {
#ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
jet_defs.push_back( JetDefinition(new ATLASConePlugin(ktR)));
#else // FASTJET_ENABLE_PLUGIN_ATLASCONE
is_unavailable("ATLASCone");
#endif // FASTJET_ENABLE_PLUGIN_ATLASCONE
}
if (all_algs || cmdline.present("-eecambridge")) {
#ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
jet_defs.push_back( JetDefinition(new EECambridgePlugin(ycut)));
#else // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
is_unavailable("EECambridge");
#endif // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
}
if (all_algs || cmdline.present("-jade")) {
#ifdef FASTJET_ENABLE_PLUGIN_JADE
JadePlugin::Strategy jade_strategy =
JadePlugin::Strategy(cmdline.value<int>("-jade-strategy",
JadePlugin::strategy_NNFJN2Plain));
jet_defs.push_back( JetDefinition(new JadePlugin(jade_strategy)));
#else // FASTJET_ENABLE_PLUGIN_JADE
is_unavailable("Jade");
#endif // FASTJET_ENABLE_PLUGIN_JADE
}
if (all_algs || cmdline.present("-cmsiterativecone")) {
#ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
jet_defs.push_back( JetDefinition(new CMSIterativeConePlugin(ktR,seed_threshold)));
#else // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
is_unavailable("CMSIterativeCone");
#endif // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
}
if (all_algs || cmdline.present("-d0runipre96cone")) {
#ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
jet_defs.push_back( JetDefinition(new D0RunIpre96ConePlugin(ktR, seed_threshold, overlap_threshold)));
#else // FASTJET_ENABLE_PLUGIN_D0RUNICONE
is_unavailable("D0RunIpre96Cone");
#endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
}
if (all_algs || cmdline.present("-d0runicone")) {
#ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
jet_defs.push_back( JetDefinition(new D0RunIConePlugin(ktR, seed_threshold, overlap_threshold)));
#else // FASTJET_ENABLE_PLUGIN_D0RUNICONE
is_unavailable("D0RunICone");
#endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
}
if (all_algs || cmdline.present("-gridjet")) {
#ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
// we want a grid_ymax of 5.0, but when using R=0.4 (i.e. grid
// spacing of 0.8), this leads to 12.5 grid cells; depending on
// whether this is 12.499999999999 or 12.5000000....1 this gets
// converted either to 12 or 13, making the results sensitive to
// rounding errors.
//
// Instead we therefore take 4.9999999999, which avoids this problem.
double grid_ymax = 4.9999999999;
jet_defs.push_back( JetDefinition(new GridJetPlugin(grid_ymax, ktR*2.0)));
#else // FASTJET_ENABLE_PLUGIN_GRIDJET
is_unavailable("GridJet");
#endif // FASTJET_ENABLE_PLUGIN_GRIDJET
// end of checking if one asks to run a plugin (don't delete this line)
}
if (all_algs ||
cmdline.present("-kt") ||
(jet_defs.size() == 0 && !found_unavailable)) {
jet_defs.push_back( JetDefinition(kt_algorithm, ktR, E_scheme, strategy));
}
string filename = cmdline.value<string>("-file", "");
if (!cmdline.all_options_used()) {cerr <<
"Error: some options were not recognized"<<endl;
exit(-1);}
for (unsigned idef = 0; idef < jet_defs.size(); idef++) {
JetDefinition & jet_def = jet_defs[idef];
istream * istr;
if (filename == "") istr = &cin;
else istr = new ifstream(filename.c_str());
for (int iev = 0; iev < nev; iev++) {
vector<PseudoJet> jets;
vector<PseudoJet> particles;
string line;
int ndone = 0;
while (getline(*istr, line)) {
//cout << line<<endl;
istringstream linestream(line);
if (line == "#END") {
ndone += 1;
if (ndone == combine) {break;}
}
if (line.substr(0,1) == "#") {continue;}
valarray<double> fourvec(4);
if (hydjet) {
// special reading from hydjet.txt event record (though actually
// this is supposed to be a standard pythia event record, so
// being able to read from it is perhaps not so bad an idea...)
int ii, istat,id,m1,m2,d1,d2;
double mass;
linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
>> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
// current file contains mass of particle as 4th entry
if (istat == 1) {
fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
+pow2(fourvec[2])+pow2(mass));
}
} else {
if (massless) {
linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
else {
linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
}
}
PseudoJet psjet(fourvec);
if (abs(psjet.rap()) < rapmax) {particles.push_back(psjet);}
}
#ifndef __FJCORE__
// add a fake underlying event which is very soft, uniformly distributed
// in eta,phi so as to allow one to reconstruct the area that is associated
// with each jet.
if (add_dense_coverage) {
GhostedAreaSpec ghosted_area_spec(ghost_maxrap);
//GhostedAreaSpec ghosted_area_spec(-2.0,4.0); // asymmetric range
// for plots, reduce the scatter default of 1, to avoid "holes"
// in the subsequent calorimeter view
ghosted_area_spec.set_grid_scatter(0.5);
ghosted_area_spec.add_ghosts(particles);
//----- old code ------------------
// srand(2);
// int nphi = 60;
// int neta = 100;
// double kt = 1e-1;
// for (int iphi = 0; iphi<nphi; iphi++) {
// for (int ieta = -neta; ieta<neta+1; ieta++) {
// double phi = (iphi+0.5) * (twopi/nphi) + rand()*0.001/RAND_MAX;
// double eta = ieta * (10.0/neta) + rand()*0.001/RAND_MAX;
// kt = 1e-20*(1+rand()*0.1/RAND_MAX);
// double pminus = kt*exp(-eta);
// double pplus = kt*exp(+eta);
// double px = kt*sin(phi);
// double py = kt*cos(phi);
// //cout << kt<<" "<<eta<<" "<<phi<<"\n";
// PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
// particles.push_back(mom);
// }
// }
}
#else
add_dense_coverage = false;
#endif
// select the particles that pass the selection cut
particles = particles_sel(particles);
// allow user to skip some number of events (e.g. for easier bug-chasing)
if (iev < skip) continue;
for (int irepeat = 0; irepeat < repeat ; irepeat++) {
int nparticles = particles.size();
try {
// one could use a unique_ptr here, but SharedPtr is available independently of C++ standard
SharedPtr<ClusterSequence> clust_seq;
if (do_areas) {
#ifndef __FJCORE__
clust_seq.reset(new ClusterSequenceArea(particles,jet_def,area_def));
#endif
} else {
clust_seq.reset(new ClusterSequence(particles,jet_def,write));
}
if (compare_strategy != plugin_strategy) {
do_compare_strategy(iev, particles, jet_def, *clust_seq, compare_strategy);
}
// repetitive output
if (repeatinclkt >= 0.0) {
vector<PseudoJet> jets_local = sorted_by_pt(clust_seq->inclusive_jets(repeatinclkt));
}
if (irepeat != 0) {continue;}
cout << "iev "<<iev<< ": number of particles = "<< nparticles << endl;
cout << "strategy used = "<< clust_seq->strategy_string()<< endl;
if (iev == 0) cout << "Jet Definition: " << jet_def.description() << " (" << fastjet_version_string() << ")" << endl;
#ifndef __FJCORE__
if (do_areas && iev == 0) cout << "Area definition: " << area_def.description() << endl;
#endif
// now provide some nice output...
if (inclkt >= 0.0) {
vector<PseudoJet> jets_local = sorted_by_pt(clust_seq->inclusive_jets(inclkt));
print_jets(jets_local, show_constituents);
}
if (excln > 0) {
cout << "Printing "<<excln<<" exclusive jets\n";
print_jets(clust_seq->exclusive_jets(excln), show_constituents);
}
if (excld > 0.0) {
cout << "Printing exclusive jets for d = "<<excld<<"\n";
print_jets(clust_seq->exclusive_jets(excld), show_constituents);
}
if (excly > 0.0) {
cout << "Printing exclusive jets for ycut = "<<excly<<"\n";
print_jets(clust_seq->exclusive_jets_ycut(excly), show_constituents);
}
if (get_all_dij) {
for (int i = nparticles-1; i >= 0; i--) {
printf("d for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_dmerge(i));
}
}
if (get_all_yij) {
for (int i = nparticles-1; i >= 0; i--) {
printf("y for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_ymerge(i));
}
}
// have the option of printing out the subjets (at scale dcut) of
// each inclusive jet
if (subdcut >= 0.0) {
print_jets_and_sub(clust_seq->inclusive_jets(), subdcut);
}
// useful for testing that recombination sequences are unique
if (unique_write) {
vector<int> unique_history = clust_seq->unique_history_order();
// construct the inverse of the above mapping
vector<int> inv_unique_history(clust_seq->history().size());
for (unsigned int i = 0; i < unique_history.size(); i++) {
inv_unique_history[unique_history[i]] = i;}
for (unsigned int i = 0; i < unique_history.size(); i++) {
ClusterSequence::history_element el =
clust_seq->history()[unique_history[i]];
int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
}
}
#ifdef FASTJET_ENABLE_PLUGIN_SISCONE
// provide some complementary information for SISCone
if (show_cones) {
const SISConeExtras * extras =
dynamic_cast<const SISConeExtras *>(clust_seq->extras());
if (extras == 0)
throw fastjet::Error("extras object for SISCone was null (this can happen with certain area types)");
cout << "most ambiguous split (difference in squared dist) = "
<< extras->most_ambiguous_split() << endl;
vector<fastjet::PseudoJet> stable_cones(extras->stable_cones());
stable_cones = sorted_by_rapidity(stable_cones);
for (unsigned int i = 0; i < stable_cones.size(); i++) {
//if (stable_cones[i].phi() < 5.0 && stable_cones[i].phi() > 4.0) {
printf("%5u %15.8f %15.8f %15.8f\n",
i,stable_cones[i].rap(),stable_cones[i].phi(),
stable_cones[i].perp() );
//}
}
// also show passes for jets
vector<PseudoJet> sisjets = clust_seq->inclusive_jets();
printf("\n%15s %15s %15s %12s %8s %8s\n","rap","phi","pt","user-index","pass","nconst");
for (unsigned i = 0; i < sisjets.size(); i++) {
printf("%15.8f %15.8f %15.8f %12d %8d %8u\n",
sisjets[i].rap(), sisjets[i].phi(), sisjets[i].perp(),
sisjets[i].user_index(), extras->pass(sisjets[i]),
(unsigned int) clust_seq->constituents(sisjets[i]).size()
);
}
}
#endif // FASTJET_ENABLE_PLUGIN_SISCONE
#ifndef __FJCORE__
if (do_bkgd) {
double rho, sigma, mean_area, empty_area, n_empty_jets;
ClusterSequenceAreaBase * csab =
dynamic_cast<ClusterSequenceAreaBase *>(clust_seq.get());
BackgroundEstimatorBase * bge_ptr = 0;
if (do_bkgd_csab) {
csab->get_median_rho_and_sigma(bkgd_range, true, rho, sigma, mean_area);
empty_area = csab->empty_area(bkgd_range);
n_empty_jets = csab->n_empty_jets(bkgd_range);
} else if (do_bkgd_jetmedian) {
JetMedianBackgroundEstimator * bge = new JetMedianBackgroundEstimator(bkgd_range);
bge_ptr = bge;
// may be null
bge->set_rescaling_class(bkgd_rescaling);
bge->set_provide_fj2_sigma(do_bkgd_fj2);
if (bkgd_alt_ktR > 0) {
ClusterSequenceAreaBase * clust_seq_bkgd =
new ClusterSequenceArea(particles, JetDefinition(kt_algorithm, bkgd_alt_ktR), area_def);
cout << "Alt JetDef for background-estimation CSAB: " << clust_seq_bkgd->jet_def().description() << endl;
bge->set_cluster_sequence(*clust_seq_bkgd);
clust_seq_bkgd->delete_self_when_unused();
} else {
bge->set_cluster_sequence(*csab);
}
if (!do_bkgd_localrange) {
rho = bge->rho();
sigma = bge->sigma();
mean_area = bge->mean_area();
empty_area = bge->empty_area();
n_empty_jets = bge->n_empty_jets();
}
} else {
assert(do_bkgd_gridmedian);
double grid_rapmin, grid_rapmax;
bkgd_range.get_rapidity_extent(grid_rapmin, grid_rapmax);
GridMedianBackgroundEstimator * bge = new GridMedianBackgroundEstimator(grid_rapmax, 2*ktR);
bge_ptr = bge;
bge->set_rescaling_class(bkgd_rescaling);
bge->set_particles(particles);
rho = bge->rho();
sigma = bge->sigma();
mean_area = bge->mean_area();
empty_area = 0;
n_empty_jets = 0;
}
if (do_bkgd_localrange || do_subtractor) {
assert(bge_ptr != 0);
cout << "Background estimator: " << bge_ptr->description() << endl;
//vector<PseudoJet>
jets = SelectorAbsRapMax(3.0)(sorted_by_pt(csab->inclusive_jets()));
vector<PseudoJet> subjets;
if (do_subtractor) {
Subtractor subtractor(bge_ptr);
subtractor.set_use_rho_m(true);
subtractor.set_safe_mass(true);
cout << "Subtractor: " << subtractor.description() << endl;
subjets = subtractor(jets);
}
print_jets_bkgd(jets, subjets, bge_ptr, do_subtractor);
// cout << "i pt rap phi m rho rho_m sigma sigma_m" << endl;
// if (do_subtractor) cout << "isub ptsub rapsub phisub msub area" << endl;
// for (unsigned i = 0; i < jets.size(); i++) {
// const PseudoJet & jet = jets[i];
// cout << i << " "
// << " " << jet.pt() << " " << jet.rap() << " " << jet.phi() << " " << jet.m()
// << " " << bge_ptr->rho(jet) << " " << bge_ptr->rho_m(jet)
// << " " << bge_ptr->sigma(jet) << " " << bge_ptr->sigma_m(jet) << endl;
// if (do_subtractor) {
// const PseudoJet & subjet = subjets[i];
// cout << i << "sub"
// << " " << subjet.pt() << " " << subjet.rap() << " " << subjet.phi() << " " << subjet.m()
// << " " << jet.area() << endl;
// }
// }
} else {
cout << " rho = " << rho
<< ", sigma = " << sigma
<< ", mean_area = " << mean_area
<< ", empty_area = " << empty_area
<< ", n_empty_jets = " << n_empty_jets
<< endl;
}
if (bge_ptr != 0) delete bge_ptr;
}
#endif
} // try
catch (Error &fjerr) {
cout << "Caught fastjet error, exiting gracefully" << endl;
exit(0);
}
} // irepeat
} // iev
// if we've instantiated a plugin, delete it
if (jet_def.strategy()==plugin_strategy){
delete jet_def.plugin();
}
// close any file that we've opened
if (istr != &cin) delete istr;
} //
}
//------ HELPER ROUTINES -----------------------------------------------
/// print a single jet
void print_jet (const PseudoJet & jet) {
unsigned int n_constituents = jet.constituents().size();
printf("%15.8f %15.8f %15.8f %8u\n",
jet.rap(), jet.phi(), jet.perp(), n_constituents);
}
//----------------------------------------------------------------------
void print_jets(const vector<PseudoJet> & jets_in, bool show_constituents) {
vector<PseudoJet> jets;
if (ee_print) {
jets = sorted_by_E(jets_in);
for (unsigned int j = 0; j < jets.size(); j++) {
printf("%5u %15.8f %15.8f %15.8f %15.8f\n",
j,jets[j].px(),jets[j].py(),jets[j].pz(),jets[j].E());
if (show_constituents) {
vector<PseudoJet> const_jets = jets[j].constituents();
for (unsigned int k = 0; k < const_jets.size(); k++) {
printf(" jet%03u %15.8f %15.8f %15.8f %15.8f\n",j,const_jets[k].px(),
const_jets[k].py(),const_jets[k].pz(),const_jets[k].E());
}
cout << "\n\n";
}
}
} else {
jets = sorted_by_pt(jets_in);
for (unsigned int j = 0; j < jets.size(); j++) {
printf("%5u %15.8f %15.8f %15.8f",
j,jets[j].rap(),jets[j].phi(),jets[j].perp());
// also print out the scalar area and the perp component of the
// 4-vector (just enough to check a reasonable 4-vector?)
#ifndef __FJCORE__
if (do_areas) printf(" %15.8f %15.8f", jets[j].area(),
jets[j].area_4vector().perp());
cout << "\n";
#endif
if (show_constituents) {
vector<PseudoJet> const_jets = jets[j].constituents();
for (unsigned int k = 0; k < const_jets.size(); k++) {
printf(" jet%03u %15.8f %15.8f %15.8f %5d\n",j,const_jets[k].rap(),
const_jets[k].phi(),sqrt(const_jets[k].kt2()), const_jets[k].cluster_hist_index());
}
cout << "\n\n";
}
}
}
if (rootfile != "") {
ofstream ostr(rootfile.c_str());
ostr << "# " << cmdline_p->command_line() << endl;
ostr << "# output for root" << endl;
assert(jets.size() > 0);
jets[0].validated_cs()->print_jets_for_root(jets,ostr);
}
}
//----- SUBJETS --------------------------------------------------------
/// a function that pretty prints a list of jets and the subjets for each
/// one
void print_jets_and_sub (const vector<PseudoJet> & jets, double dcut) {
// sort jets into increasing pt
vector<PseudoJet> sorted_jets = sorted_by_pt(jets);
// label the columns
printf("Printing jets and their subjets with subdcut = %10.5f\n",dcut);
printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity",
"phi", "pt", "n constituents");
// the kind of subjet finding used to test consistency among them
SubType sub_type = subtype_internal;
//SubType sub_type = subtype_newclust_dcut;
//SubType sub_type = subtype_newclust_R;
// print out the details for each jet
//for (unsigned int i = 0; i < sorted_jets.size(); i++) {
for (vector<PseudoJet>::const_iterator jet = sorted_jets.begin();
jet != sorted_jets.end(); jet++) {
const JetDefinition & jet_def = jet->validated_cs()->jet_def();
// if jet pt^2 < dcut with kt alg, then some methods of
// getting subjets will return nothing -- so skip the jet
if (jet_def.jet_algorithm() == kt_algorithm
&& jet->perp2() < dcut) continue;
printf("%5u ",(unsigned int) (jet - sorted_jets.begin()));
print_jet(*jet);
vector<PseudoJet> subjets;
ClusterSequence * cspoint;
if (sub_type == subtype_internal) {
cspoint = 0;
subjets = jet->exclusive_subjets(dcut);
double ddnp1 = jet->exclusive_subdmerge_max(subjets.size());
double ddn = jet->exclusive_subdmerge_max(subjets.size()-1);
cout << " for " << ddnp1 << " < d < " << ddn << " one has " << endl;
} else if (sub_type == subtype_newclust_dcut) {
cspoint = new ClusterSequence(jet->constituents(), jet_def);
subjets = cspoint->exclusive_jets(dcut);
} else if (sub_type == subtype_newclust_R) {
assert(jet_def.jet_algorithm() == cambridge_algorithm);
JetDefinition subjd(jet_def.jet_algorithm(), jet_def.R()*sqrt(dcut));
cspoint = new ClusterSequence(jet->constituents(), subjd);
subjets = cspoint->inclusive_jets();
} else {
cerr << "unrecognized subtype for subjet finding" << endl;
exit(-1);
}
subjets = sorted_by_pt(subjets);
for (unsigned int j = 0; j < subjets.size(); j++) {
printf(" -sub-%02u ",j);
print_jet(subjets[j]);
}
if (cspoint != 0) delete cspoint;
//ClusterSequence subseq(clust_seq->constituents(sorted_jets[i]),
// JetDefinition(cambridge_algorithm, 0.4));
//vector<PseudoJet> subjets = sorted_by_pt(subseq.inclusive_jets());
//for (unsigned int j = 0; j < subjets.size(); j++) {
// printf(" -sub-%02u ",j);
// print_jet(subseq, subjets[j]);
//}
}
}
/// if abs(x)<precision/2, return 0
double make_safe_zero_truncation(double x, double precision){
return std::abs(x)<0.5*precision ? 0.0 : x;
}
#ifndef __FJCORE__
void print_jets_bkgd(const vector<PseudoJet> &jets,
const vector<PseudoJet> &subtracted_jets,
BackgroundEstimatorBase * bge_ptr,
bool do_subtractor){
printf("Printing jets, background information");
if (do_subtractor)
printf(" and subtracted jets\n");
printf("%5s %15s %15s %15s %15s %15s %15s %15s %15s\n","jet #",
"rapidity", "phi", "pt", "pt^2+m^2",
"rho", "rho_m", "sigma", "sigma_m");
if (do_subtractor)
printf("%5s %15s %15s %15s %15s %15ss\n","jet #",
"rapidity", "phi", "pt", "sqrt(pt^2+m^2)", "area");
for (unsigned i = 0; i < jets.size(); i++) {
const PseudoJet & jet = jets[i];
BackgroundEstimate estimate = bge_ptr->estimate(jet);
// Note that the values of rho_m sometimes comes out as +- a very
// small number and the format can produce either 0.00000000 or
// -0.00000000. The call to "make_safe_zero_truncation" makes sure it is
// printed wo the - sign in each case
printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n", i,
jet.rap(), jet.phi(), jet.perp(), jet.mt(),
estimate.rho(), make_safe_zero_truncation(estimate.rho_m(),1e-8),
estimate.sigma(), estimate.sigma_m());
if (do_subtractor) {
const PseudoJet & subjet = subtracted_jets[i];
printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f\n", i,
subjet.rap(), subjet.phi(), subjet.perp(), subjet.mt(), jet.area());
}
}
}
#endif// __FJCORE__
//----------------------------------------------------------------------
void signal_failed_comparison(int iev,
const string & message,
const vector<PseudoJet> & particles) {
cout << "# failed comparison, reason is " << message << endl;
cout << "# iev = " << iev << endl;
cout << setprecision(16);
for (unsigned i = 0; i < particles.size(); i++) {
const PseudoJet & p = particles[i];
cout << p.px() << " "
<< p.py() << " "
<< p.pz() << " "
<< p.E () << endl;
}
cout << "#END" << endl;
}
//----------------------------------------------------------------------
void do_compare_strategy(int iev,
const vector<PseudoJet> & particles,
const JetDefinition & jet_def,
const ClusterSequence & cs,
int compare_strategy) {
// create a jet def with the reference comparison strategy
JetDefinition jet_def_ref(jet_def.jet_algorithm(),
jet_def.R(),
jet_def.recombination_scheme(),
Strategy(compare_strategy));
// do the clustering
ClusterSequence cs_ref(particles, jet_def_ref);
// now compare the outputs. At this stage just based on the clustering
// sequence - get more sophisticated later...
const vector<ClusterSequence::history_element> & history_in = cs.history();
const vector<ClusterSequence::history_element> & history_ref = cs_ref.history();
if (history_in.size() != history_ref.size()) {
signal_failed_comparison(iev, "history sizes do not match", particles);
}
// now run over each clustering step to do the comparison
for (unsigned i = cs.n_particles(); i < history_in.size(); i++) {
bool fail_parents = (history_in[i].parent1 != history_ref[i].parent1 ||
history_in[i].parent2 != history_ref[i].parent2);
bool fail_dij = (history_in[i].dij != history_ref[i].dij);
if (fail_parents || fail_dij) {
ostringstream ostr;
ostr << "at step " << i << ", ";
if (fail_parents) ostr << "parents do not match, ";
if (fail_dij) ostr << "dij does not match, ";
ostr << "history in (p1, p2, dij) = "
<< history_in[i].parent1 << " " << history_in[i].parent2 << " " << history_in[i].dij;
ostr << ", history ref (p1, p2, dij) = "
<< history_ref[i].parent1 << " " << history_ref[i].parent2 << " " << history_ref[i].dij;
signal_failed_comparison(iev, ostr.str(), particles);
break;
}
}
}
main
int main()
an example program showing how to use fastjet
Definition 01-basic.cc:50
CmdLine
Definition CmdLine.hh:53
fastjet::ATLASConePlugin
Implementation of the ATLAS Cone (plugin for fastjet v2.4 upwards).
Definition ATLASConePlugin.hh:52
fastjet::AreaDefinition
class that holds a generic area definition
Definition AreaDefinition.hh:82
fastjet::BackgroundEstimate
/// a class that holds the result of the calculation
Definition BackgroundEstimatorBase.hh:53
fastjet::BackgroundEstimate::sigma_m
double sigma_m() const
fluctuations in the purely longitudinal (particle-mass-induced) component of the background density p...
Definition BackgroundEstimatorBase.hh:82
fastjet::BackgroundEstimate::rho
double rho() const
background density per unit area
Definition BackgroundEstimatorBase.hh:66
fastjet::BackgroundEstimate::sigma
double sigma() const
background fluctuations per unit square-root area must be multipled by sqrt(area) to get fluctuations...
Definition BackgroundEstimatorBase.hh:71
fastjet::BackgroundEstimate::rho_m
double rho_m() const
purely longitudinal (particle-mass-induced) component of the background density per unit area
Definition BackgroundEstimatorBase.hh:78
fastjet::BackgroundEstimatorBase
Abstract base class that provides the basic interface for classes that estimate levels of background ...
Definition BackgroundEstimatorBase.hh:183
fastjet::BackgroundEstimatorBase::description
virtual std::string description() const =0
returns a textual description of the background estimator
fastjet::BackgroundEstimatorBase::estimate
virtual BackgroundEstimate estimate() const =0
get the full set of background properties
fastjet::BackgroundRescalingYPolynomial
A background rescaling that is a simple polynomial in y.
Definition BackgroundEstimatorBase.hh:412
fastjet::CMSIterativeConePlugin
Implementation of the CMS Iterative Cone (plugin for fastjet v2.4 upwards).
Definition CMSIterativeConePlugin.hh:49
fastjet::ClusterSequenceAreaBase
base class that sets interface for extensions of ClusterSequence that provide information about the a...
Definition ClusterSequenceAreaBase.hh:48
fastjet::ClusterSequenceAreaBase::get_median_rho_and_sigma
virtual void get_median_rho_and_sigma(const Selector &selector, bool use_area_4vector, double &median, double &sigma, double &mean_area) const
using jets withing the selector range (and with 4-vector areas if use_area_4vector),...
Definition ClusterSequenceAreaBase.cc:200
fastjet::ClusterSequenceAreaBase::n_empty_jets
virtual double n_empty_jets(const Selector &selector) const
return something similar to the number of pure ghost jets in the given selector's range in an active ...
Definition ClusterSequenceAreaBase.hh:136
fastjet::ClusterSequenceAreaBase::empty_area
virtual double empty_area(const Selector &selector) const
return the total area, corresponding to the given Selector, that is free of jets, in general based on...
Definition ClusterSequenceAreaBase.cc:57
fastjet::ClusterSequenceArea
General class for user to obtain ClusterSequence with additional area information.
Definition ClusterSequenceArea.hh:51
fastjet::ClusterSequence
deals with clustering
Definition ClusterSequence.hh:63
fastjet::ClusterSequence::jet_def
const JetDefinition & jet_def() const
return a reference to the jet definition
Definition ClusterSequence.hh:283
fastjet::ClusterSequence::history
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
Definition ClusterSequence.hh:1007
fastjet::ClusterSequence::delete_self_when_unused
void delete_self_when_unused()
by calling this routine you tell the ClusterSequence to delete itself when all the Pseudojets associa...
Definition ClusterSequence.cc:1778
fastjet::ClusterSequence::n_particles
unsigned int n_particles() const
returns the number of particles that were provided to the clustering algorithm (helps the user find t...
Definition ClusterSequence.hh:1011
fastjet::ClusterSequence::inclusive_jets
std::vector< PseudoJet > inclusive_jets(const double ptmin=0.0) const
return a vector of all jets (in the sense of the inclusive algorithm) with pt >= ptmin.
Definition ClusterSequence.cc:916
fastjet::ClusterSequence::print_banner
static void print_banner()
This is the function that is automatically called during clustering to print the FastJet banner.
Definition ClusterSequence.cc:464
fastjet::ClusterSequence::exclusive_jets
std::vector< PseudoJet > exclusive_jets(const double dcut) const
return a vector of all jets (in the sense of the exclusive algorithm) that would be obtained when run...
Definition ClusterSequence.cc:987
fastjet::D0RunIConePlugin
A plugin for FastJet (v3.0 or later) that provides an interface to the D0 version of Run-I cone algor...
Definition D0RunIConePlugin.hh:63
fastjet::D0RunIIConePlugin
Implementation of the D0 Run II Cone (plugin for fastjet v2.1 upwards).
Definition D0RunIIConePlugin.hh:61
fastjet::D0RunIpre96ConePlugin
A plugin for FastJet (v3.0 or later) that provides an interface to the pre 1996 D0 version of Run-I c...
Definition D0RunIpre96ConePlugin.hh:63
fastjet::EECambridgePlugin
Implementation of the e+e- Cambridge algorithm (plugin for fastjet v2.4 upwards).
Definition EECambridgePlugin.hh:59
fastjet::Error
base class corresponding to errors that can be thrown by FastJet
Definition Error.hh:52
fastjet::GhostedAreaSpec
Parameters to configure the computation of jet areas using ghosts.
Definition GhostedAreaSpec.hh:103
fastjet::GridJetPlugin
plugin for fastjet (v3.0 upwards) that clusters particles such that all particles in a given cell of ...
Definition GridJetPlugin.hh:65
fastjet::GridMedianBackgroundEstimator
Background Estimator based on the median pt/area of a set of grid cells.
Definition GridMedianBackgroundEstimator.hh:69
fastjet::GridMedianBackgroundEstimator::sigma
double sigma() const override
returns sigma, the background fluctuations per unit area; must be multipled by sqrt(area) to get fluc...
Definition GridMedianBackgroundEstimator.cc:168
fastjet::GridMedianBackgroundEstimator::set_rescaling_class
virtual void set_rescaling_class(const FunctionOfPseudoJet< double > *rescaling_class) override
Set a pointer to a class that calculates the rescaling factor as a function of the jet (position).
Definition GridMedianBackgroundEstimator.cc:268
fastjet::GridMedianBackgroundEstimator::set_particles
void set_particles(const std::vector< PseudoJet > &particles) override
tell the background estimator that it has a new event, composed of the specified particles.
Definition GridMedianBackgroundEstimator.cc:43
fastjet::GridMedianBackgroundEstimator::mean_area
double mean_area() const
returns the area of the grid cells (all identical, but referred to as "mean" area for uniformity with...
Definition GridMedianBackgroundEstimator.hh:191
fastjet::GridMedianBackgroundEstimator::rho
double rho() const override
returns rho, the median background density per unit area
Definition GridMedianBackgroundEstimator.cc:158
fastjet::JadePlugin
Implementation of the e+e- Jade algorithm (plugin for fastjet v2.4 upwards).
Definition JadePlugin.hh:77
fastjet::JadePlugin::Strategy
Strategy
enum that contains the two clustering strategy options; for higher multiplicities,...
Definition JadePlugin.hh:82
fastjet::JetDefinition
class that is intended to hold a full definition of the jet clusterer
Definition JetDefinition.hh:250
fastjet::JetDefinition::description
std::string description() const
return a textual description of the current jet definition
Definition JetDefinition.cc:106
fastjet::JetDefinition::plugin
const Plugin * plugin() const
return a pointer to the plugin
Definition JetDefinition.hh:425
fastjet::JetDefinition::jet_algorithm
JetAlgorithm jet_algorithm() const
return information about the definition...
Definition JetDefinition.hh:432
fastjet::JetMedianBackgroundEstimator
Class to estimate the pt density of the background per unit area, using the median of the distributio...
Definition JetMedianBackgroundEstimator.hh:86
fastjet::JetMedianBackgroundEstimator::set_provide_fj2_sigma
void set_provide_fj2_sigma(bool provide_fj2_sigma=true)
The FastJet v2.X sigma calculation had a small spurious offset in the limit of a small number of jets...
Definition JetMedianBackgroundEstimator.hh:351
fastjet::JetMedianBackgroundEstimator::n_empty_jets
double n_empty_jets() const
Returns the number of empty jets used when computing the background properties.
Definition JetMedianBackgroundEstimator.cc:547
fastjet::JetMedianBackgroundEstimator::empty_area
double empty_area() const
Returns the estimate of the area (within the range defined by the selector) that is not occupied by j...
Definition JetMedianBackgroundEstimator.cc:518
fastjet::JetMedianBackgroundEstimator::sigma
double sigma() const override
get sigma, the background fluctuations per unit area
Definition JetMedianBackgroundEstimator.cc:314
fastjet::JetMedianBackgroundEstimator::mean_area
double mean_area() const
Returns the mean area of the jets used to actually compute the background properties in the last call...
Definition JetMedianBackgroundEstimator.cc:437
fastjet::JetMedianBackgroundEstimator::rho
double rho() const override
get rho, the median background density per unit area
Definition JetMedianBackgroundEstimator.cc:303
fastjet::JetMedianBackgroundEstimator::set_rescaling_class
virtual void set_rescaling_class(const FunctionOfPseudoJet< double > *rescaling_class_in) override
Set a pointer to a class that calculates the rescaling factor as a function of the jet (position).
Definition JetMedianBackgroundEstimator.hh:378
fastjet::JetMedianBackgroundEstimator::set_cluster_sequence
void set_cluster_sequence(const ClusterSequenceAreaBase &csa)
(re)set the cluster sequence (with area support) to be used by future calls to rho() etc.
Definition JetMedianBackgroundEstimator.cc:183
fastjet::PseudoJet
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition PseudoJet.hh:68
fastjet::PseudoJet::constituents
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition PseudoJet.cc:692
fastjet::PseudoJet::mt
double mt() const
returns the transverse mass = sqrt(kt^2+m^2)
Definition PseudoJet.hh:175
fastjet::PseudoJet::rap
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition PseudoJet.hh:138
fastjet::PseudoJet::phi
double phi() const
returns phi (in the range 0..2pi)
Definition PseudoJet.hh:123
fastjet::PseudoJet::perp
double perp() const
returns the scalar transverse momentum
Definition PseudoJet.hh:158
fastjet::PseudoJet::perp2
double perp2() const
returns the squared transverse momentum
Definition PseudoJet.hh:156
fastjet::PseudoJet::exclusive_subdmerge_max
double exclusive_subdmerge_max(int nsub) const
Returns the maximum dij that occurred in the whole event at the stage that the nsub+1 -> nsub merge o...
Definition PseudoJet.cc:774
fastjet::PseudoJet::area
virtual double area() const
return the jet (scalar) area.
Definition PseudoJet.cc:829
fastjet::PseudoJet::exclusive_subjets
std::vector< PseudoJet > exclusive_subjets(const double dcut) const
return a vector of all subjets of the current jet (in the sense of the exclusive algorithm) that woul...
Definition PseudoJet.cc:715
fastjet::PseudoJet::validated_cs
const ClusterSequence * validated_cs() const
shorthand for validated_cluster_sequence()
Definition PseudoJet.cc:565
fastjet::SISConeBaseExtras::pass
int pass(const PseudoJet &jet) const
return the # of the pass at which a given jet was found; will return -1 if the pass is invalid
Definition SISConeBasePlugin.hh:246
fastjet::SISConeBaseExtras::most_ambiguous_split
double most_ambiguous_split() const
return the smallest difference in squared distance encountered during splitting between a particle an...
Definition SISConeBasePlugin.hh:260
fastjet::SISConeBaseExtras::stable_cones
const std::vector< PseudoJet > & stable_cones() const
returns a reference to the vector of stable cones (aka protocones)
Definition SISConeBasePlugin.hh:239
fastjet::SISConeExtras
Class that provides extra information about a SISCone clustering.
Definition SISConePlugin.hh:250
fastjet::Selector
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition Selector.hh:149
fastjet::Selector::get_rapidity_extent
void get_rapidity_extent(double &rapmin, double &rapmax) const
returns the rapidity range for which it may return "true"
Definition Selector.hh:232
fastjet::SharedPtr
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
Definition SharedPtr.hh:341
fastjet::SharedPtr::get
T * get() const
get the stored pointer
Definition SharedPtr.hh:473
fastjet::SharedPtr::reset
void reset()
reset the pointer to default value (NULL)
Definition SharedPtr.hh:381
fastjet::Subtractor
Class that helps perform jet background subtraction.
Definition Subtractor.hh:62
fastjet::TrackJetPlugin
Implementation of the TrackJet algorithm (plugin for fastjet v2.4 upwards).
Definition TrackJetPlugin.hh:49
fastjet::VoronoiAreaSpec
Specification for the computation of the Voronoi jet area.
Definition AreaDefinition.hh:50
fastjet::SelectorNHardest
Selector SelectorNHardest(unsigned int n)
select the n hardest objects
Definition Selector.cc:1074
fastjet::SelectorStrip
Selector SelectorStrip(const double half_width)
select objets within a rapidity distance 'half_width' from the location of the reference jet,...
Definition Selector.cc:1267
fastjet::SelectorAbsRapMax
Selector SelectorAbsRapMax(double absrapmax)
select objects with |rap| <= absrapmax
Definition Selector.cc:880
fastjet
the FastJet namespace
Definition AreaDefinition.cc:37
fastjet::sorted_by_rapidity
vector< PseudoJet > sorted_by_rapidity(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing rapidity
Definition PseudoJet.cc:890
fastjet::Strategy
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
Definition JetDefinition.hh:51
fastjet::Best
@ Best
automatic selection of the best (based on N), including the LazyTiled strategies that are new to FJ3....
Definition JetDefinition.hh:102
fastjet::plugin_strategy
@ plugin_strategy
the plugin has been used...
Definition JetDefinition.hh:127
fastjet::RecombinationScheme
RecombinationScheme
The various recombination schemes.
Definition JetDefinition.hh:193
fastjet::E_scheme
@ E_scheme
summing the 4-momenta
Definition JetDefinition.hh:195
fastjet::sorted_by_E
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
Definition PseudoJet.cc:898
fastjet::cambridge_algorithm
@ cambridge_algorithm
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
Definition JetDefinition.hh:142
fastjet::kt_algorithm
@ kt_algorithm
the longitudinally invariant kt algorithm
Definition JetDefinition.hh:139
fastjet::fastjet_version_string
string fastjet_version_string()
return a string containing information about the release
Definition ClusterSequence.cc:457
fastjet::sorted_by_pt
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition PseudoJet.cc:882
fastjet::ClusterSequence::history_element
a single element in the clustering history
Definition ClusterSequence.hh:416
fastjet::ClusterSequence::history_element::parent2
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
Definition ClusterSequence.hh:421
fastjet::ClusterSequence::history_element::dij
double dij
index in the _jets vector where we will find the
Definition ClusterSequence.hh:439