42#include "fastjet/PseudoJet.hh"
43#include "fastjet/ClusterSequenceArea.hh"
44#include "fastjet/Selector.hh"
45#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
46#include "fastjet/tools/Subtractor.hh"
48#include <fastjet/config.h>
61 vector<PseudoJet> hard_event, full_event;
67 double particle_maxrap = 5.0;
71 while (getline(cin, line)) {
72 istringstream linestream(line);
75 if (line.substr(0,4) ==
"#END") {
break;}
76 if (line.substr(0,9) ==
"#SUBSTART") {
78 if (nsub == 1) hard_event = full_event;
81 if (line.substr(0,1) ==
"#") {
continue;}
83 linestream >> px >> py >> pz >> E;
88 if (abs(particle.rap()) <= particle_maxrap) full_event.push_back(particle);
92 if (nsub == 1) hard_event = full_event;
96 cerr <<
"Error: read empty event\n";
112 double ghost_maxrap = 6.0;
127 vector<PseudoJet> hard_jets =
sorted_by_pt(clust_seq_hard.inclusive_jets(ptmin));
128 vector<PseudoJet> full_jets =
sorted_by_pt(clust_seq_full.inclusive_jets(ptmin));
183#if FASTJET_VERSION_NUMBER >= 30100
184 subtractor.set_use_rho_m(
true);
185 subtractor.set_safe_mass(
true);
194 bkgd_estimator.set_particles(full_event);
201 cout <<
"Main clustering:" << endl;
202 cout <<
" Ran: " << jet_def.description() << endl;
203 cout <<
" Area: " << area_def.description() << endl;
204 cout <<
" Particles up to |y|=" << particle_maxrap << endl;
207 cout <<
"Background estimation:" << endl;
208 cout <<
" " << bkgd_estimator.description() << endl << endl;;
209 cout <<
" Giving, for the full event" << endl;
211 cout <<
" rho = " << bkgd_estimate.
rho() << endl;
212 cout <<
" sigma = " << bkgd_estimate.
sigma() << endl;
213 cout <<
" rho_m = " << bkgd_estimate.
rho_m() << endl;
214 cout <<
" sigma_m = " << bkgd_estimate.
sigma_m() << endl;
217 cout <<
"Jets above " << ptmin <<
" GeV in the hard event (" << hard_event.size() <<
" particles)" << endl;
218 cout <<
"---------------------------------------\n";
219 printf(
"%5s %15s %15s %15s %15s %15s\n",
"jet #",
"rapidity",
"phi",
"pt",
"m",
"area");
220 for (
unsigned int i = 0; i < hard_jets.size(); i++) {
221 printf(
"%5u %15.8f %15.8f %15.8f %15.8f %15.8f\n", i,
222 hard_jets[i].rap(), hard_jets[i].phi(), hard_jets[i].pt(), hard_jets[i].m(),
223 hard_jets[i].area());
233 cout <<
"Jets above " << ptmin <<
" GeV in the full event (" << full_event.size() <<
" particles)" << endl;
234 cout <<
"---------------------------------------\n";
235 printf(
"%5s %15s %15s %15s %15s %15s %15s %15s %15s %15s\n",
"jet #",
"rapidity",
"phi",
"pt",
"m",
"area",
"rap_sub",
"phi_sub",
"pt_sub",
"m_sub");
239 vector<PseudoJet> subtracted_jets = subtractor(full_jets);
241 for (
unsigned int i=0; i<full_jets.size(); i++){
243 if (subtracted_jets[i].pt2() >= ptmin*ptmin){
244 printf(
"%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n", idx,
245 full_jets[i].rap(), full_jets[i].phi(), full_jets[i].pt(), full_jets[i].m(),
247 subtracted_jets[i].rap(), subtracted_jets[i].phi(),
248 subtracted_jets[i].pt(),
249 subtracted_jets[i].m());
int main()
an example program showing how to use fastjet
class that holds a generic area definition
/// a class that holds the result of the calculation
double sigma_m() const
fluctuations in the purely longitudinal (particle-mass-induced) component of the background density p...
double rho() const
background density per unit area
double sigma() const
background fluctuations per unit square-root area must be multipled by sqrt(area) to get fluctuations...
double rho_m() const
purely longitudinal (particle-mass-induced) component of the background density per unit area
General class for user to obtain ClusterSequence with additional area information.
Parameters to configure the computation of jet areas using ghosts.
class that is intended to hold a full definition of the jet clusterer
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Class that helps perform jet background subtraction.
Selector SelectorNHardest(unsigned int n)
select the n hardest objects
Selector SelectorAbsRapMax(double absrapmax)
select objects with |rap| <= absrapmax
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2