fastjet_timing_plugins.cc File Reference

#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/GhostedAreaSpec.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"
Include dependency graph for fastjet_timing_plugins.cc:

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 ee_print = false
 sort and pretty print jets, with exact behaviour depending on whether ee_print is true or not

Function Documentation

void is_unavailable ( const string &  algname  ) 

Definition at line 229 of file fastjet_timing_plugins.cc.

Referenced by main().

00229                                             {
00230   cerr << algname << " requested, but not available for this compilation";
00231   exit(0);
00232 }

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 237 of file fastjet_timing_plugins.cc.

References CmdLine::all_options_used(), antikt_algorithm, Best, cambridge_algorithm, CmdLine::double_val(), E_scheme, ee_genkt_algorithm, ee_kt_algorithm, ee_print, fastjet_version_string(), genkt_algorithm, CmdLine::int_val(), is_unavailable(), kt_algorithm, fastjet::SISConeBaseExtras::most_ambiguous_split(), fastjet::SISConeBaseExtras::pass(), fastjet::d0::inline_maths::phi(), plugin_strategy, pow2(), CmdLine::present(), print_jets(), print_jets_and_sub(), rootfile, fastjet::SISConeBasePlugin::set_ghost_separation_scale(), sorted_by_pt(), sorted_by_rapidity(), fastjet::SISConeBaseExtras::stable_cones(), and CmdLine::value().

00237                                   {
00238 
00239   CmdLine cmdline(argc,argv);
00240   cmdline_p = &cmdline;
00241   // allow the use to specify the fj::Strategy either through the
00242   // -clever or the -strategy options (both will take numerical
00243   // values); the latter will override the former.
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); // allow -r and -R
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"); // only works for siscone
00269 
00270   // for cone algorithms
00271   // allow -f and -overlap
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   // for ee algorithms, allow to specify ycut
00277   double ycut = cmdline.double_val("-ycut",0.08);
00278 
00279   // for printing jets to a file for reading by root
00280   rootfile = cmdline.value<string>("-root","");
00281 
00282   // out default scheme is the E_scheme
00283   fj::RecombinationScheme scheme = fj::E_scheme;
00284 
00285   // The following option causes the Cambridge algo to be used.
00286   // Note that currently the only output that works sensibly here is
00287   // "-incl 0"
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 // checking if one asks to run a plugin (don't delete this line)
00303   } else if (cmdline.present("-midpoint")) {
00304 #ifdef ENABLE_PLUGIN_CDFCONES
00305     typedef fj::CDFMidPointPlugin MPPlug; // for brevity
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; // default
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; // for brevity
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       // cause it to use the jet-definition's own recombiner
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; // was 8 GeV in earlier work
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 // end of checking if one asks to run a plugin (don't delete this line)
00402   } else {
00403     cmdline.present("-kt"); // kt is default, but allow user to specify it too [and ignore return value!]
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       //cout << line<<endl;
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       // special reading from hydjet.txt event record (though actually
00429       // this is supposed to be a standard pythia event record, so
00430       // being able to read from it is perhaps not so bad an idea...)
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       // current file contains mass of particle as 4th entry
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   // add a fake underlying event which is very soft, uniformly distributed
00453   // in eta,phi so as to allow one to reconstruct the area that is associated
00454   // with each jet.
00455   if (add_dense_coverage) {
00456     fj::GhostedAreaSpec ghosted_area_spec(ghost_maxrap);
00457     //fj::GhostedAreaSpec ghosted_area_spec(-2.0,4.0); // asymmetric range
00458     // for plots, reduce the scatter default of 1, to avoid "holes"
00459     // in the subsequent calorimeter view
00460     ghosted_area_spec.set_grid_scatter(0.5); 
00461     ghosted_area_spec.add_ghosts(jets);
00462     //----- old code ------------------
00463     // srand(2);
00464     // int nphi = 60;
00465     // int neta = 100;
00466     // double kt = 1e-1;
00467     // for (int iphi = 0; iphi<nphi; iphi++) {
00468     //   for (int ieta = -neta; ieta<neta+1; ieta++) {
00469     //  double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00470     //  double eta = ieta * (10.0/neta)  + rand()*0.001/RAND_MAX;
00471     //  kt = 1e-20*(1+rand()*0.1/RAND_MAX);
00472     //  double pminus = kt*exp(-eta);
00473     //  double pplus  = kt*exp(+eta);
00474     //  double px = kt*sin(phi);
00475     //  double py = kt*cos(phi);
00476     //  //cout << kt<<" "<<eta<<" "<<phi<<"\n";
00477     //  fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00478     //  jets.push_back(mom);
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     // now provide some nice output...
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     // have the option of printing out the subjets (at scale dcut) of
00526     // each inclusive jet
00527     if (subdcut >= 0.0) {
00528       print_jets_and_sub(clust_seq, clust_seq.inclusive_jets(), subdcut);
00529     }
00530     
00531     // useful for testing that recombination sequences are unique
00532     if (unique_write) {
00533       vector<int> unique_history = clust_seq.unique_history_order();
00534       // construct the inverse of the above mapping
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     // provide some complementary information for SISCone 
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       //if (stable_cones[i].phi() < 5.0 && stable_cones[i].phi() > 4.0) {
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       // also show passes for jets
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   } // try
00580   catch (fastjet::Error fjerr) {
00581     cout << "Caught fastjet error, exiting gracefully" << endl;
00582     exit(0);
00583   }
00584 
00585   } // irepeat
00586   } // iev
00587 
00588   // if we've instantiated a plugin, delete it
00589   if (jet_def.strategy()==fj::plugin_strategy){
00590     delete jet_def.plugin();
00591   }
00592 }

