00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00160
00161
00162 #include "fastjet/PseudoJet.hh"
00163 #include "fastjet/ClusterSequence.hh"
00164 #include "fastjet/GhostedAreaSpec.hh"
00165 #include<iostream>
00166 #include<sstream>
00167 #include<fstream>
00168 #include<valarray>
00169 #include<vector>
00170 #include <cstdlib>
00171 #include<cstddef>
00172 #include "CmdLine.hh"
00173
00174
00175 #include "fastjet/config.h"
00176
00177
00178 #ifdef ENABLE_PLUGIN_SISCONE
00179 #include "fastjet/SISConePlugin.hh"
00180 #include "fastjet/SISConeSphericalPlugin.hh"
00181 #endif
00182 #ifdef ENABLE_PLUGIN_CDFCONES
00183 #include "fastjet/CDFMidPointPlugin.hh"
00184 #include "fastjet/CDFJetCluPlugin.hh"
00185 #endif
00186 #ifdef ENABLE_PLUGIN_PXCONE
00187 #include "fastjet/PxConePlugin.hh"
00188 #endif
00189 #ifdef ENABLE_PLUGIN_D0RUNIICONE
00190 #include "fastjet/D0RunIIConePlugin.hh"
00191 #endif
00192 #ifdef ENABLE_PLUGIN_TRACKJET
00193 #include "fastjet/TrackJetPlugin.hh"
00194 #endif
00195 #ifdef ENABLE_PLUGIN_ATLASCONE
00196 #include "fastjet/ATLASConePlugin.hh"
00197 #endif
00198 #ifdef ENABLE_PLUGIN_EECAMBRIDGE
00199 #include "fastjet/EECambridgePlugin.hh"
00200 #endif
00201 #ifdef ENABLE_PLUGIN_JADE
00202 #include "fastjet/JadePlugin.hh"
00203 #endif
00204 #ifdef ENABLE_PLUGIN_CMSITERATIVECONE
00205 #include "fastjet/CMSIterativeConePlugin.hh"
00206 #endif
00207
00208
00209 using namespace std;
00210
00211
00212
00213 namespace fj = fastjet;
00214
00215 inline double pow2(const double x) {return x*x;}
00216
00217
00218 void print_jets_and_sub (fj::ClusterSequence & clust_seq,
00219 const vector<fj::PseudoJet> & jets, double dcut);
00220
00221 string rootfile;
00222 CmdLine * cmdline_p;
00223
00226 bool ee_print = false;
00227 void print_jets(const vector<fj::PseudoJet> & jets, const fj::ClusterSequence & cs, bool show_const = false);
00228
00229 void is_unavailable(const string & algname) {
00230 cerr << algname << " requested, but not available for this compilation";
00231 exit(0);
00232 }
00233
00234
00237 int main (int argc, char ** argv) {
00238
00239 CmdLine cmdline(argc,argv);
00240 cmdline_p = &cmdline;
00241
00242
00243
00244 fj::Strategy strategy = fj::Strategy(cmdline.int_val("-strategy",
00245 cmdline.int_val("-clever", fj::Best)));
00246 int repeat = cmdline.int_val("-repeat",1);
00247 int combine = cmdline.int_val("-combine",1);
00248 bool write = cmdline.present("-write");
00249 bool unique_write = cmdline.present("-unique_write");
00250 bool hydjet = cmdline.present("-hydjet");
00251 double ktR = cmdline.double_val("-r",1.0);
00252 ktR = cmdline.double_val("-R",ktR);
00253 double inclkt = cmdline.double_val("-incl",-1.0);
00254 int excln = cmdline.int_val ("-excln",-1);
00255 double excld = cmdline.double_val("-excld",-1.0);
00256 double excly = cmdline.double_val("-excly",-1.0);
00257 ee_print = cmdline.present("-ee-print");
00258 bool get_all_dij = cmdline.present("-get-all-dij");
00259 bool get_all_yij = cmdline.present("-get-all-yij");
00260 double subdcut = cmdline.double_val("-subdcut",-1.0);
00261 double etamax = cmdline.double_val("-etamax",1.0e305);
00262 bool show_constituents = cmdline.present("-const");
00263 bool massless = cmdline.present("-massless");
00264 int nev = cmdline.int_val("-nev",1);
00265 bool add_dense_coverage = cmdline.present("-dense");
00266 double ghost_maxrap = cmdline.value("-ghost-maxrap",5.0);
00267
00268 bool show_cones = cmdline.present("-cones");
00269
00270
00271
00272 double overlap_threshold = cmdline.double_val("-overlap",0.5);
00273 overlap_threshold = cmdline.double_val("-f",overlap_threshold);
00274 double seed_threshold = cmdline.double_val("-seed",1.0);
00275
00276
00277 double ycut = cmdline.double_val("-ycut",0.08);
00278
00279
00280 rootfile = cmdline.value<string>("-root","");
00281
00282
00283 fj::RecombinationScheme scheme = fj::E_scheme;
00284
00285
00286
00287
00288 fj::JetDefinition jet_def;
00289 if (cmdline.present("-cam") || cmdline.present("-CA")) {
00290 jet_def = fj::JetDefinition(fj::cambridge_algorithm, ktR, scheme, strategy);
00291 } else if (cmdline.present("-antikt")) {
00292 jet_def = fj::JetDefinition(fj::antikt_algorithm, ktR, scheme, strategy);
00293 } else if (cmdline.present("-genkt")) {
00294 double p = cmdline.value<double>("-genkt");
00295 jet_def = fj::JetDefinition(fj::genkt_algorithm, ktR, p, scheme, strategy);
00296 } else if (cmdline.present("-eekt")) {
00297 jet_def = fj::JetDefinition(fj::ee_kt_algorithm);
00298 } else if (cmdline.present("-eegenkt")) {
00299 double p = cmdline.value<double>("-eegenkt");
00300 jet_def = fj::JetDefinition(fj::ee_genkt_algorithm, ktR, p, scheme, strategy);
00301
00302
00303 } else if (cmdline.present("-midpoint")) {
00304 #ifdef ENABLE_PLUGIN_CDFCONES
00305 typedef fj::CDFMidPointPlugin MPPlug;
00306 double cone_area_fraction = 1.0;
00307 int max_pair_size = 2;
00308 int max_iterations = 100;
00309 MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
00310 if (cmdline.present("-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
00311 if (cmdline.present("-sm-pt")) sm_scale = MPPlug::SM_pt;
00312 if (cmdline.present("-sm-mt")) sm_scale = MPPlug::SM_mt;
00313 if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et;
00314 jet_def = fj::JetDefinition( new fj::CDFMidPointPlugin (
00315 seed_threshold, ktR,
00316 cone_area_fraction, max_pair_size,
00317 max_iterations, overlap_threshold,
00318 sm_scale));
00319 #else // ENABLE_PLUGIN_CDFCONES
00320 is_unavailable("midpoint");
00321 #endif // ENABLE_PLUGIN_CDFCONES
00322 } else if (cmdline.present("-pxcone")) {
00323 #ifdef ENABLE_PLUGIN_PXCONE
00324 double min_jet_energy = 5.0;
00325 jet_def = fj::JetDefinition( new fj::PxConePlugin (
00326 ktR, min_jet_energy,
00327 overlap_threshold));
00328 #else // ENABLE_PLUGIN_PXCONE
00329 is_unavailable("pxcone");
00330 #endif // ENABLE_PLUGIN_PXCONE
00331 } else if (cmdline.present("-jetclu")) {
00332 #ifdef ENABLE_PLUGIN_CDFCONES
00333 jet_def = fj::JetDefinition( new fj::CDFJetCluPlugin (
00334 ktR, overlap_threshold, seed_threshold));
00335 #else // ENABLE_PLUGIN_CDFCONES
00336 is_unavailable("pxcone");
00337 #endif // ENABLE_PLUGIN_CDFCONES
00338 } else if (cmdline.present("-siscone") || cmdline.present("-sisconespheri")) {
00339 #ifdef ENABLE_PLUGIN_SISCONE
00340 typedef fj::SISConePlugin SISPlug;
00341 int npass = cmdline.value("-npass",0);
00342 if (cmdline.present("-siscone")) {
00343 double sisptmin = cmdline.value("-sisptmin",0.0);
00344 SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass,sisptmin);
00345 if (cmdline.present("-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
00346 if (cmdline.present("-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
00347 if (cmdline.present("-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
00348 if (cmdline.present("-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
00349
00350 plugin->set_use_jet_def_recombiner(true);
00351 jet_def = fj::JetDefinition(plugin);
00352 } else {
00353 double sisEmin = cmdline.value("-sisEmin",0.0);
00354 fj::SISConeSphericalPlugin * plugin =
00355 new fj::SISConeSphericalPlugin(ktR, overlap_threshold,npass,sisEmin);
00356 if (cmdline.present("-ghost-sep")) {
00357 plugin->set_ghost_separation_scale(cmdline.value<double>("-ghost-sep"));
00358 }
00359 jet_def = fj::JetDefinition(plugin);
00360 }
00361 #else // ENABLE_PLUGIN_SISCONE
00362 is_unavailable("siscone");
00363 #endif // ENABLE_PLUGIN_SISCONE
00364 } else if (cmdline.present("-d0runiicone")) {
00365 #ifdef ENABLE_PLUGIN_D0RUNIICONE
00366 double min_jet_Et = 6.0;
00367 jet_def = fj::JetDefinition(new fj::D0RunIIConePlugin(ktR,min_jet_Et));
00368 #else // ENABLE_PLUGIN_D0RUNIICONE
00369 is_unavailable("D0RunIICone");
00370 #endif // ENABLE_PLUGIN_D0RUNIICONE
00371 } else if (cmdline.present("-trackjet")) {
00372 #ifdef ENABLE_PLUGIN_TRACKJET
00373 jet_def = fj::JetDefinition(new fj::TrackJetPlugin(ktR));
00374 #else // ENABLE_PLUGIN_TRACKJET
00375 is_unavailable("TrackJet");
00376 #endif // ENABLE_PLUGIN_TRACKJET
00377 } else if (cmdline.present("-atlascone")) {
00378 #ifdef ENABLE_PLUGIN_ATLASCONE
00379 jet_def = fj::JetDefinition(new fj::ATLASConePlugin(ktR));
00380 #else // ENABLE_PLUGIN_ATLASCONE
00381 is_unavailable("ATLASCone");
00382 #endif // ENABLE_PLUGIN_ATLASCONE
00383 } else if (cmdline.present("-eecambridge")) {
00384 #ifdef ENABLE_PLUGIN_EECAMBRIDGE
00385 jet_def = fj::JetDefinition(new fj::EECambridgePlugin(ycut));
00386 #else // ENABLE_PLUGIN_EECAMBRIDGE
00387 is_unavailable("EECambridge");
00388 #endif // ENABLE_PLUGIN_EECAMBRIDGE
00389 } else if (cmdline.present("-jade")) {
00390 #ifdef ENABLE_PLUGIN_JADE
00391 jet_def = fj::JetDefinition(new fj::JadePlugin());
00392 #else // ENABLE_PLUGIN_JADE
00393 is_unavailable("Jade");
00394 #endif // ENABLE_PLUGIN_JADE
00395 } else if (cmdline.present("-cmsiterativecone")) {
00396 #ifdef ENABLE_PLUGIN_CMSITERATIVECONE
00397 jet_def = fj::JetDefinition(new fj::CMSIterativeConePlugin(ktR,seed_threshold));
00398 #else // ENABLE_PLUGIN_CMSITERATIVECONE
00399 is_unavailable("CMSIterativeCone");
00400 #endif // ENABLE_PLUGIN_CMSITERATIVECONE
00401
00402 } else {
00403 cmdline.present("-kt");
00404 jet_def = fj::JetDefinition(fj::kt_algorithm, ktR, strategy);
00405 }
00406
00407
00408
00409 if (!cmdline.all_options_used()) {cerr <<
00410 "Error: some options were not recognized"<<endl;
00411 exit(-1);}
00412
00413
00414 for (int iev = 0; iev < nev; iev++) {
00415 vector<fj::PseudoJet> jets;
00416 string line;
00417 int ndone = 0;
00418 while (getline(cin, line)) {
00419
00420 istringstream linestream(line);
00421 if (line == "#END") {
00422 ndone += 1;
00423 if (ndone == combine) {break;}
00424 }
00425 if (line.substr(0,1) == "#") {continue;}
00426 valarray<double> fourvec(4);
00427 if (hydjet) {
00428
00429
00430
00431 int ii, istat,id,m1,m2,d1,d2;
00432 double mass;
00433 linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00434 >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00435
00436 if (istat == 1) {
00437 fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00438 +pow2(fourvec[2])+pow2(mass));
00439 }
00440 } else {
00441 if (massless) {
00442 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00443 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00444 else {
00445 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00446 }
00447 }
00448 fj::PseudoJet psjet(fourvec);
00449 if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
00450 }
00451
00452
00453
00454
00455 if (add_dense_coverage) {
00456 fj::GhostedAreaSpec ghosted_area_spec(ghost_maxrap);
00457
00458
00459
00460 ghosted_area_spec.set_grid_scatter(0.5);
00461 ghosted_area_spec.add_ghosts(jets);
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481 }
00482
00483 for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00484 int nparticles = jets.size();
00485 try {
00486 fj::ClusterSequence clust_seq(jets,jet_def,write);
00487 if (irepeat != 0) {continue;}
00488 cout << "iev "<<iev<< ": number of particles = "<< nparticles << endl;
00489 cout << "strategy used = "<< clust_seq.strategy_string()<< endl;
00490 cout << "Algorithm: " << jet_def.description() << " (" << fj::fastjet_version_string() << ")" << endl;
00491
00492
00493 if (inclkt >= 0.0) {
00494 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt));
00495 print_jets(jets, clust_seq, show_constituents);
00496
00497 }
00498
00499 if (excln > 0) {
00500 cout << "Printing "<<excln<<" exclusive jets\n";
00501 print_jets(clust_seq.exclusive_jets(excln), clust_seq, show_constituents);
00502 }
00503
00504 if (excld > 0.0) {
00505 cout << "Printing exclusive jets for d = "<<excld<<"\n";
00506 print_jets(clust_seq.exclusive_jets(excld), clust_seq, show_constituents);
00507 }
00508
00509 if (excly > 0.0) {
00510 cout << "Printing exclusive jets for ycut = "<<excly<<"\n";
00511 print_jets(clust_seq.exclusive_jets_ycut(excly), clust_seq, show_constituents);
00512 }
00513
00514 if (get_all_dij) {
00515 for (int i = nparticles-1; i >= 0; i--) {
00516 printf("d for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq.exclusive_dmerge(i));
00517 }
00518 }
00519 if (get_all_yij) {
00520 for (int i = nparticles-1; i >= 0; i--) {
00521 printf("y for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq.exclusive_ymerge(i));
00522 }
00523 }
00524
00525
00526
00527 if (subdcut >= 0.0) {
00528 print_jets_and_sub(clust_seq, clust_seq.inclusive_jets(), subdcut);
00529 }
00530
00531
00532 if (unique_write) {
00533 vector<int> unique_history = clust_seq.unique_history_order();
00534
00535 vector<int> inv_unique_history(clust_seq.history().size());
00536 for (unsigned int i = 0; i < unique_history.size(); i++) {
00537 inv_unique_history[unique_history[i]] = i;}
00538
00539 for (unsigned int i = 0; i < unique_history.size(); i++) {
00540 fj::ClusterSequence::history_element el =
00541 clust_seq.history()[unique_history[i]];
00542 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00543 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00544 printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00545 }
00546 }
00547
00548
00549 #ifdef ENABLE_PLUGIN_SISCONE
00550
00551 if (show_cones) {
00552 const fj::SISConeExtras * extras =
00553 dynamic_cast<const fj::SISConeExtras *>(clust_seq.extras());
00554 cout << "most ambiguous split (difference in squared dist) = "
00555 << extras->most_ambiguous_split() << endl;
00556 vector<fastjet::PseudoJet> stable_cones(extras->stable_cones());
00557 stable_cones = sorted_by_rapidity(stable_cones);
00558 for (unsigned int i = 0; i < stable_cones.size(); i++) {
00559
00560 printf("%5u %15.8f %15.8f %15.8f\n",
00561 i,stable_cones[i].rap(),stable_cones[i].phi(),
00562 stable_cones[i].perp() );
00563
00564 }
00565
00566
00567 vector<fj::PseudoJet> sisjets = clust_seq.inclusive_jets();
00568 printf("\n%15s %15s %15s %12s %8s %8s\n","rap","phi","pt","user-index","pass","nconst");
00569 for (unsigned i = 0; i < sisjets.size(); i++) {
00570 printf("%15.8f %15.8f %15.8f %12d %8d %8d\n",
00571 sisjets[i].rap(), sisjets[i].phi(), sisjets[i].perp(),
00572 sisjets[i].user_index(), extras->pass(sisjets[i]),
00573 clust_seq.constituents(sisjets[i]).size()
00574 );
00575
00576 }
00577 }
00578 #endif // ENABLE_PLUGIN_SISCONE
00579 }
00580 catch (fastjet::Error fjerr) {
00581 cout << "Caught fastjet error, exiting gracefully" << endl;
00582 exit(0);
00583 }
00584
00585 }
00586 }
00587
00588
00589 if (jet_def.strategy()==fj::plugin_strategy){
00590 delete jet_def.plugin();
00591 }
00592 }
00593
00594
00595
00596
00597
00599 void print_jet (const fj::ClusterSequence & clust_seq,
00600 const fj::PseudoJet & jet) {
00601 int n_constituents = clust_seq.constituents(jet).size();
00602 printf("%15.8f %15.8f %15.8f %8u\n",
00603 jet.rap(), jet.phi(), jet.perp(), n_constituents);
00604 }
00605
00606
00607
00608 void print_jets(const vector<fj::PseudoJet> & jets_in, const fj::ClusterSequence & cs, bool show_constituents) {
00609 vector<fj::PseudoJet> jets;
00610 if (ee_print) {
00611 jets = sorted_by_E(jets_in);
00612 for (size_t j = 0; j < jets.size(); j++) {
00613 printf("%5u %15.8f %15.8f %15.8f %15.8f\n",
00614 j,jets[j].px(),jets[j].py(),jets[j].pz(),jets[j].E());
00615 if (show_constituents) {
00616 vector<fj::PseudoJet> const_jets = cs.constituents(jets[j]);
00617 for (size_t k = 0; k < const_jets.size(); k++) {
00618 printf(" jet%03u %15.8f %15.8f %15.8f %15.8f\n",j,const_jets[k].px(),
00619 const_jets[k].py(),const_jets[k].pz(),const_jets[k].E());
00620 }
00621 cout << "\n\n";
00622 }
00623
00624 }
00625 } else {
00626 jets = sorted_by_pt(jets_in);
00627 for (size_t j = 0; j < jets.size(); j++) {
00628 printf("%5u %15.8f %15.8f %15.8f\n",
00629 j,jets[j].rap(),jets[j].phi(),jets[j].perp());
00630
00631 if (show_constituents) {
00632 vector<fj::PseudoJet> const_jets = cs.constituents(jets[j]);
00633 for (size_t k = 0; k < const_jets.size(); k++) {
00634 printf(" jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
00635 const_jets[k].phi(),sqrt(const_jets[k].kt2()));
00636 }
00637 cout << "\n\n";
00638 }
00639 }
00640 }
00641
00642 if (rootfile != "") {
00643 ofstream ostr(rootfile.c_str());
00644 ostr << "# " << cmdline_p->command_line() << endl;
00645 ostr << "# output for root" << endl;
00646 cs.print_jets_for_root(jets,ostr);
00647 }
00648
00649 }
00650
00651
00652
00655 void print_jets_and_sub (fj::ClusterSequence & clust_seq,
00656 const vector<fj::PseudoJet> & jets, double dcut) {
00657
00658
00659 vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);
00660
00661
00662 printf("Printing jets and their subjets with subdcut = %10.5f\n",dcut);
00663 printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity",
00664 "phi", "pt", "n constituents");
00665
00666
00667 enum SubType {internal, newclust_dcut, newclust_R};
00668 SubType subtype = internal;
00669
00670
00671
00672 for (unsigned int i = 0; i < sorted_jets.size(); i++) {
00673
00674
00675 if (clust_seq.jet_def().jet_algorithm() == fj::kt_algorithm
00676 && sorted_jets[i].perp2() < dcut) continue;
00677
00678 printf("%5u ",i);
00679 print_jet(clust_seq, sorted_jets[i]);
00680 vector<fj::PseudoJet> subjets;
00681 fj::ClusterSequence * cspoint;
00682 if (subtype == internal) {
00683 cspoint = &clust_seq;
00684 subjets = clust_seq.exclusive_subjets(sorted_jets[i], dcut);
00685
00686 double ddnp1 = clust_seq.exclusive_subdmerge_max(sorted_jets[i], subjets.size());
00687 double ddn = clust_seq.exclusive_subdmerge_max(sorted_jets[i], subjets.size()-1);
00688 cout << " for " << ddnp1 << " < d < " << ddn << " one has " << endl;
00689
00690 } else if (subtype == newclust_dcut) {
00691 cspoint = new fj::ClusterSequence(clust_seq.constituents(sorted_jets[i]),
00692 clust_seq.jet_def());
00693 subjets = cspoint->exclusive_jets(dcut);
00694
00695 } else if (subtype == newclust_R) {
00696 assert(clust_seq.jet_def().jet_algorithm() == fj::cambridge_algorithm);
00697 fj::JetDefinition subjd(clust_seq.jet_def().jet_algorithm(),
00698 clust_seq.jet_def().R()*sqrt(dcut));
00699 cspoint = new fj::ClusterSequence(clust_seq.constituents(sorted_jets[i]),
00700 subjd);
00701 subjets = cspoint->inclusive_jets();
00702 } else {
00703 cerr << "unrecognized subtype for subjet finding" << endl;
00704 exit(-1);
00705 }
00706
00707 subjets = sorted_by_pt(subjets);
00708 for (unsigned int j = 0; j < subjets.size(); j++) {
00709 printf(" -sub-%02u ",j);
00710 print_jet(*cspoint, subjets[j]);
00711 }
00712
00713 if (cspoint != &clust_seq) delete cspoint;
00714
00715
00716
00717
00718
00719
00720
00721
00722 }
00723
00724 }
00725