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"
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 rapmax = cmdline.double_val(
"-rapmax",1.0e305);
307 if (cmdline.present(
"-etamax")) {
308 cerr <<
"WARNING: -etamax options actually sets maximum rapidity (and overrides -rapmax)\n";
309 rapmax = cmdline.double_val(
"-etamax");
311 bool show_constituents = cmdline.present(
"-const");
312 bool massless = cmdline.present(
"-massless");
313 int nev = cmdline.int_val(
"-nev",1);
314 bool add_dense_coverage = cmdline.present(
"-dense");
315 double ghost_maxrap = cmdline.value(
"-ghost-maxrap",5.0);
316 bool all_algs = cmdline.present(
"-all-algs");
318 fj::Selector particles_sel = (cmdline.present(
"-nhardest"))
320 : fj::SelectorIdentity();
322 do_areas = cmdline.present(
"-area");
327 cmdline.value(
"-area:repeat", 1),
328 cmdline.value(
"-ghost-area", 0.01));
329 if (cmdline.present(
"-area:fj2")) ghost_spec.set_fj2_placement(
true);
330 if (cmdline.present(
"-area:explicit")) {
332 }
else if (cmdline.present(
"-area:passive")) {
334 }
else if (cmdline.present(
"-area:voronoi")) {
335 double Rfact = cmdline.value<
double>(
"-area:voronoi");
339 cmdline.present(
"-area:active");
343 bool do_bkgd = cmdline.present(
"-bkgd");
344 bool do_bkgd_csab =
false, do_bkgd_jetmedian =
false, do_bkgd_fj2 =
false;
345 bool do_bkgd_gridmedian =
false;
349 if (cmdline.present(
"-bkgd:csab")) {do_bkgd_csab =
true;}
350 else if (cmdline.present(
"-bkgd:jetmedian")) {do_bkgd_jetmedian =
true;
351 do_bkgd_fj2 = cmdline.present(
"-bkgd:fj2");
352 }
else if (cmdline.present(
"-bkgd:gridmedian")) {do_bkgd_gridmedian =
true;
354 throw fj::Error(
"with the -bkgd option, some particular background must be specified (csab or jetmedian)");
356 assert(do_areas || do_bkgd_gridmedian);
359 bool show_cones = cmdline.present(
"-cones");
363 double overlap_threshold = cmdline.double_val(
"-overlap",0.5);
364 overlap_threshold = cmdline.double_val(
"-f",overlap_threshold);
365 double seed_threshold = cmdline.double_val(
"-seed",1.0);
368 double ycut = cmdline.double_val(
"-ycut",0.08);
371 rootfile = cmdline.value<
string>(
"-root",
"");
379 vector<fj::JetDefinition> jet_defs;
380 if (all_algs || cmdline.present(
"-cam") || cmdline.present(
"-CA")) {
381 jet_defs.push_back(
fj::JetDefinition(fj::cambridge_algorithm, ktR, scheme, strategy));
383 if (all_algs || cmdline.present(
"-antikt")) {
384 jet_defs.push_back(
fj::JetDefinition(fj::antikt_algorithm, ktR, scheme, strategy));
386 if (all_algs || cmdline.present(
"-genkt")) {
388 if (cmdline.present(
"-genkt")) p = cmdline.value<
double>(
"-genkt");
390 jet_defs.push_back(
fj::JetDefinition(fj::genkt_algorithm, ktR, p, scheme, strategy));
392 if (all_algs || cmdline.present(
"-eekt")) {
395 if (all_algs || cmdline.present(
"-eegenkt")) {
397 if (cmdline.present(
"-eegenkt")) p = cmdline.value<
double>(
"-eegenkt");
399 jet_defs.push_back(
fj::JetDefinition(fj::ee_genkt_algorithm, ktR, p, scheme, strategy));
403 if (all_algs || cmdline.present(
"-midpoint")) {
404 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
405 typedef fj::CDFMidPointPlugin MPPlug;
406 double cone_area_fraction = 1.0;
407 int max_pair_size = 2;
408 int max_iterations = 100;
409 MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
410 if (cmdline.present(
"-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
411 if (cmdline.present(
"-sm-pt")) sm_scale = MPPlug::SM_pt;
412 if (cmdline.present(
"-sm-mt")) sm_scale = MPPlug::SM_mt;
413 if (cmdline.present(
"-sm-Et")) sm_scale = MPPlug::SM_Et;
416 cone_area_fraction, max_pair_size,
417 max_iterations, overlap_threshold,
419 #else // FASTJET_ENABLE_PLUGIN_CDFCONES
420 is_unavailable(
"MidPoint");
421 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
423 if (all_algs || cmdline.present(
"-pxcone")) {
424 #ifdef FASTJET_ENABLE_PLUGIN_PXCONE
425 double min_jet_energy = 5.0;
428 overlap_threshold)));
429 #else // FASTJET_ENABLE_PLUGIN_PXCONE
430 is_unavailable(
"PxCone");
431 #endif // FASTJET_ENABLE_PLUGIN_PXCONE
433 if (all_algs || cmdline.present(
"-jetclu")) {
434 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
436 ktR, overlap_threshold, seed_threshold)));
437 #else // FASTJET_ENABLE_PLUGIN_CDFCONES
438 is_unavailable(
"JetClu");
439 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
441 if (all_algs || cmdline.present(
"-siscone") || cmdline.present(
"-sisconespheri")) {
442 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
444 int npass = cmdline.value(
"-npass",0);
445 if (all_algs || cmdline.present(
"-siscone")) {
446 double sisptmin = cmdline.value(
"-sisptmin",0.0);
447 SISPlug * plugin =
new SISPlug (ktR, overlap_threshold,npass,sisptmin);
448 if (cmdline.present(
"-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
449 if (cmdline.present(
"-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
450 if (cmdline.present(
"-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
451 if (cmdline.present(
"-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
453 plugin->set_use_jet_def_recombiner(
true);
456 if (all_algs || cmdline.present(
"-sisconespheri")) {
457 double sisEmin = cmdline.value(
"-sisEmin",0.0);
460 if (cmdline.present(
"-ghost-sep")) {
461 plugin->set_ghost_separation_scale(cmdline.value<
double>(
"-ghost-sep"));
465 #else // FASTJET_ENABLE_PLUGIN_SISCONE
466 is_unavailable(
"SISCone");
467 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
469 if (all_algs || cmdline.present(
"-d0runiicone")) {
470 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
471 double min_jet_Et = 6.0;
473 #else // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
474 is_unavailable(
"D0RunIICone");
475 #endif // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
477 if (all_algs || cmdline.present(
"-trackjet")) {
478 #ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
480 #else // FASTJET_ENABLE_PLUGIN_TRACKJET
481 is_unavailable(
"TrackJet");
482 #endif // FASTJET_ENABLE_PLUGIN_TRACKJET
484 if (all_algs || cmdline.present(
"-atlascone")) {
485 #ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
487 #else // FASTJET_ENABLE_PLUGIN_ATLASCONE
488 is_unavailable(
"ATLASCone");
489 #endif // FASTJET_ENABLE_PLUGIN_ATLASCONE
491 if (all_algs || cmdline.present(
"-eecambridge")) {
492 #ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
494 #else // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
495 is_unavailable(
"EECambridge");
496 #endif // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
498 if (all_algs || cmdline.present(
"-jade")) {
499 #ifdef FASTJET_ENABLE_PLUGIN_JADE
501 #else // FASTJET_ENABLE_PLUGIN_JADE
502 is_unavailable(
"Jade");
503 #endif // FASTJET_ENABLE_PLUGIN_JADE
505 if (all_algs || cmdline.present(
"-cmsiterativecone")) {
506 #ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
508 #else // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
509 is_unavailable(
"CMSIterativeCone");
510 #endif // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
512 if (all_algs || cmdline.present(
"-d0runipre96cone")) {
513 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
514 jet_defs.push_back(
fj::JetDefinition(
new fj::D0RunIpre96ConePlugin(ktR, seed_threshold, overlap_threshold)));
515 #else // FASTJET_ENABLE_PLUGIN_D0RUNICONE
516 is_unavailable(
"D0RunIpre96Cone");
517 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
519 if (all_algs || cmdline.present(
"-d0runicone")) {
520 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
521 jet_defs.push_back(
fj::JetDefinition(
new fj::D0RunIConePlugin(ktR, seed_threshold, overlap_threshold)));
522 #else // FASTJET_ENABLE_PLUGIN_D0RUNICONE
523 is_unavailable(
"D0RunICone");
524 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
526 if (all_algs || cmdline.present(
"-gridjet")) {
527 #ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
535 double grid_ymax = 4.9999999999;
537 #else // FASTJET_ENABLE_PLUGIN_GRIDJET
538 is_unavailable(
"GridJet");
539 #endif // FASTJET_ENABLE_PLUGIN_GRIDJET
543 cmdline.present(
"-kt") ||
544 (jet_defs.size() == 0 && !found_unavailable)) {
548 string filename = cmdline.value<
string>(
"-file",
"");
551 if (!cmdline.all_options_used()) {cerr <<
552 "Error: some options were not recognized"<<endl;
555 for (
unsigned idef = 0; idef < jet_defs.size(); idef++) {
558 if (filename ==
"") istr = &cin;
559 else istr =
new ifstream(filename.c_str());
561 for (
int iev = 0; iev < nev; iev++) {
562 vector<fj::PseudoJet> jets;
563 vector<fj::PseudoJet> particles;
566 while (getline(*istr, line)) {
568 istringstream linestream(line);
569 if (line ==
"#END") {
571 if (ndone == combine) {
break;}
573 if (line.substr(0,1) ==
"#") {
continue;}
574 valarray<double> fourvec(4);
579 int ii, istat,id,m1,m2,d1,d2;
581 linestream >> ii>> istat >>
id >> m1 >> m2 >> d1 >> d2
582 >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
585 fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
586 +pow2(fourvec[2])+pow2(mass));
590 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
591 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
593 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
597 if (abs(psjet.rap()) < rapmax) {particles.push_back(psjet);}
603 if (add_dense_coverage) {
608 ghosted_area_spec.set_grid_scatter(0.5);
609 ghosted_area_spec.add_ghosts(particles);
632 particles = particles_sel(particles);
634 for (
int irepeat = 0; irepeat < repeat ; irepeat++) {
635 int nparticles = particles.size();
637 auto_ptr<fj::ClusterSequence> clust_seq;
645 if (repeatinclkt >= 0.0) {
646 vector<fj::PseudoJet> jets_local =
sorted_by_pt(clust_seq->inclusive_jets(repeatinclkt));
649 if (irepeat != 0) {
continue;}
650 cout <<
"iev "<<iev<<
": number of particles = "<< nparticles << endl;
651 cout <<
"strategy used = "<< clust_seq->strategy_string()<< endl;
653 if (do_areas && iev == 0) cout <<
"Area definition: " << area_def.description() << endl;
657 vector<fj::PseudoJet> jets_local =
sorted_by_pt(clust_seq->inclusive_jets(inclkt));
663 cout <<
"Printing "<<excln<<
" exclusive jets\n";
664 print_jets(clust_seq->exclusive_jets(excln), show_constituents);
668 cout <<
"Printing exclusive jets for d = "<<excld<<
"\n";
669 print_jets(clust_seq->exclusive_jets(excld), show_constituents);
673 cout <<
"Printing exclusive jets for ycut = "<<excly<<
"\n";
674 print_jets(clust_seq->exclusive_jets_ycut(excly), show_constituents);
678 for (
int i = nparticles-1; i >= 0; i--) {
679 printf(
"d for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_dmerge(i));
683 for (
int i = nparticles-1; i >= 0; i--) {
684 printf(
"y for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_ymerge(i));
690 if (subdcut >= 0.0) {
691 print_jets_and_sub(clust_seq->inclusive_jets(), subdcut);
696 vector<int> unique_history = clust_seq->unique_history_order();
698 vector<int> inv_unique_history(clust_seq->history().size());
699 for (
unsigned int i = 0; i < unique_history.size(); i++) {
700 inv_unique_history[unique_history[i]] = i;}
702 for (
unsigned int i = 0; i < unique_history.size(); i++) {
703 fj::ClusterSequence::history_element el =
704 clust_seq->history()[unique_history[i]];
705 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
706 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
707 printf(
"%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
712 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
718 throw fastjet::Error(
"extras object for SISCone was null (this can happen with certain area types)");
719 cout <<
"most ambiguous split (difference in squared dist) = "
721 vector<fastjet::PseudoJet> stable_cones(extras->
stable_cones());
723 for (
unsigned int i = 0; i < stable_cones.size(); i++) {
725 printf(
"%5u %15.8f %15.8f %15.8f\n",
726 i,stable_cones[i].rap(),stable_cones[i].phi(),
727 stable_cones[i].perp() );
732 vector<fj::PseudoJet> sisjets = clust_seq->inclusive_jets();
733 printf(
"\n%15s %15s %15s %12s %8s %8s\n",
"rap",
"phi",
"pt",
"user-index",
"pass",
"nconst");
734 for (
unsigned i = 0; i < sisjets.size(); i++) {
735 printf(
"%15.8f %15.8f %15.8f %12d %8d %8u\n",
736 sisjets[i].rap(), sisjets[i].phi(), sisjets[i].perp(),
737 sisjets[i].user_index(), extras->
pass(sisjets[i]),
738 (
unsigned int) clust_seq->constituents(sisjets[i]).size()
743 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
746 double rho, sigma, mean_area, empty_area, n_empty_jets;
753 }
else if (do_bkgd_jetmedian) {
755 bge.set_provide_fj2_sigma(do_bkgd_fj2);
756 bge.set_cluster_sequence(*csab);
759 mean_area = bge.mean_area();
760 empty_area = bge.empty_area();
761 n_empty_jets = bge.n_empty_jets();
763 assert(do_bkgd_gridmedian);
764 double rapmin, rapmax;
766 fj::GridMedianBackgroundEstimator bge(rapmax, 2*ktR);
767 bge.set_particles(particles);
770 mean_area = bge.mean_area();
774 cout <<
" rho = " << rho
775 <<
", sigma = " << sigma
776 <<
", mean_area = " << mean_area
777 <<
", empty_area = " << empty_area
778 <<
", n_empty_jets = " << n_empty_jets
783 cout <<
"Caught fastjet error, exiting gracefully" << endl;
790 if (jet_def.strategy()==fj::plugin_strategy){
794 if (istr != &cin)
delete istr;
805 unsigned int n_constituents = jet.
constituents().size();
806 printf(
"%15.8f %15.8f %15.8f %8u\n",
807 jet.
rap(), jet.
phi(), jet.
perp(), n_constituents);
812 void print_jets(
const vector<fj::PseudoJet> & jets_in,
bool show_constituents) {
813 vector<fj::PseudoJet> jets;
816 for (
unsigned int j = 0; j < jets.size(); j++) {
817 printf(
"%5u %15.8f %15.8f %15.8f %15.8f\n",
818 j,jets[j].px(),jets[j].py(),jets[j].pz(),jets[j].E());
819 if (show_constituents) {
820 vector<fj::PseudoJet> const_jets = jets[j].constituents();
821 for (
unsigned int k = 0; k < const_jets.size(); k++) {
822 printf(
" jet%03u %15.8f %15.8f %15.8f %15.8f\n",j,const_jets[k].px(),
823 const_jets[k].py(),const_jets[k].pz(),const_jets[k].E());
831 for (
unsigned int j = 0; j < jets.size(); j++) {
832 printf(
"%5u %15.8f %15.8f %15.8f",
833 j,jets[j].rap(),jets[j].phi(),jets[j].perp());
836 if (do_areas) printf(
" %15.8f %15.8f", jets[j].area(),
837 jets[j].area_4vector().perp());
840 if (show_constituents) {
841 vector<fj::PseudoJet> const_jets = jets[j].constituents();
842 for (
unsigned int k = 0; k < const_jets.size(); k++) {
843 printf(
" jet%03u %15.8f %15.8f %15.8f %5d\n",j,const_jets[k].rap(),
844 const_jets[k].phi(),sqrt(const_jets[k].kt2()), const_jets[k].cluster_hist_index());
851 if (rootfile !=
"") {
852 ofstream ostr(rootfile.c_str());
853 ostr <<
"# " << cmdline_p->command_line() << endl;
854 ostr <<
"# output for root" << endl;
855 assert(jets.size() > 0);
856 jets[0].validated_cs()->print_jets_for_root(jets,ostr);
865 void print_jets_and_sub (
const vector<fj::PseudoJet> & jets,
double dcut) {
871 printf(
"Printing jets and their subjets with subdcut = %10.5f\n",dcut);
872 printf(
"%5s %15s %15s %15s %15s\n",
"jet #",
"rapidity",
873 "phi",
"pt",
"n constituents");
876 enum SubType {
internal, newclust_dcut, newclust_R};
877 SubType subtype =
internal;
883 for (vector<fj::PseudoJet>::const_iterator jet = sorted_jets.begin();
884 jet != sorted_jets.end(); jet++) {
890 && jet->
perp2() < dcut)
continue;
893 printf(
"%5u ",(
unsigned int) (jet - sorted_jets.begin()));
895 vector<fj::PseudoJet> subjets;
897 if (subtype ==
internal) {
902 cout <<
" for " << ddnp1 <<
" < d < " << ddn <<
" one has " << endl;
903 }
else if (subtype == newclust_dcut) {
906 }
else if (subtype == newclust_R) {
912 cerr <<
"unrecognized subtype for subjet finding" << endl;
917 for (
unsigned int j = 0; j < subjets.size(); j++) {
918 printf(
" -sub-%02u ",j);
919 print_jet(subjets[j]);
922 if (cspoint != 0)
delete cspoint;