double pow2 ( const double  x  )  [inline]

Definition at line 215 of file fastjet_timing_plugins.cc.

00215 {return x*x;}

void print_jet ( const fj::ClusterSequence &  clust_seq,
const fj::PseudoJet &  jet 
)

print a single jet

Definition at line 599 of file fastjet_timing_plugins.cc.

00600                                          {
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 }

void print_jets ( const vector< fj::PseudoJet > &  jets,
const fj::ClusterSequence &  cs,
bool  show_const = false 
)

Definition at line 608 of file fastjet_timing_plugins.cc.

References CmdLine::command_line(), ee_print, fastjet::d0::inline_maths::phi(), rootfile, sorted_by_E(), and sorted_by_pt().

00608                                                                                                            {
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 }

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 655 of file fastjet_timing_plugins.cc.

References cambridge_algorithm, kt_algorithm, print_jet(), and sorted_by_pt().

00656                                                                         {
00657 
00658   // sort jets into increasing pt
00659   vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);  
00660 
00661   // label the columns
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   // have various kinds of subjet finding, to test consistency among them
00667   enum SubType {internal, newclust_dcut, newclust_R};
00668   SubType subtype = internal;
00669   //SubType subtype = newclust_dcut;
00670 
00671   // print out the details for each jet
00672   for (unsigned int i = 0; i < sorted_jets.size(); i++) {
00673     // if jet pt^2 < dcut with kt alg, then some methods of
00674     // getting subjets will return nothing -- so skip the jet
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       //subjets = clust_seq.exclusive_subjets(sorted_jets[i], 5);
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       //subjets = clust_seq.exclusive_subjets(sorted_jets[i], dd*1.0000001);
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       //subjets = cspoint->exclusive_jets(int(min(5U,cspoint->n_particles())));
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     //fj::ClusterSequence subseq(clust_seq.constituents(sorted_jets[i]),
00716     //                          fj::JetDefinition(fj::cambridge_algorithm, 0.4));
00717     //vector<fj::PseudoJet> subjets = sorted_by_pt(subseq.inclusive_jets());
00718     //for (unsigned int j = 0; j < subjets.size(); j++) {
00719     //  printf("    -sub-%02u ",j);
00720     //  print_jet(subseq, subjets[j]);
00721     //}
00722   }
00723 
00724 }

Variable Documentation

CmdLine* cmdline_p

Definition at line 222 of file fastjet_timing_plugins.cc.

bool ee_print = false

sort and pretty print jets, with exact behaviour depending on whether ee_print is true or not

Definition at line 226 of file fastjet_timing_plugins.cc.

Referenced by main(), and print_jets().

string rootfile

Definition at line 221 of file fastjet_timing_plugins.cc.

Referenced by main(), and print_jets().

Generated on 26 Feb 2010 for fastjet by  doxygen 1.6.1