199#include "fastjet/ClusterSequenceArea.hh"
200#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
201#include "fastjet/tools/GridMedianBackgroundEstimator.hh"
202#include "fastjet/tools/Subtractor.hh"
203#include "fastjet/Selector.hh"
219#include "fastjet/config.h"
223#ifdef FASTJET_ENABLE_PLUGIN_SISCONE
224#include "fastjet/SISConePlugin.hh"
225#include "fastjet/SISConeSphericalPlugin.hh"
227#ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
228#include "fastjet/CDFMidPointPlugin.hh"
229#include "fastjet/CDFJetCluPlugin.hh"
231#ifdef FASTJET_ENABLE_PLUGIN_PXCONE
232#include "fastjet/PxConePlugin.hh"
234#ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
235#include "fastjet/D0RunIIConePlugin.hh"
237#ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
238#include "fastjet/TrackJetPlugin.hh"
240#ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
241#include "fastjet/ATLASConePlugin.hh"
243#ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
244#include "fastjet/EECambridgePlugin.hh"
246#ifdef FASTJET_ENABLE_PLUGIN_JADE
247#include "fastjet/JadePlugin.hh"
249#ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
250#include "fastjet/CMSIterativeConePlugin.hh"
252#ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
253#include "fastjet/D0RunIpre96ConePlugin.hh"
254#include "fastjet/D0RunIConePlugin.hh"
256#ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
257#include "fastjet/GridJetPlugin.hh"
267using namespace fjcore;
270inline double pow2(
const double x) {
return x*x;}
273void print_jets_and_sub (
const vector<PseudoJet> & jets,
double dcut);
276void print_jets_bkgd(
const vector<PseudoJet> &jets,
277 const vector<PseudoJet> &subtracted_jets,
288enum SubType {subtype_internal, subtype_newclust_dcut, subtype_newclust_R};
290void do_compare_strategy(
int iev,
291 const vector<PseudoJet> & particles,
294 int compare_strategy);
304bool ee_print =
false;
305void print_jets(
const vector<PseudoJet> & jets,
bool show_const =
false);
307bool found_unavailable =
false;
308void is_unavailable(
const string & algname) {
309 cerr << algname <<
" requested, but not available for this compilation" << endl;
310 found_unavailable =
true;
317int main (
int argc,
char ** argv) {
320 cout <<
"# " << cmdline.command_line() << endl;
323 cmdline_p = &cmdline;
328 cmdline.int_val(
"-clever",
Best)));
329 int repeat = cmdline.int_val(
"-repeat",1);
330 int combine = cmdline.int_val(
"-combine",1);
331 bool write = cmdline.present(
"-write");
332 bool unique_write = cmdline.present(
"-unique_write");
333 bool hydjet = cmdline.present(
"-hydjet");
334 double ktR = cmdline.double_val(
"-r",1.0);
335 ktR = cmdline.double_val(
"-R",ktR);
336 double inclkt = cmdline.double_val(
"-incl",-1.0);
337 double repeatinclkt = cmdline.double_val(
"-repeat-incl",-1.0);
338 int excln = cmdline.int_val (
"-excln",-1);
339 double excld = cmdline.double_val(
"-excld",-1.0);
340 double excly = cmdline.double_val(
"-excly",-1.0);
341 ee_print = cmdline.present(
"-ee-print");
342 bool get_all_dij = cmdline.present(
"-get-all-dij");
343 bool get_all_yij = cmdline.present(
"-get-all-yij");
344 double subdcut = cmdline.double_val(
"-subdcut",-1.0);
345 double rapmax = cmdline.double_val(
"-rapmax",1.0e305);
346 if (cmdline.present(
"-etamax")) {
347 cerr <<
"WARNING: -etamax options actually sets maximum rapidity (and overrides -rapmax)\n";
348 rapmax = cmdline.double_val(
"-etamax");
350 bool show_constituents = cmdline.present(
"-const");
351 bool massless = cmdline.present(
"-massless");
352 int nev = cmdline.int_val(
"-nev",1);
353 int skip = cmdline.int_val(
"-skip",0);
354 bool add_dense_coverage = cmdline.present(
"-dense");
355 double ghost_maxrap = cmdline.value(
"-ghost-maxrap",5.0);
356 bool all_algs = cmdline.present(
"-all-algs");
361 int compare_strategy = cmdline.value<
int>(
"-compare-strategy",
plugin_strategy);
364 Selector particles_sel = (cmdline.present(
"-nhardest"))
366 : SelectorIdentity();
368 do_areas = cmdline.present(
"-area");
374 cmdline.value(
"-area:repeat", 1),
375 cmdline.value(
"-ghost-area", 0.01));
376 if (cmdline.present(
"-area:fj2")) ghost_spec.set_fj2_placement(
true);
377 if (cmdline.present(
"-area:explicit")) {
378 area_def =
AreaDefinition(active_area_explicit_ghosts, ghost_spec);
379 }
else if (cmdline.present(
"-area:passive")) {
381 }
else if (cmdline.present(
"-area:voronoi")) {
382 double Rfact = cmdline.value<
double>(
"-area:voronoi");
386 cmdline.present(
"-area:active");
394 bool do_bkgd = cmdline.present(
"-bkgd");
396 bool do_bkgd_csab =
false, do_bkgd_jetmedian =
false, do_bkgd_fj2 =
false;
397 bool do_bkgd_gridmedian =
false;
398 bool do_bkgd_localrange =
false;
399 bool do_subtractor =
false;
400 double bkgd_alt_ktR = -1.0;
405 if (cmdline.present(
"-bkgd:csab")) {do_bkgd_csab =
true;}
406 else if (cmdline.present(
"-bkgd:jetmedian")) {do_bkgd_jetmedian =
true;
407 do_bkgd_fj2 = cmdline.present(
"-bkgd:fj2");
408 do_bkgd_localrange = cmdline.present(
"-bkgd:localrange");
409 bkgd_alt_ktR = cmdline.value(
"-bkgd:alt-ktR", bkgd_alt_ktR);
411 }
else if (cmdline.present(
"-bkgd:gridmedian")) {
412 do_bkgd_gridmedian =
true;
414 throw Error(
"with the -bkgd option, some particular background must be specified (csab or jetmedian)");
416 if (cmdline.present(
"-bkgd:rescaling")) {
419 assert(do_areas || do_bkgd_gridmedian);
420 do_subtractor = cmdline.present(
"-subtractor");
421 if (do_subtractor) assert(do_areas);
427 bool show_cones = cmdline.present(
"-cones");
431 double overlap_threshold = cmdline.double_val(
"-overlap",0.5);
432 overlap_threshold = cmdline.double_val(
"-f",overlap_threshold);
433 double seed_threshold = cmdline.double_val(
"-seed",1.0);
440 double ycut = cmdline.double_val(
"-ycut",0.08);
443 rootfile = cmdline.value<
string>(
"-root",
"");
451 vector<JetDefinition> jet_defs;
452 if (all_algs || cmdline.present(
"-cam") || cmdline.present(
"-CA")) {
455 if (all_algs || cmdline.present(
"-antikt")) {
458 if (all_algs || cmdline.present(
"-genkt")) {
460 if (cmdline.present(
"-genkt")) p = cmdline.value<
double>(
"-genkt");
464 if (all_algs || cmdline.present(
"-eekt")) {
467 if (all_algs || cmdline.present(
"-eegenkt")) {
469 if (cmdline.present(
"-eegenkt")) p = cmdline.value<
double>(
"-eegenkt");
475 if (all_algs || cmdline.present(
"-midpoint")) {
476#ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
478 double cone_area_fraction = 1.0;
479 int max_pair_size = 2;
480 int max_iterations = 100;
481 MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
482 if (cmdline.present(
"-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
483 if (cmdline.present(
"-sm-pt")) sm_scale = MPPlug::SM_pt;
484 if (cmdline.present(
"-sm-mt")) sm_scale = MPPlug::SM_mt;
485 if (cmdline.present(
"-sm-Et")) sm_scale = MPPlug::SM_Et;
488 cone_area_fraction, max_pair_size,
489 max_iterations, overlap_threshold,
492 is_unavailable(
"MidPoint");
495 if (all_algs || cmdline.present(
"-pxcone")) {
496#ifdef FASTJET_ENABLE_PLUGIN_PXCONE
497 double min_jet_energy = 5.0;
499 int mode = cmdline.value(
"-pxcone-mode", 2);
500 bool E_scheme_jets =
false;
503 overlap_threshold, E_scheme_jets, mode)));
505 is_unavailable(
"PxCone");
508 if (all_algs || cmdline.present(
"-jetclu")) {
509#ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
511 ktR, overlap_threshold, seed_threshold)));
513 is_unavailable(
"JetClu");
516 if (all_algs || cmdline.present(
"-siscone") || cmdline.present(
"-sisconespheri")) {
517#ifdef FASTJET_ENABLE_PLUGIN_SISCONE
519 int npass = cmdline.value(
"-npass",0);
520 if (all_algs || cmdline.present(
"-siscone")) {
521 double sisptmin = cmdline.value(
"-sisptmin",0.0);
522 bool cache = cmdline.present(
"-cache");
523 SISPlug * plugin =
new SISPlug (ktR, overlap_threshold,npass,sisptmin,cache);
524 if (cmdline.present(
"-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
525 if (cmdline.present(
"-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
526 if (cmdline.present(
"-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
527 if (cmdline.present(
"-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
529 plugin->set_use_jet_def_recombiner(
true);
532 if (all_algs || cmdline.present(
"-sisconespheri")) {
533 double sisEmin = cmdline.value(
"-sisEmin",0.0);
536 if (cmdline.present(
"-ghost-sep")) {
542 is_unavailable(
"SISCone");
545 if (all_algs || cmdline.present(
"-d0runiicone")) {
546#ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
547 double min_jet_Et = 6.0;
550 is_unavailable(
"D0RunIICone");
553 if (all_algs || cmdline.present(
"-trackjet")) {
554#ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
557 is_unavailable(
"TrackJet");
560 if (all_algs || cmdline.present(
"-atlascone")) {
561#ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
564 is_unavailable(
"ATLASCone");
567 if (all_algs || cmdline.present(
"-eecambridge")) {
568#ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
571 is_unavailable(
"EECambridge");
574 if (all_algs || cmdline.present(
"-jade")) {
575#ifdef FASTJET_ENABLE_PLUGIN_JADE
578 JadePlugin::strategy_NNFJN2Plain));
581 is_unavailable(
"Jade");
584 if (all_algs || cmdline.present(
"-cmsiterativecone")) {
585#ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
588 is_unavailable(
"CMSIterativeCone");
591 if (all_algs || cmdline.present(
"-d0runipre96cone")) {
592#ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
595 is_unavailable(
"D0RunIpre96Cone");
598 if (all_algs || cmdline.present(
"-d0runicone")) {
599#ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
602 is_unavailable(
"D0RunICone");
605 if (all_algs || cmdline.present(
"-gridjet")) {
606#ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
614 double grid_ymax = 4.9999999999;
617 is_unavailable(
"GridJet");
622 cmdline.present(
"-kt") ||
623 (jet_defs.size() == 0 && !found_unavailable)) {
627 string filename = cmdline.value<
string>(
"-file",
"");
630 if (!cmdline.all_options_used()) {cerr <<
631 "Error: some options were not recognized"<<endl;
634 for (
unsigned idef = 0; idef < jet_defs.size(); idef++) {
637 if (filename ==
"") istr = &cin;
638 else istr =
new ifstream(filename.c_str());
640 for (
int iev = 0; iev < nev; iev++) {
641 vector<PseudoJet> jets;
642 vector<PseudoJet> particles;
645 while (getline(*istr, line)) {
647 istringstream linestream(line);
648 if (line ==
"#END") {
650 if (ndone == combine) {
break;}
652 if (line.substr(0,1) ==
"#") {
continue;}
653 valarray<double> fourvec(4);
658 int ii, istat,id,m1,m2,d1,d2;
660 linestream >> ii>> istat >>
id >> m1 >> m2 >> d1 >> d2
661 >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
664 fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
665 +pow2(fourvec[2])+pow2(mass));
669 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
670 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
672 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
676 if (abs(psjet.rap()) < rapmax) {particles.push_back(psjet);}
683 if (add_dense_coverage) {
688 ghosted_area_spec.set_grid_scatter(0.5);
689 ghosted_area_spec.add_ghosts(particles);
711 add_dense_coverage =
false;
715 particles = particles_sel(particles);
718 if (iev < skip)
continue;
720 for (
int irepeat = 0; irepeat < repeat ; irepeat++) {
721 int nparticles = particles.size();
733 do_compare_strategy(iev, particles, jet_def, *clust_seq, compare_strategy);
737 if (repeatinclkt >= 0.0) {
738 vector<PseudoJet> jets_local =
sorted_by_pt(clust_seq->inclusive_jets(repeatinclkt));
741 if (irepeat != 0) {
continue;}
742 cout <<
"iev "<<iev<<
": number of particles = "<< nparticles << endl;
743 cout <<
"strategy used = "<< clust_seq->strategy_string()<< endl;
746 if (do_areas && iev == 0) cout <<
"Area definition: " << area_def.description() << endl;
751 vector<PseudoJet> jets_local =
sorted_by_pt(clust_seq->inclusive_jets(inclkt));
752 print_jets(jets_local, show_constituents);
757 cout <<
"Printing "<<excln<<
" exclusive jets\n";
758 print_jets(clust_seq->exclusive_jets(excln), show_constituents);
762 cout <<
"Printing exclusive jets for d = "<<excld<<
"\n";
763 print_jets(clust_seq->exclusive_jets(excld), show_constituents);
767 cout <<
"Printing exclusive jets for ycut = "<<excly<<
"\n";
768 print_jets(clust_seq->exclusive_jets_ycut(excly), show_constituents);
772 for (
int i = nparticles-1; i >= 0; i--) {
773 printf(
"d for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_dmerge(i));
777 for (
int i = nparticles-1; i >= 0; i--) {
778 printf(
"y for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_ymerge(i));
784 if (subdcut >= 0.0) {
785 print_jets_and_sub(clust_seq->inclusive_jets(), subdcut);
790 vector<int> unique_history = clust_seq->unique_history_order();
792 vector<int> inv_unique_history(clust_seq->history().size());
793 for (
unsigned int i = 0; i < unique_history.size(); i++) {
794 inv_unique_history[unique_history[i]] = i;}
796 for (
unsigned int i = 0; i < unique_history.size(); i++) {
798 clust_seq->history()[unique_history[i]];
801 printf(
"%7d u %15.8e %7d u %7d u\n",i,el.
dij,uhp1, uhp2);
806#ifdef FASTJET_ENABLE_PLUGIN_SISCONE
812 throw fastjet::Error(
"extras object for SISCone was null (this can happen with certain area types)");
813 cout <<
"most ambiguous split (difference in squared dist) = "
815 vector<fastjet::PseudoJet> stable_cones(extras->
stable_cones());
817 for (
unsigned int i = 0; i < stable_cones.size(); i++) {
819 printf(
"%5u %15.8f %15.8f %15.8f\n",
820 i,stable_cones[i].rap(),stable_cones[i].phi(),
821 stable_cones[i].perp() );
826 vector<PseudoJet> sisjets = clust_seq->inclusive_jets();
827 printf(
"\n%15s %15s %15s %12s %8s %8s\n",
"rap",
"phi",
"pt",
"user-index",
"pass",
"nconst");
828 for (
unsigned i = 0; i < sisjets.size(); i++) {
829 printf(
"%15.8f %15.8f %15.8f %12d %8d %8u\n",
830 sisjets[i].rap(), sisjets[i].phi(), sisjets[i].perp(),
831 sisjets[i].user_index(), extras->
pass(sisjets[i]),
832 (
unsigned int) clust_seq->constituents(sisjets[i]).size()
841 double rho, sigma, mean_area, empty_area, n_empty_jets;
849 }
else if (do_bkgd_jetmedian) {
855 if (bkgd_alt_ktR > 0) {
858 cout <<
"Alt JetDef for background-estimation CSAB: " << clust_seq_bkgd->
jet_def().
description() << endl;
864 if (!do_bkgd_localrange) {
866 sigma = bge->
sigma();
872 assert(do_bkgd_gridmedian);
873 double grid_rapmin, grid_rapmax;
880 sigma = bge->
sigma();
885 if (do_bkgd_localrange || do_subtractor) {
886 assert(bge_ptr != 0);
887 cout <<
"Background estimator: " << bge_ptr->
description() << endl;
890 vector<PseudoJet> subjets;
893 subtractor.set_use_rho_m(
true);
894 subtractor.set_safe_mass(
true);
895 cout <<
"Subtractor: " << subtractor.description() << endl;
896 subjets = subtractor(jets);
898 print_jets_bkgd(jets, subjets, bge_ptr, do_subtractor);
915 cout <<
" rho = " << rho
916 <<
", sigma = " << sigma
917 <<
", mean_area = " << mean_area
918 <<
", empty_area = " << empty_area
919 <<
", n_empty_jets = " << n_empty_jets
922 if (bge_ptr != 0)
delete bge_ptr;
926 catch (
Error &fjerr) {
927 cout <<
"Caught fastjet error, exiting gracefully" << endl;
938 if (istr != &cin)
delete istr;
949 unsigned int n_constituents = jet.
constituents().size();
950 printf(
"%15.8f %15.8f %15.8f %8u\n",
951 jet.
rap(), jet.
phi(), jet.
perp(), n_constituents);
956void print_jets(
const vector<PseudoJet> & jets_in,
bool show_constituents) {
957 vector<PseudoJet> jets;
960 for (
unsigned int j = 0; j < jets.size(); j++) {
961 printf(
"%5u %15.8f %15.8f %15.8f %15.8f\n",
962 j,jets[j].px(),jets[j].py(),jets[j].pz(),jets[j].E());
963 if (show_constituents) {
964 vector<PseudoJet> const_jets = jets[j].constituents();
965 for (
unsigned int k = 0; k < const_jets.size(); k++) {
966 printf(
" jet%03u %15.8f %15.8f %15.8f %15.8f\n",j,const_jets[k].px(),
967 const_jets[k].py(),const_jets[k].pz(),const_jets[k].E());
975 for (
unsigned int j = 0; j < jets.size(); j++) {
976 printf(
"%5u %15.8f %15.8f %15.8f",
977 j,jets[j].rap(),jets[j].phi(),jets[j].perp());
981 if (do_areas) printf(
" %15.8f %15.8f", jets[j].area(),
982 jets[j].area_4vector().perp());
986 if (show_constituents) {
987 vector<PseudoJet> const_jets = jets[j].constituents();
988 for (
unsigned int k = 0; k < const_jets.size(); k++) {
989 printf(
" jet%03u %15.8f %15.8f %15.8f %5d\n",j,const_jets[k].rap(),
990 const_jets[k].phi(),sqrt(const_jets[k].kt2()), const_jets[k].cluster_hist_index());
997 if (rootfile !=
"") {
998 ofstream ostr(rootfile.c_str());
1000 ostr <<
"# output for root" << endl;
1001 assert(jets.size() > 0);
1002 jets[0].validated_cs()->print_jets_for_root(jets,ostr);
1011void print_jets_and_sub (
const vector<PseudoJet> & jets,
double dcut) {
1017 printf(
"Printing jets and their subjets with subdcut = %10.5f\n",dcut);
1018 printf(
"%5s %15s %15s %15s %15s\n",
"jet #",
"rapidity",
1019 "phi",
"pt",
"n constituents");
1022 SubType sub_type = subtype_internal;
1028 for (vector<PseudoJet>::const_iterator jet = sorted_jets.begin();
1029 jet != sorted_jets.end(); jet++) {
1035 && jet->
perp2() < dcut)
continue;
1038 printf(
"%5u ",(
unsigned int) (jet - sorted_jets.begin()));
1040 vector<PseudoJet> subjets;
1042 if (sub_type == subtype_internal) {
1047 cout <<
" for " << ddnp1 <<
" < d < " << ddn <<
" one has " << endl;
1048 }
else if (sub_type == subtype_newclust_dcut) {
1051 }
else if (sub_type == subtype_newclust_R) {
1057 cerr <<
"unrecognized subtype for subjet finding" << endl;
1062 for (
unsigned int j = 0; j < subjets.size(); j++) {
1063 printf(
" -sub-%02u ",j);
1064 print_jet(subjets[j]);
1067 if (cspoint != 0)
delete cspoint;
1081double make_safe_zero_truncation(
double x,
double precision){
1082 return std::abs(x)<0.5*precision ? 0.0 : x;
1086void print_jets_bkgd(
const vector<PseudoJet> &jets,
1087 const vector<PseudoJet> &subtracted_jets,
1089 bool do_subtractor){
1090 printf(
"Printing jets, background information");
1092 printf(
" and subtracted jets\n");
1093 printf(
"%5s %15s %15s %15s %15s %15s %15s %15s %15s\n",
"jet #",
1094 "rapidity",
"phi",
"pt",
"pt^2+m^2",
1095 "rho",
"rho_m",
"sigma",
"sigma_m");
1097 printf(
"%5s %15s %15s %15s %15s %15ss\n",
"jet #",
1098 "rapidity",
"phi",
"pt",
"sqrt(pt^2+m^2)",
"area");
1100 for (
unsigned i = 0; i < jets.size(); i++) {
1107 printf(
"%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n", i,
1109 estimate.
rho(), make_safe_zero_truncation(estimate.
rho_m(),1e-8),
1111 if (do_subtractor) {
1112 const PseudoJet & subjet = subtracted_jets[i];
1113 printf(
"%5u %15.8f %15.8f %15.8f %15.8f %15.8f\n", i,
1121void signal_failed_comparison(
int iev,
1122 const string & message,
1123 const vector<PseudoJet> & particles) {
1124 cout <<
"# failed comparison, reason is " << message << endl;
1125 cout <<
"# iev = " << iev << endl;
1126 cout << setprecision(16);
1127 for (
unsigned i = 0; i < particles.size(); i++) {
1129 cout << p.px() <<
" "
1134 cout <<
"#END" << endl;
1138void do_compare_strategy(
int iev,
1139 const vector<PseudoJet> & particles,
1142 int compare_strategy) {
1147 jet_def.recombination_scheme(),
1154 const vector<ClusterSequence::history_element> & history_in = cs.
history();
1155 const vector<ClusterSequence::history_element> & history_ref = cs_ref.history();
1157 if (history_in.size() != history_ref.size()) {
1158 signal_failed_comparison(iev,
"history sizes do not match", particles);
1162 for (
unsigned i = cs.
n_particles(); i < history_in.size(); i++) {
1163 bool fail_parents = (history_in[i].parent1 != history_ref[i].parent1 ||
1164 history_in[i].parent2 != history_ref[i].parent2);
1165 bool fail_dij = (history_in[i].dij != history_ref[i].dij);
1166 if (fail_parents || fail_dij) {
1168 ostr <<
"at step " << i <<
", ";
1169 if (fail_parents) ostr <<
"parents do not match, ";
1170 if (fail_dij) ostr <<
"dij does not match, ";
1171 ostr <<
"history in (p1, p2, dij) = "
1172 << history_in[i].parent1 <<
" " << history_in[i].parent2 <<
" " << history_in[i].dij;
1173 ostr <<
", history ref (p1, p2, dij) = "
1174 << history_ref[i].parent1 <<
" " << history_ref[i].parent2 <<
" " << history_ref[i].dij;
1175 signal_failed_comparison(iev, ostr.str(), particles);
int main()
an example program showing how to use fastjet
string command_line() const
return the full command line
Implementation of the ATLAS Cone (plugin for fastjet v2.4 upwards)
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
Abstract base class that provides the basic interface for classes that estimate levels of background ...
virtual std::string description() const =0
returns a textual description of the background estimator
virtual BackgroundEstimate estimate() const =0
get the full set of background properties
A background rescaling that is a simple polynomial in y.
Implementation of the JetClu algorithm from CDF (plugin for fastjet-v2.1 upwards)
Implementation of the MidPoint algorithm from CDF (plugin for fastjet-v2.1 upwards)
Implementation of the CMS Iterative Cone (plugin for fastjet v2.4 upwards)
base class that sets interface for extensions of ClusterSequence that provide information about the a...
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),...
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 ...
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...
General class for user to obtain ClusterSequence with additional area information.
const JetDefinition & jet_def() const
return a reference to the jet definition
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
void delete_self_when_unused()
by calling this routine you tell the ClusterSequence to delete itself when all the Pseudojets associa...
unsigned int n_particles() const
returns the number of particles that were provided to the clustering algorithm (helps the user find t...
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.
static void print_banner()
This is the function that is automatically called during clustering to print the FastJet banner.
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...
A plugin for FastJet (v3.0 or later) that provides an interface to the D0 version of Run-I cone algor...
Implementation of the D0 Run II Cone (plugin for fastjet v2.1 upwards)
A plugin for FastJet (v3.0 or later) that provides an interface to the pre 1996 D0 version of Run-I c...
Implementation of the e+e- Cambridge algorithm (plugin for fastjet v2.4 upwards)
base class corresponding to errors that can be thrown by FastJet
Parameters to configure the computation of jet areas using ghosts.
plugin for fastjet (v3.0 upwards) that clusters particles such that all particles in a given cell of ...
Implementation of the e+e- Jade algorithm (plugin for fastjet v2.4 upwards)
Strategy
enum that contains the two clustering strategy options; for higher multiplicities,...
class that is intended to hold a full definition of the jet clusterer
std::string description() const
return a textual description of the current jet definition
const Plugin * plugin() const
return a pointer to the plugin
JetAlgorithm jet_algorithm() const
return information about the definition...
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
double mt() const
returns the transverse mass = sqrt(kt^2+m^2)
double rap() const
returns the rapidity or some large value when the rapidity is infinite
double phi() const
returns phi (in the range 0..2pi)
double perp() const
returns the scalar transverse momentum
double perp2() const
returns the squared transverse momentum
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...
virtual double area() const
return the jet (scalar) area.
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...
const ClusterSequence * validated_cs() const
shorthand for validated_cluster_sequence()
Implementation of the PxCone algorithm (plugin for fastjet v2.1 upwards)
virtual void set_ghost_separation_scale(double scale) const
set the ghost separation scale for passive area determinations just in the next run (strictly speakin...
Implementation of the SISCone algorithm (plugin for fastjet v2.1 upwards)
Implementation of the spherical version of the SISCone algorithm (plugin for fastjet v2....
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
void get_rapidity_extent(double &rapmin, double &rapmax) const
returns the rapidity range for which it may return "true"
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
T * get() const
get the stored pointer
void reset()
reset the pointer to default value (NULL)
Class that helps perform jet background subtraction.
Implementation of the TrackJet algorithm (plugin for fastjet v2.4 upwards)
Specification for the computation of the Voronoi jet area.
Selector SelectorNHardest(unsigned int n)
select the n hardest objects
Selector SelectorStrip(const double half_width)
select objets within a rapidity distance 'half_width' from the location of the reference jet,...
Selector SelectorAbsRapMax(double absrapmax)
select objects with |rap| <= absrapmax
vector< PseudoJet > sorted_by_rapidity(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing rapidity
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
@ Best
automatic selection of the best (based on N), including the LazyTiled strategies that are new to FJ3....
@ plugin_strategy
the plugin has been used...
RecombinationScheme
The various recombination schemes.
@ E_scheme
summing the 4-momenta
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
@ ee_genkt_algorithm
the e+e- genkt algorithm (R > 2 and p=1 gives ee_kt)
@ genkt_algorithm
like the k_t but with distance measures dij = min(kti^{2p},ktj^{2p}) Delta R_{ij}^2 / R^2 diB = 1/kti...
@ ee_kt_algorithm
the e+e- kt algorithm
@ cambridge_algorithm
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
@ antikt_algorithm
like the k_t but with distance measures dij = min(1/kti^2,1/ktj^2) Delta R_{ij}^2 / R^2 diB = 1/kti^2
@ kt_algorithm
the longitudinally invariant kt algorithm
string fastjet_version_string()
return a string containing information about the release
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
a single element in the clustering history
int parent1
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
int parent2
index in _history where second parent of this jet was created (InexistentParent if this jet is an ori...
double dij
the distance corresponding to the recombination at this stage of the clustering.