Main Page | Namespace List | Class Hierarchy | Class List | Directories | File List | Namespace Members | Class Members | File Members

fastjet_timing_plugins.cc File Reference

#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include <iostream>
#include <sstream>
#include <valarray>
#include <vector>
#include <cstdlib>
#include <cstddef>
#include "CmdLine.hh"
#include "PxConePlugin.hh"
#include "SISConePlugin.hh"
#include "CDFMidPointPlugin.hh"
#include "CDFJetCluPlugin.hh"

Include dependency graph for fastjet_timing_plugins.cc:

Go to the source code of this file.

Functions

double pow2 (const double x)
int main (int argc, char **argv)
 a program to test and time the kt algorithm as implemented in fastjet


Function Documentation

int main int  argc,
char **  argv
 

a program to test and time the kt algorithm as implemented in fastjet

Definition at line 114 of file fastjet_timing_plugins.cc.

References CmdLine::all_options_used(), Best, fastjet::cambridge_algorithm, fastjet::ClusterSequence::constituents(), fastjet::JetDefinition::description(), fastjet::ClusterSequence::history_element::dij, CmdLine::double_val(), fastjet::ClusterSequence::exclusive_jets(), fastjet::ClusterSequence::extras(), fastjet::ClusterSequence::history(), fastjet::ClusterSequence::inclusive_jets(), CmdLine::int_val(), fastjet::kt_algorithm, fastjet::SISConeExtras::most_ambiguous_split(), fastjet::ClusterSequence::history_element::parent1, fastjet::ClusterSequence::history_element::parent2, pow2(), CmdLine::present(), fastjet::PseudoJet::rap(), fastjet::sorted_by_E(), fastjet::sorted_by_pt(), fastjet::ClusterSequence::strategy_string(), twopi, fastjet::ClusterSequence::unique_history_order(), and CmdLine::value().

