#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 |
|
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 }
|
|
Definition at line 111 of file fastjet_timing_plugins.cc. 00111 {return x*x;}
|