193 #include "fastjet/ClusterSequenceArea.hh"
194 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
195 #include "fastjet/tools/GridMedianBackgroundEstimator.hh"
196 #include "fastjet/Selector.hh"
204 #include "CmdLine.hh"
207 #include "fastjet/config.h"
210 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
211 #include "fastjet/SISConePlugin.hh"
212 #include "fastjet/SISConeSphericalPlugin.hh"
214 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
215 #include "fastjet/CDFMidPointPlugin.hh"
216 #include "fastjet/CDFJetCluPlugin.hh"
218 #ifdef FASTJET_ENABLE_PLUGIN_PXCONE
219 #include "fastjet/PxConePlugin.hh"
221 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
222 #include "fastjet/D0RunIIConePlugin.hh"
224 #ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
225 #include "fastjet/TrackJetPlugin.hh"
227 #ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
228 #include "fastjet/ATLASConePlugin.hh"
230 #ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
231 #include "fastjet/EECambridgePlugin.hh"
233 #ifdef FASTJET_ENABLE_PLUGIN_JADE
234 #include "fastjet/JadePlugin.hh"
236 #ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
237 #include "fastjet/CMSIterativeConePlugin.hh"
239 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
240 #include "fastjet/D0RunIpre96ConePlugin.hh"
241 #include "fastjet/D0RunIConePlugin.hh"
243 #ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
244 #include "fastjet/GridJetPlugin.hh"
252 namespace fj = fastjet;
254 inline double pow2(
const double x) {
return x*x;}
257 void print_jets_and_sub (
const vector<fj::PseudoJet> & jets,
double dcut);
266 bool ee_print =
false;
267 void print_jets(
const vector<fj::PseudoJet> & jets,
bool show_const =
false);
269 bool found_unavailable =
false;
270 void is_unavailable(
const string & algname) {
271 cerr << algname <<
" requested, but not available for this compilation" << endl;
272 found_unavailable =
true;
279 int main (
int argc,
char ** argv) {
283 CmdLine cmdline(argc,argv);
284 cmdline_p = &cmdline;
289 cmdline.int_val(
"-clever", fj::Best)));
290 int repeat = cmdline.int_val(
"-repeat",1);
291 int combine = cmdline.int_val(
"-combine",1);
292 bool write = cmdline.present(
"-write");
293 bool unique_write = cmdline.present(
"-unique_write");
294 bool hydjet = cmdline.present(
"-hydjet");
295 double ktR = cmdline.double_val(
"-r",1.0);
296 ktR = cmdline.double_val(
"-R",ktR);
297 double inclkt = cmdline.double_val(
"-incl",-1.0);
298 double repeatinclkt = cmdline.double_val(
"-repeat-incl",-1.0);
299 int excln = cmdline.int_val (
"-excln",-1);
300 double excld = cmdline.double_val(
"-excld",-1.0);
301 double excly = cmdline.double_val(
"-excly",-1.0);
302 ee_print = cmdline.present(
"-ee-print");
303 bool get_all_dij = cmdline.present(
"-get-all-dij");
304 bool get_all_yij = cmdline.present(
"-get-all-yij");
305 double subdcut = cmdline.double_val(
"-subdcut",-1.0);
306 double etamax = cmdline.double_val(
"-etamax",1.0e305);
307 bool show_constituents = cmdline.present(
"-const");
308 bool massless = cmdline.present(
"-massless");
309 int nev = cmdline.int_val(
"-nev",1);
310 bool add_dense_coverage = cmdline.present(
"-dense");
311 double ghost_maxrap = cmdline.value(
"-ghost-maxrap",5.0);
312 bool all_algs = cmdline.present(
"-all-algs");
314 fj::Selector particles_sel = (cmdline.present(
"-nhardest"))
316 : fj::SelectorIdentity();
318 do_areas = cmdline.present(
"-area");
323 cmdline.value(
"-area:repeat", 1),
324 cmdline.value(
"-ghost-area", 0.01));
326 if (cmdline.present(
"-area:explicit")) {
328 }
else if (cmdline.present(
"-area:passive")) {
330 }
else if (cmdline.present(
"-area:voronoi")) {
331 double Rfact = cmdline.value<
double>(
"-area:voronoi");
335 cmdline.present(
"-area:active");
339 bool do_bkgd = cmdline.present(
"-bkgd");
340 bool do_bkgd_csab =
false, do_bkgd_jetmedian =
false, do_bkgd_fj2 =
false;
341 bool do_bkgd_gridmedian =
false;
345 if (cmdline.present(
"-bkgd:csab")) {do_bkgd_csab =
true;}
346 else if (cmdline.present(
"-bkgd:jetmedian")) {do_bkgd_jetmedian =
true;
347 do_bkgd_fj2 = cmdline.present(
"-bkgd:fj2");
348 }
else if (cmdline.present(
"-bkgd:gridmedian")) {do_bkgd_gridmedian =
true;
350 throw fj::Error(
"with the -bkgd option, some particular background must be specified (csab or jetmedian)");
352 assert(do_areas || do_bkgd_gridmedian);
355 bool show_cones = cmdline.present(
"-cones");
359 double overlap_threshold = cmdline.double_val(
"-overlap",0.5);
360 overlap_threshold = cmdline.double_val(
"-f",overlap_threshold);
361 double seed_threshold = cmdline.double_val(
"-seed",1.0);
364 double ycut = cmdline.double_val(
"-ycut",0.08);
367 rootfile = cmdline.value<
string>(
"-root",
"");
375 vector<fj::JetDefinition> jet_defs;
376 if (all_algs || cmdline.present(
"-cam") || cmdline.present(
"-CA")) {
377 jet_defs.push_back(
fj::JetDefinition(fj::cambridge_algorithm, ktR, scheme, strategy));
379 if (all_algs || cmdline.present(
"-antikt")) {
380 jet_defs.push_back(
fj::JetDefinition(fj::antikt_algorithm, ktR, scheme, strategy));
382 if (all_algs || cmdline.present(
"-genkt")) {
384 if (cmdline.present(
"-genkt")) p = cmdline.value<
double>(
"-genkt");
386 jet_defs.push_back(
fj::JetDefinition(fj::genkt_algorithm, ktR, p, scheme, strategy));
388 if (all_algs || cmdline.present(
"-eekt")) {
391 if (all_algs || cmdline.present(
"-eegenkt")) {
393 if (cmdline.present(
"-eegenkt")) p = cmdline.value<
double>(
"-eegenkt");
395 jet_defs.push_back(
fj::JetDefinition(fj::ee_genkt_algorithm, ktR, p, scheme, strategy));
399 if (all_algs || cmdline.present(
"-midpoint")) {
400 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
401 typedef fj::CDFMidPointPlugin MPPlug;
402 double cone_area_fraction = 1.0;
403 int max_pair_size = 2;
404 int max_iterations = 100;
405 MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
406 if (cmdline.present(
"-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
407 if (cmdline.present(
"-sm-pt")) sm_scale = MPPlug::SM_pt;
408 if (cmdline.present(
"-sm-mt")) sm_scale = MPPlug::SM_mt;
409 if (cmdline.present(
"-sm-Et")) sm_scale = MPPlug::SM_Et;
412 cone_area_fraction, max_pair_size,
413 max_iterations, overlap_threshold,
415 #else // FASTJET_ENABLE_PLUGIN_CDFCONES
416 is_unavailable(
"midpoint");
417 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
419 if (all_algs || cmdline.present(
"-pxcone")) {
420 #ifdef FASTJET_ENABLE_PLUGIN_PXCONE
421 double min_jet_energy = 5.0;
424 overlap_threshold)));
425 #else // FASTJET_ENABLE_PLUGIN_PXCONE
426 is_unavailable(
"pxcone");
427 #endif // FASTJET_ENABLE_PLUGIN_PXCONE
429 if (all_algs || cmdline.present(
"-jetclu")) {
430 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
432 ktR, overlap_threshold, seed_threshold)));
433 #else // FASTJET_ENABLE_PLUGIN_CDFCONES
434 is_unavailable(
"pxcone");
435 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
437 if (all_algs || cmdline.present(
"-siscone") || cmdline.present(
"-sisconespheri")) {
438 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
440 int npass = cmdline.value(
"-npass",0);
441 if (all_algs || cmdline.present(
"-siscone")) {
442 double sisptmin = cmdline.value(
"-sisptmin",0.0);
443 SISPlug * plugin =
new SISPlug (ktR, overlap_threshold,npass,sisptmin);
444 if (cmdline.present(
"-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
445 if (cmdline.present(
"-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
446 if (cmdline.present(
"-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
447 if (cmdline.present(
"-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
449 plugin->set_use_jet_def_recombiner(
true);
452 if (all_algs || cmdline.present(
"-sisconespheri")) {
453 double sisEmin = cmdline.value(
"-sisEmin",0.0);
456 if (cmdline.present(
"-ghost-sep")) {
457 plugin->set_ghost_separation_scale(cmdline.value<
double>(
"-ghost-sep"));
461 #else // FASTJET_ENABLE_PLUGIN_SISCONE
462 is_unavailable(
"siscone");
463 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
465 if (all_algs || cmdline.present(
"-d0runiicone")) {
466 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
467 double min_jet_Et = 6.0;
469 #else // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
470 is_unavailable(
"D0RunIICone");
471 #endif // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
473 if (all_algs || cmdline.present(
"-trackjet")) {
474 #ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
476 #else // FASTJET_ENABLE_PLUGIN_TRACKJET
477 is_unavailable(
"TrackJet");
478 #endif // FASTJET_ENABLE_PLUGIN_TRACKJET
480 if (all_algs || cmdline.present(
"-atlascone")) {
481 #ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
483 #else // FASTJET_ENABLE_PLUGIN_ATLASCONE
484 is_unavailable(
"ATLASCone");
485 #endif // FASTJET_ENABLE_PLUGIN_ATLASCONE
487 if (all_algs || cmdline.present(
"-eecambridge")) {
488 #ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
490 #else // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
491 is_unavailable(
"EECambridge");
492 #endif // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
494 if (all_algs || cmdline.present(
"-jade")) {
495 #ifdef FASTJET_ENABLE_PLUGIN_JADE
497 #else // FASTJET_ENABLE_PLUGIN_JADE
498 is_unavailable(
"Jade");
499 #endif // FASTJET_ENABLE_PLUGIN_JADE
501 if (all_algs || cmdline.present(
"-cmsiterativecone")) {
502 #ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
504 #else // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
505 is_unavailable(
"CMSIterativeCone");
506 #endif // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
508 if (all_algs || cmdline.present(
"-d0runipre96cone")) {
509 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
510 jet_defs.push_back(
fj::JetDefinition(
new fj::D0RunIpre96ConePlugin(ktR, seed_threshold, overlap_threshold)));
511 #else // FASTJET_ENABLE_PLUGIN_D0RUNICONE
512 is_unavailable(
"D0RunICone");
513 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
515 if (all_algs || cmdline.present(
"-d0runicone")) {
516 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
517 jet_defs.push_back(
fj::JetDefinition(
new fj::D0RunIConePlugin(ktR, seed_threshold, overlap_threshold)));
518 #else // FASTJET_ENABLE_PLUGIN_D0RUNICONE
519 is_unavailable(
"D0RunICone");
520 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
522 if (all_algs || cmdline.present(
"-gridjet")) {
523 #ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
531 double grid_ymax = 4.9999999999;
533 #else // FASTJET_ENABLE_PLUGIN_GRIDJET
534 is_unavailable(
"GridJet");
535 #endif // FASTJET_ENABLE_PLUGIN_GRIDJET
539 cmdline.present(
"-kt") ||
540 (jet_defs.size() == 0 && !found_unavailable)) {
544 string filename = cmdline.value<
string>(
"-file",
"");
547 if (!cmdline.all_options_used()) {cerr <<
548 "Error: some options were not recognized"<<endl;
551 for (
unsigned idef = 0; idef < jet_defs.size(); idef++) {
554 if (filename ==
"") istr = &cin;
555 else istr =
new ifstream(filename.c_str());
557 for (
int iev = 0; iev < nev; iev++) {
558 vector<fj::PseudoJet> jets;
559 vector<fj::PseudoJet> particles;
562 while (getline(*istr, line)) {
564 istringstream linestream(line);
565 if (line ==
"#END") {
567 if (ndone == combine) {
break;}
569 if (line.substr(0,1) ==
"#") {
continue;}
570 valarray<double> fourvec(4);
575 int ii, istat,id,m1,m2,d1,d2;
577 linestream >> ii>> istat >>
id >> m1 >> m2 >> d1 >> d2
578 >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
581 fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
582 +pow2(fourvec[2])+pow2(mass));
586 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
587 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
589 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
593 if (abs(psjet.rap() < etamax)) {particles.push_back(psjet);}
599 if (add_dense_coverage) {
604 ghosted_area_spec.set_grid_scatter(0.5);
605 ghosted_area_spec.add_ghosts(particles);
628 particles = particles_sel(particles);
630 for (
int irepeat = 0; irepeat < repeat ; irepeat++) {
631 int nparticles = particles.size();
633 auto_ptr<fj::ClusterSequence> clust_seq;
641 if (repeatinclkt >= 0.0) {
642 vector<fj::PseudoJet> jets_local =
sorted_by_pt(clust_seq->inclusive_jets(repeatinclkt));
645 if (irepeat != 0) {
continue;}
646 cout <<
"iev "<<iev<<
": number of particles = "<< nparticles << endl;
647 cout <<
"strategy used = "<< clust_seq->strategy_string()<< endl;
649 if (do_areas && iev == 0) cout <<
"Area definition: " << area_def.description() << endl;
653 vector<fj::PseudoJet> jets_local =
sorted_by_pt(clust_seq->inclusive_jets(inclkt));
659 cout <<
"Printing "<<excln<<
" exclusive jets\n";
660 print_jets(clust_seq->exclusive_jets(excln), show_constituents);
664 cout <<
"Printing exclusive jets for d = "<<excld<<
"\n";
665 print_jets(clust_seq->exclusive_jets(excld), show_constituents);
669 cout <<
"Printing exclusive jets for ycut = "<<excly<<
"\n";
670 print_jets(clust_seq->exclusive_jets_ycut(excly), show_constituents);
674 for (
int i = nparticles-1; i >= 0; i--) {
675 printf(
"d for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_dmerge(i));
679 for (
int i = nparticles-1; i >= 0; i--) {
680 printf(
"y for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_ymerge(i));
686 if (subdcut >= 0.0) {
687 print_jets_and_sub(clust_seq->inclusive_jets(), subdcut);
692 vector<int> unique_history = clust_seq->unique_history_order();
694 vector<int> inv_unique_history(clust_seq->history().size());
695 for (
unsigned int i = 0; i < unique_history.size(); i++) {
696 inv_unique_history[unique_history[i]] = i;}
698 for (
unsigned int i = 0; i < unique_history.size(); i++) {
699 fj::ClusterSequence::history_element el =
700 clust_seq->history()[unique_history[i]];
701 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
702 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
703 printf(
"%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
708 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
713 cout <<
"most ambiguous split (difference in squared dist) = "
715 vector<fastjet::PseudoJet> stable_cones(extras->
stable_cones());
717 for (
unsigned int i = 0; i < stable_cones.size(); i++) {
719 printf(
"%5u %15.8f %15.8f %15.8f\n",
720 i,stable_cones[i].rap(),stable_cones[i].phi(),
721 stable_cones[i].perp() );
726 vector<fj::PseudoJet> sisjets = clust_seq->inclusive_jets();
727 printf(
"\n%15s %15s %15s %12s %8s %8s\n",
"rap",
"phi",
"pt",
"user-index",
"pass",
"nconst");
728 for (
unsigned i = 0; i < sisjets.size(); i++) {
729 printf(
"%15.8f %15.8f %15.8f %12d %8d %8u\n",
730 sisjets[i].rap(), sisjets[i].phi(), sisjets[i].perp(),
731 sisjets[i].user_index(), extras->
pass(sisjets[i]),
732 (
unsigned int) clust_seq->constituents(sisjets[i]).size()
737 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
740 double rho, sigma, mean_area, empty_area, n_empty_jets;
747 }
else if (do_bkgd_jetmedian) {
749 bge.set_provide_fj2_sigma(do_bkgd_fj2);
750 bge.set_cluster_sequence(*csab);
753 mean_area = bge.mean_area();
754 empty_area = bge.empty_area();
755 n_empty_jets = bge.n_empty_jets();
757 assert(do_bkgd_gridmedian);
758 double rapmin, rapmax;
760 fj::GridMedianBackgroundEstimator bge(rapmax, 2*ktR);
761 bge.set_particles(particles);
764 mean_area = bge.mean_area();
768 cout <<
" rho = " << rho
769 <<
", sigma = " << sigma
770 <<
", mean_area = " << mean_area
771 <<
", empty_area = " << empty_area
772 <<
", n_empty_jets = " << n_empty_jets
777 cout <<
"Caught fastjet error, exiting gracefully" << endl;
788 if (istr != &cin)
delete istr;
799 unsigned int n_constituents = jet.
constituents().size();
800 printf(
"%15.8f %15.8f %15.8f %8u\n",
801 jet.
rap(), jet.
phi(), jet.
perp(), n_constituents);
806 void print_jets(
const vector<fj::PseudoJet> & jets_in,
bool show_constituents) {
807 vector<fj::PseudoJet> jets;
810 for (
unsigned int j = 0; j < jets.size(); j++) {
811 printf(
"%5u %15.8f %15.8f %15.8f %15.8f\n",
812 j,jets[j].px(),jets[j].py(),jets[j].pz(),jets[j].E());
813 if (show_constituents) {
814 vector<fj::PseudoJet> const_jets = jets[j].constituents();
815 for (
unsigned int k = 0; k < const_jets.size(); k++) {
816 printf(
" jet%03u %15.8f %15.8f %15.8f %15.8f\n",j,const_jets[k].px(),
817 const_jets[k].py(),const_jets[k].pz(),const_jets[k].E());
825 for (
unsigned int j = 0; j < jets.size(); j++) {
826 printf(
"%5u %15.8f %15.8f %15.8f",
827 j,jets[j].rap(),jets[j].phi(),jets[j].perp());
830 if (do_areas) printf(
" %15.8f %15.8f", jets[j].area(),
831 jets[j].area_4vector().perp());
834 if (show_constituents) {
835 vector<fj::PseudoJet> const_jets = jets[j].constituents();
836 for (
unsigned int k = 0; k < const_jets.size(); k++) {
837 printf(
" jet%03u %15.8f %15.8f %15.8f %5d\n",j,const_jets[k].rap(),
838 const_jets[k].phi(),sqrt(const_jets[k].kt2()), const_jets[k].cluster_hist_index());
845 if (rootfile !=
"") {
846 ofstream ostr(rootfile.c_str());
847 ostr <<
"# " << cmdline_p->command_line() << endl;
848 ostr <<
"# output for root" << endl;
849 assert(jets.size() > 0);
850 jets[0].validated_cs()->print_jets_for_root(jets,ostr);
859 void print_jets_and_sub (
const vector<fj::PseudoJet> & jets,
double dcut) {
865 printf(
"Printing jets and their subjets with subdcut = %10.5f\n",dcut);
866 printf(
"%5s %15s %15s %15s %15s\n",
"jet #",
"rapidity",
867 "phi",
"pt",
"n constituents");
870 enum SubType {
internal, newclust_dcut, newclust_R};
871 SubType subtype =
internal;
877 for (vector<fj::PseudoJet>::const_iterator jet = sorted_jets.begin();
878 jet != sorted_jets.end(); jet++) {
884 && jet->
perp2() < dcut)
continue;
887 printf(
"%5u ",(
unsigned int) (jet - sorted_jets.begin()));
889 vector<fj::PseudoJet> subjets;
891 if (subtype ==
internal) {
896 cout <<
" for " << ddnp1 <<
" < d < " << ddn <<
" one has " << endl;
897 }
else if (subtype == newclust_dcut) {
900 }
else if (subtype == newclust_R) {
906 cerr <<
"unrecognized subtype for subjet finding" << endl;
911 for (
unsigned int j = 0; j < subjets.size(); j++) {
912 printf(
" -sub-%02u ",j);
913 print_jet(subjets[j]);
916 if (cspoint != 0)
delete cspoint;