|
fastjet 2.4.5
|
|
fastjet 2.4.5
|
#include "fastjet/PseudoJet.hh"#include "fastjet/ClusterSequenceArea.hh"#include <iostream>#include <sstream>#include <fstream>#include <valarray>#include <vector>#include <cstdlib>#include <cstddef>#include "CmdLine.hh"#include "fastjet/config.h"#include "fastjet/SISConePlugin.hh"#include "fastjet/SISConeSphericalPlugin.hh"#include "fastjet/CDFMidPointPlugin.hh"#include "fastjet/CDFJetCluPlugin.hh"#include "fastjet/D0RunIIConePlugin.hh"#include "fastjet/TrackJetPlugin.hh"#include "fastjet/ATLASConePlugin.hh"#include "fastjet/EECambridgePlugin.hh"#include "fastjet/JadePlugin.hh"#include "fastjet/CMSIterativeConePlugin.hh"
Go to the source code of this file.
Functions | |
| double | pow2 (const double x) |
| void | print_jets_and_sub (fj::ClusterSequence &clust_seq, const vector< fj::PseudoJet > &jets, double dcut) |
| a function that pretty prints a list of jets and the subjets for each one | |
| void | print_jets (const vector< fj::PseudoJet > &jets, const fj::ClusterSequence &cs, bool show_const=false) |
| void | is_unavailable (const string &algname) |
| int | main (int argc, char **argv) |
| a program to test and time a range of algorithms as implemented or wrapped in fastjet | |
| void | print_jet (const fj::ClusterSequence &clust_seq, const fj::PseudoJet &jet) |
| print a single jet | |
Variables | |
| string | rootfile |
| CmdLine * | cmdline_p |
| bool | do_areas |
| bool | ee_print = false |
| sort and pretty print jets, with exact behaviour depending on whether ee_print is true or not | |
| void is_unavailable | ( | const string & | algname | ) |
Definition at line 240 of file fastjet_timing_plugins.cc.
Referenced by main().
{
cerr << algname << " requested, but not available for this compilation";
exit(0);
}
| int main | ( | int | argc, |
| char ** | argv | ||
| ) |
a program to test and time a range of algorithms as implemented or wrapped in fastjet
Definition at line 248 of file fastjet_timing_plugins.cc.
References fastjet::active_area, fastjet::active_area_explicit_ghosts, fastjet::GhostedAreaSpec::add_ghosts(), CmdLine::all_options_used(), fastjet::antikt_algorithm, Best, fastjet::cambridge_algorithm, cmdline_p, fastjet::AreaDefinition::description(), fastjet::JetDefinition::description(), do_areas, CmdLine::double_val(), fastjet::E_scheme, fastjet::ee_genkt_algorithm, fastjet::ee_kt_algorithm, ee_print, fastjet::fastjet_version_string(), fastjet::genkt_algorithm, CmdLine::int_val(), is_unavailable(), fastjet::kt_algorithm, fastjet::SISConeBaseExtras::most_ambiguous_split(), fastjet::SISConeBaseExtras::pass(), fastjet::passive_area, fastjet::d0::inline_maths::phi(), fastjet::JetDefinition::plugin(), fastjet::plugin_strategy, pow2(), CmdLine::present(), print_jets(), print_jets_and_sub(), fastjet::PseudoJet::rap(), rootfile, fastjet::SISConeBasePlugin::set_ghost_separation_scale(), fastjet::GhostedAreaSpec::set_grid_scatter(), fastjet::sorted_by_pt(), fastjet::sorted_by_rapidity(), fastjet::SISConeBaseExtras::stable_cones(), fastjet::JetDefinition::strategy(), CmdLine::value(), and fastjet::voronoi_area.
{
CmdLine cmdline(argc,argv);
cmdline_p = &cmdline;
// allow the use to specify the fj::Strategy either through the
// -clever or the -strategy options (both will take numerical
// values); the latter will override the former.
fj::Strategy strategy = fj::Strategy(cmdline.int_val("-strategy",
cmdline.int_val("-clever", fj::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);
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 etamax = cmdline.double_val("-etamax",1.0e305);
bool show_constituents = cmdline.present("-const");
bool massless = cmdline.present("-massless");
int nev = cmdline.int_val("-nev",1);
bool add_dense_coverage = cmdline.present("-dense");
double ghost_maxrap = cmdline.value("-ghost-maxrap",5.0);
do_areas = cmdline.present("-area");
fj::AreaDefinition area_def;
if (do_areas) {
assert(!write); // it's incompatible
fj::GhostedAreaSpec ghost_spec(ghost_maxrap,
cmdline.value("-area:repeat", 1),
cmdline.value("-ghost-area", 0.01));
cmdline.present("-area:fj2"); // we only have fj2 areas, but allow the option anyway
if (cmdline.present("-area:explicit")) {
area_def = fj::AreaDefinition(fj::active_area_explicit_ghosts, ghost_spec);
} else if (cmdline.present("-area:passive")) {
area_def = fj::AreaDefinition(fj::passive_area, ghost_spec);
} else if (cmdline.present("-area:voronoi")) {
double Rfact = cmdline.value<double>("-area:voronoi");
area_def = fj::AreaDefinition(fj::voronoi_area,
fj::VoronoiAreaSpec(Rfact));
} else {
cmdline.present("-area:active"); // allow, but do not require, arg
area_def = fj::AreaDefinition(fj::active_area, ghost_spec);
}
}
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);
// 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
fj::RecombinationScheme scheme = fj::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"
fj::JetDefinition jet_def;
if (cmdline.present("-cam") || cmdline.present("-CA")) {
jet_def = fj::JetDefinition(fj::cambridge_algorithm, ktR, scheme, strategy);
} else if (cmdline.present("-antikt")) {
jet_def = fj::JetDefinition(fj::antikt_algorithm, ktR, scheme, strategy);
} else if (cmdline.present("-genkt")) {
double p = cmdline.value<double>("-genkt");
jet_def = fj::JetDefinition(fj::genkt_algorithm, ktR, p, scheme, strategy);
} else if (cmdline.present("-eekt")) {
jet_def = fj::JetDefinition(fj::ee_kt_algorithm);
} else if (cmdline.present("-eegenkt")) {
double p = cmdline.value<double>("-eegenkt");
jet_def = fj::JetDefinition(fj::ee_genkt_algorithm, ktR, p, scheme, strategy);
// checking if one asks to run a plugin (don't delete this line)
} else if (cmdline.present("-midpoint")) {
#ifdef ENABLE_PLUGIN_CDFCONES
typedef fj::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_def = fj::JetDefinition( new fj::CDFMidPointPlugin (
seed_threshold, ktR,
cone_area_fraction, max_pair_size,
max_iterations, overlap_threshold,
sm_scale));
#else // ENABLE_PLUGIN_CDFCONES
is_unavailable("midpoint");
#endif // ENABLE_PLUGIN_CDFCONES
} else if (cmdline.present("-pxcone")) {
#ifdef ENABLE_PLUGIN_PXCONE
double min_jet_energy = 5.0;
jet_def = fj::JetDefinition( new fj::PxConePlugin (
ktR, min_jet_energy,
overlap_threshold));
#else // ENABLE_PLUGIN_PXCONE
is_unavailable("pxcone");
#endif // ENABLE_PLUGIN_PXCONE
} else if (cmdline.present("-jetclu")) {
#ifdef ENABLE_PLUGIN_CDFCONES
jet_def = fj::JetDefinition( new fj::CDFJetCluPlugin (
ktR, overlap_threshold, seed_threshold));
#else // ENABLE_PLUGIN_CDFCONES
is_unavailable("pxcone");
#endif // ENABLE_PLUGIN_CDFCONES
} else if (cmdline.present("-siscone") || cmdline.present("-sisconespheri")) {
#ifdef ENABLE_PLUGIN_SISCONE
typedef fj::SISConePlugin SISPlug; // for brevity
int npass = cmdline.value("-npass",0);
if (cmdline.present("-siscone")) {
double sisptmin = cmdline.value("-sisptmin",0.0);
SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass,sisptmin);
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_def = fj::JetDefinition(plugin);
} else {
double sisEmin = cmdline.value("-sisEmin",0.0);
fj::SISConeSphericalPlugin * plugin =
new fj::SISConeSphericalPlugin(ktR, overlap_threshold,npass,sisEmin);
if (cmdline.present("-ghost-sep")) {
plugin->set_ghost_separation_scale(cmdline.value<double>("-ghost-sep"));
}
jet_def = fj::JetDefinition(plugin);
}
#else // ENABLE_PLUGIN_SISCONE
is_unavailable("siscone");
#endif // ENABLE_PLUGIN_SISCONE
} else if (cmdline.present("-d0runiicone")) {
#ifdef ENABLE_PLUGIN_D0RUNIICONE
double min_jet_Et = 6.0; // was 8 GeV in earlier work
jet_def = fj::JetDefinition(new fj::D0RunIIConePlugin(ktR,min_jet_Et));
#else // ENABLE_PLUGIN_D0RUNIICONE
is_unavailable("D0RunIICone");
#endif // ENABLE_PLUGIN_D0RUNIICONE
} else if (cmdline.present("-trackjet")) {
#ifdef ENABLE_PLUGIN_TRACKJET
jet_def = fj::JetDefinition(new fj::TrackJetPlugin(ktR));
#else // ENABLE_PLUGIN_TRACKJET
is_unavailable("TrackJet");
#endif // ENABLE_PLUGIN_TRACKJET
} else if (cmdline.present("-atlascone")) {
#ifdef ENABLE_PLUGIN_ATLASCONE
jet_def = fj::JetDefinition(new fj::ATLASConePlugin(ktR));
#else // ENABLE_PLUGIN_ATLASCONE
is_unavailable("ATLASCone");
#endif // ENABLE_PLUGIN_ATLASCONE
} else if (cmdline.present("-eecambridge")) {
#ifdef ENABLE_PLUGIN_EECAMBRIDGE
jet_def = fj::JetDefinition(new fj::EECambridgePlugin(ycut));
#else // ENABLE_PLUGIN_EECAMBRIDGE
is_unavailable("EECambridge");
#endif // ENABLE_PLUGIN_EECAMBRIDGE
} else if (cmdline.present("-jade")) {
#ifdef ENABLE_PLUGIN_JADE
jet_def = fj::JetDefinition(new fj::JadePlugin());
#else // ENABLE_PLUGIN_JADE
is_unavailable("Jade");
#endif // ENABLE_PLUGIN_JADE
} else if (cmdline.present("-cmsiterativecone")) {
#ifdef ENABLE_PLUGIN_CMSITERATIVECONE
jet_def = fj::JetDefinition(new fj::CMSIterativeConePlugin(ktR,seed_threshold));
#else // ENABLE_PLUGIN_CMSITERATIVECONE
is_unavailable("CMSIterativeCone");
#endif // ENABLE_PLUGIN_CMSITERATIVECONE
// end of checking if one asks to run a plugin (don't delete this line)
} else {
cmdline.present("-kt"); // kt is default, but allow user to specify it too [and ignore return value!]
jet_def = fj::JetDefinition(fj::kt_algorithm, ktR, strategy);
}
if (!cmdline.all_options_used()) {cerr <<
"Error: some options were not recognized"<<endl;
exit(-1);}
for (int iev = 0; iev < nev; iev++) {
vector<fj::PseudoJet> jets;
string line;
int ndone = 0;
while (getline(cin, 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];
}
}
fj::PseudoJet psjet(fourvec);
if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
}
// 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) {
fj::GhostedAreaSpec ghosted_area_spec(ghost_maxrap);
//fj::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(jets);
//----- 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) * (fj::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";
// fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
// jets.push_back(mom);
// }
// }
}
for (int irepeat = 0; irepeat < repeat ; irepeat++) {
int nparticles = jets.size();
try {
auto_ptr<fj::ClusterSequence> clust_seq;
if (do_areas) {
clust_seq.reset(new fj::ClusterSequenceArea(jets,jet_def,area_def));
} else {
clust_seq.reset(new fj::ClusterSequence(jets,jet_def,write));
}
if (irepeat != 0) {continue;}
cout << "iev "<<iev<< ": number of particles = "<< nparticles << endl;
cout << "strategy used = "<< clust_seq->strategy_string()<< endl;
cout << "Algorithm: " << jet_def.description() << " (" << fj::fastjet_version_string() << ")" << endl;
if (do_areas && iev == 0) cout << "Area definition: " << area_def.description() << endl;
// now provide some nice output...
if (inclkt >= 0.0) {
vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq->inclusive_jets(inclkt));
print_jets(jets, *clust_seq, show_constituents);
}
if (excln > 0) {
cout << "Printing "<<excln<<" exclusive jets\n";
print_jets(clust_seq->exclusive_jets(excln), *clust_seq, show_constituents);
}
if (excld > 0.0) {
cout << "Printing exclusive jets for d = "<<excld<<"\n";
print_jets(clust_seq->exclusive_jets(excld), *clust_seq, show_constituents);
}
if (excly > 0.0) {
cout << "Printing exclusive jets for ycut = "<<excly<<"\n";
print_jets(clust_seq->exclusive_jets_ycut(excly), *clust_seq, 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, 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++) {
fj::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 ENABLE_PLUGIN_SISCONE
// provide some complementary information for SISCone
if (show_cones) {
const fj::SISConeExtras * extras =
dynamic_cast<const fj::SISConeExtras *>(clust_seq->extras());
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<fj::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 %8d\n",
sisjets[i].rap(), sisjets[i].phi(), sisjets[i].perp(),
sisjets[i].user_index(), extras->pass(sisjets[i]),
clust_seq->constituents(sisjets[i]).size()
);
}
}
#endif // ENABLE_PLUGIN_SISCONE
} // try
catch (fastjet::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()==fj::plugin_strategy){
delete jet_def.plugin();
}
}
| double pow2 | ( | const double | x | ) | [inline] |
Definition at line 224 of file fastjet_timing_plugins.cc.
{return x*x;}
| void print_jet | ( | const fj::ClusterSequence & | clust_seq, |
| const fj::PseudoJet & | jet | ||
| ) |
print a single jet
Definition at line 637 of file fastjet_timing_plugins.cc.
References fastjet::ClusterSequence::constituents(), fastjet::PseudoJet::perp(), fastjet::PseudoJet::phi(), and fastjet::PseudoJet::rap().
{
int n_constituents = clust_seq.constituents(jet).size();
printf("%15.8f %15.8f %15.8f %8u\n",
jet.rap(), jet.phi(), jet.perp(), n_constituents);
}
| void print_jets | ( | const vector< fj::PseudoJet > & | jets, |
| const fj::ClusterSequence & | cs, | ||
| bool | show_const = false |
||
| ) |
Definition at line 646 of file fastjet_timing_plugins.cc.
References fastjet::ClusterSequenceAreaBase::area(), fastjet::ClusterSequenceAreaBase::area_4vector(), cmdline_p, CmdLine::command_line(), fastjet::ClusterSequence::constituents(), do_areas, ee_print, fastjet::d0::inline_maths::phi(), fastjet::ClusterSequence::print_jets_for_root(), rootfile, fastjet::sorted_by_E(), and fastjet::sorted_by_pt().
{
vector<fj::PseudoJet> jets;
if (ee_print) {
jets = sorted_by_E(jets_in);
for (size_t 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<fj::PseudoJet> const_jets = cs.constituents(jets[j]);
for (size_t 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 (size_t 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?)
if (do_areas) {
const fj::ClusterSequenceAreaBase * csab =
dynamic_cast<const fj::ClusterSequenceAreaBase *>(&cs);
printf(" %15.8f %15.8f", csab->area(jets[j]),
csab->area_4vector(jets[j]).perp());
}
cout << "\n";
if (show_constituents) {
vector<fj::PseudoJet> const_jets = cs.constituents(jets[j]);
for (size_t k = 0; k < const_jets.size(); k++) {
printf(" jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
const_jets[k].phi(),sqrt(const_jets[k].kt2()));
}
cout << "\n\n";
}
}
}
if (rootfile != "") {
ofstream ostr(rootfile.c_str());
ostr << "# " << cmdline_p->command_line() << endl;
ostr << "# output for root" << endl;
cs.print_jets_for_root(jets,ostr);
}
}
| void print_jets_and_sub | ( | fj::ClusterSequence & | clust_seq, |
| const vector< fj::PseudoJet > & | jets, | ||
| double | dcut | ||
| ) |
a function that pretty prints a list of jets and the subjets for each one
Definition at line 702 of file fastjet_timing_plugins.cc.
References fastjet::cambridge_algorithm, fastjet::ClusterSequence::constituents(), fastjet::ClusterSequence::exclusive_jets(), fastjet::ClusterSequence::exclusive_subdmerge_max(), fastjet::ClusterSequence::exclusive_subjets(), fastjet::ClusterSequence::inclusive_jets(), fastjet::ClusterSequence::jet_def(), fastjet::kt_algorithm, print_jet(), and fastjet::sorted_by_pt().
{
// sort jets into increasing pt
vector<fj::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");
// have various kinds of subjet finding, to test consistency among them
enum SubType {internal, newclust_dcut, newclust_R};
SubType subtype = internal;
//SubType subtype = newclust_dcut;
// print out the details for each jet
for (unsigned int i = 0; i < sorted_jets.size(); i++) {
// if jet pt^2 < dcut with kt alg, then some methods of
// getting subjets will return nothing -- so skip the jet
if (clust_seq.jet_def().jet_algorithm() == fj::kt_algorithm
&& sorted_jets[i].perp2() < dcut) continue;
printf("%5u ",i);
print_jet(clust_seq, sorted_jets[i]);
vector<fj::PseudoJet> subjets;
fj::ClusterSequence * cspoint;
if (subtype == internal) {
cspoint = &clust_seq;
subjets = clust_seq.exclusive_subjets(sorted_jets[i], dcut);
//subjets = clust_seq.exclusive_subjets(sorted_jets[i], 5);
double ddnp1 = clust_seq.exclusive_subdmerge_max(sorted_jets[i], subjets.size());
double ddn = clust_seq.exclusive_subdmerge_max(sorted_jets[i], subjets.size()-1);
cout << " for " << ddnp1 << " < d < " << ddn << " one has " << endl;
//subjets = clust_seq.exclusive_subjets(sorted_jets[i], dd*1.0000001);
} else if (subtype == newclust_dcut) {
cspoint = new fj::ClusterSequence(clust_seq.constituents(sorted_jets[i]),
clust_seq.jet_def());
subjets = cspoint->exclusive_jets(dcut);
//subjets = cspoint->exclusive_jets(int(min(5U,cspoint->n_particles())));
} else if (subtype == newclust_R) {
assert(clust_seq.jet_def().jet_algorithm() == fj::cambridge_algorithm);
fj::JetDefinition subjd(clust_seq.jet_def().jet_algorithm(),
clust_seq.jet_def().R()*sqrt(dcut));
cspoint = new fj::ClusterSequence(clust_seq.constituents(sorted_jets[i]),
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(*cspoint, subjets[j]);
}
if (cspoint != &clust_seq) delete cspoint;
//fj::ClusterSequence subseq(clust_seq.constituents(sorted_jets[i]),
// fj::JetDefinition(fj::cambridge_algorithm, 0.4));
//vector<fj::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]);
//}
}
}
Definition at line 231 of file fastjet_timing_plugins.cc.
Referenced by main(), and print_jets().
| bool do_areas |
Definition at line 233 of file fastjet_timing_plugins.cc.
Referenced by main(), and print_jets().
| bool ee_print = false |
sort and pretty print jets, with exact behaviour depending on whether ee_print is true or not
Definition at line 237 of file fastjet_timing_plugins.cc.
Referenced by main(), and print_jets().
| string rootfile |
Definition at line 230 of file fastjet_timing_plugins.cc.
Referenced by main(), and print_jets().
1.7.4