00114                                   {
00115 
00116   CmdLine cmdline(argc,argv);
00117   // allow the use to specify the fj::Strategy either through the
00118   // -clever or the -strategy options (both will take numerical
00119   // values); the latter will override the former.
00120   fj::Strategy  strategy  = fj::Strategy(cmdline.int_val("-strategy",
00121                                         cmdline.int_val("-clever", fj::Best)));
00122   int  repeat  = cmdline.int_val("-repeat",1);
00123   int  combine = cmdline.int_val("-combine",1);
00124   bool write   = cmdline.present("-write");
00125   bool unique_write = cmdline.present("-unique_write");
00126   bool hydjet  = cmdline.present("-hydjet");
00127   double ktR   = cmdline.double_val("-r",1.0);
00128   ktR   = cmdline.double_val("-R",ktR); // allow -r and -R
00129   double inclkt = cmdline.double_val("-incl",-1.0);
00130   int    excln  = cmdline.int_val   ("-excln",-1);
00131   double excld  = cmdline.double_val("-excld",-1.0);
00132   double etamax = cmdline.double_val("-etamax",1.0e305);
00133   bool   show_constituents = cmdline.present("-const");
00134   bool   massless = cmdline.present("-massless");
00135   int  nev     = cmdline.int_val("-nev",1);
00136   bool add_dense_coverage = cmdline.present("-dense");
00137 
00138   bool show_cones = cmdline.present("-cones"); // only works for siscone
00139 
00140   // for cone algorithms
00141   // allow -f and -overlap
00142   double overlap_threshold = cmdline.double_val("-overlap",0.5);
00143   overlap_threshold = cmdline.double_val("-f",overlap_threshold); 
00144   double seed_threshold = cmdline.double_val("-seed",1.0);
00145 
00146   // The following option causes the Cambridge algo to be used.
00147   // Note that currently the only output that works sensibly here is
00148   // "-incl 0"
00149   fj::JetDefinition jet_def;
00150   if (cmdline.present("-cam")) {
00151     jet_def = fj::JetDefinition(fj::cambridge_algorithm, ktR, strategy);
00152   } else if (cmdline.present("-midpoint")) {
00153     typedef fj::CDFMidPointPlugin MPPlug; // for brevity
00154     double cone_area_fraction = 1.0;
00155     int    max_pair_size = 2;
00156     int    max_iterations = 100;
00157     MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
00158     if (cmdline.present("-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
00159     if (cmdline.present("-sm-pt")) sm_scale = MPPlug::SM_pt; // default
00160     if (cmdline.present("-sm-mt")) sm_scale = MPPlug::SM_mt;
00161     if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et;
00162     jet_def = fj::JetDefinition( new fj::CDFMidPointPlugin (
00163                                       seed_threshold, ktR, 
00164                                       cone_area_fraction, max_pair_size,
00165                                       max_iterations, overlap_threshold,
00166                                       sm_scale));
00167   } else if (cmdline.present("-pxcone")) {
00168     double min_jet_energy = 5.0;
00169     jet_def = fj::JetDefinition( new fj::PxConePlugin (
00170                                       ktR, min_jet_energy,
00171                                       overlap_threshold));
00172   } else if (cmdline.present("-jetclu")) {
00173     double seed_threshold = 1.0;
00174     jet_def = fj::JetDefinition( new fj::CDFJetCluPlugin (
00175                                       ktR, overlap_threshold, seed_thshold));
00176   } else if (cmdline.present("-siscone")) {
00177     typedef fj::SISConePlugin SISPlug; // for brevity
00178     int npass = cmdline.value("-npass",1);
00179     SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass);
00180     if (cmdline.present("-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
00181     if (cmdline.present("-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
00182     if (cmdline.present("-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
00183     if (cmdline.present("-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
00184     jet_def = fj::JetDefinition(plugin);
00185   } else {
00186     jet_def = fj::JetDefinition(fj::kt_algorithm, ktR, strategy);
00187   }
00188 
00189 
00190 
00191   if (!cmdline.all_options_used()) {cerr << 
00192       "Error: some options were not recognized"<<endl; 
00193     exit(-1);}
00194 
00195 
00196   for (int iev = 0; iev < nev; iev++) {
00197   vector<fj::PseudoJet> jets;
00198   string line;
00199   int  ndone = 0;
00200   while (getline(cin, line)) {
00201       //cout << line<<endl;
00202     istringstream linestream(line);
00203     if (line == "#END") {
00204       ndone += 1;
00205       if (ndone == combine) {break;}
00206     }
00207     if (line.substr(0,1) == "#") {continue;}
00208     valarray<double> fourvec(4);
00209     if (hydjet) {
00210       // special reading from hydjet.txt event record (though actually
00211       // this is supposed to be a standard pythia event record, so
00212       // being able to read from it is perhaps not so bad an idea...)
00213       int ii, istat,id,m1,m2,d1,d2;
00214       double mass;
00215       linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00216                  >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00217       // current file contains mass of particle as 4th entry
00218       if (istat == 1) {
00219         fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00220                           +pow2(fourvec[2])+pow2(mass));
00221       }
00222     } else {
00223       if (massless) {
00224         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00225         fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00226       else {
00227         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00228       }
00229     }
00230     fj::PseudoJet psjet(fourvec);
00231     if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
00232   }
00233 
00234   // add a fake underlying event which is very soft, uniformly distributed
00235   // in eta,phi so as to allow one to reconstruct the area that is associated
00236   // with each jet.
00237   if (add_dense_coverage) {
00238     srand(2);
00239     int nphi = 60;
00240     int neta = 100;
00241     double kt = 1e-1;
00242     for (int iphi = 0; iphi<nphi; iphi++) {
00243       for (int ieta = -neta; ieta<neta+1; ieta++) {
00244         double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00245         double eta = ieta * (10.0/neta)  + rand()*0.001/RAND_MAX;
00246         kt = 0.0000001*(1+rand()*0.1/RAND_MAX);
00247         double pminus = kt*exp(-eta);
00248         double pplus  = kt*exp(+eta);
00249         double px = kt*sin(phi);
00250         double py = kt*cos(phi);
00251         //cout << kt<<" "<<eta<<" "<<phi<<"\n";
00252         fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00253         jets.push_back(mom);
00254       }
00255     }
00256   }
00257   
00258   for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00259     fj::ClusterSequence clust_seq(jets,jet_def,write);
00260     if (irepeat != 0) {continue;}
00261     cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl;
00262     cout << "strategy used =  "<< clust_seq.strategy_string()<< endl;
00263     cout << "Algorithm: " << jet_def.description() << endl;
00264 
00265     // now provide some nice output...
00266     if (inclkt >= 0.0) {
00267       vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt));
00268       for (size_t j = 0; j < jets.size(); j++) {
00269         //printf("%5u %15.8f %15.8f %15.8e\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00270         printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00271         if (show_constituents) {
00272           vector<fj::PseudoJet> const_jets = clust_seq.constituents(jets[j]);
00273           for (size_t k = 0; k < const_jets.size(); k++) {
00274             printf("        jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
00275                    const_jets[k].phi(),sqrt(const_jets[k].kt2()));
00276           }
00277           cout << "\n\n";
00278         }
00279       }
00280     }
00281 
00282     if (excln > 0) {
00283       vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
00284  
00285       cout << "Printing "<<excln<<" exclusive jets\n";
00286       for (size_t j = 0; j < jets.size(); j++) {
00287         printf("%5u %15.8f %15.8f %15.8f\n",
00288                //j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00289                j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
00290       }
00291     }
00292 
00293     if (excld > 0.0) {
00294       vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
00295       cout << "Printing exclusive jets for d = "<<excld<<"\n";
00296       for (size_t j = 0; j < jets.size(); j++) {
00297         printf("%5u %15.8f %15.8f %15.8f\n",
00298                j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00299       }
00300     }
00301     
00302     // useful for testing that recombination sequences are unique
00303     if (unique_write) {
00304       vector<int> unique_history = clust_seq.unique_history_order();
00305       // construct the inverse of the above mapping
00306       vector<int> inv_unique_history(clust_seq.history().size());
00307       for (unsigned int i = 0; i < unique_history.size(); i++) {
00308         inv_unique_history[unique_history[i]] = i;}
00309 
00310       for (unsigned int i = 0; i < unique_history.size(); i++) {
00311         fj::ClusterSequence::history_element el = 
00312           clust_seq.history()[unique_history[i]];
00313         int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00314         int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00315         printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00316       }
00317     }
00318 
00319 
00320     // provide some complementary information for SISCone 
00321     if (show_cones) {
00322       const fj::SISConeExtras * extras = 
00323         dynamic_cast<const fj::SISConeExtras *>(clust_seq.extras());
00324       cout << "most ambiguous split (difference in squared dist) = "
00325            << extras->most_ambiguous_split() << endl;
00326     }
00327   } // irepeat
00328 
00329   } // iev
00330 }

double pow2 const double  x  )  [inline]
 

Definition at line 111 of file fastjet_timing_plugins.cc.

00111 {return x*x;}


Generated on Mon Apr 2 20:57:52 2007 for fastjet by  doxygen 1.4.2