#include "fastjet/PseudoJet.hh"#include "fastjet/ClusterSequence.hh"#include <iostream>#include <sstream>#include <fstream>#include <valarray>#include <vector>#include <cstdlib>#include <cstddef>#include "CmdLine.hh"#include "fastjet/config.h"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 | |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
a program to test and time the kt algorithm as implemented in fastjet
Definition at line 123 of file fastjet_timing_plugins.cc.
References CmdLine::all_options_used(), fastjet::antikt_algorithm, Best, fastjet::cambridge_algorithm, CmdLine::double_val(), fastjet::fastjet_version_string(), CmdLine::int_val(), fastjet::kt_algorithm, pow2(), CmdLine::present(), fastjet::PseudoJet::rap(), fastjet::sorted_by_E(), fastjet::sorted_by_pt(), fastjet::sorted_by_rapidity(), fastjet::ClusterSequence::strategy_string(), twopi, and CmdLine::value().
00123 { 00124 00125 CmdLine cmdline(argc,argv); 00126 // allow the use to specify the fj::Strategy either through the 00127 // -clever or the -strategy options (both will take numerical 00128 // values); the latter will override the former. 00129 fj::Strategy strategy = fj::Strategy(cmdline.int_val("-strategy", 00130 cmdline.int_val("-clever", fj::Best))); 00131 int repeat = cmdline.int_val("-repeat",1); 00132 int combine = cmdline.int_val("-combine",1); 00133 bool write = cmdline.present("-write"); 00134 bool unique_write = cmdline.present("-unique_write"); 00135 bool hydjet = cmdline.present("-hydjet"); 00136 double ktR = cmdline.double_val("-r",1.0); 00137 ktR = cmdline.double_val("-R",ktR); // allow -r and -R 00138 double inclkt = cmdline.double_val("-incl",-1.0); 00139 int excln = cmdline.int_val ("-excln",-1); 00140 double excld = cmdline.double_val("-excld",-1.0); 00141 double etamax = cmdline.double_val("-etamax",1.0e305); 00142 bool show_constituents = cmdline.present("-const"); 00143 bool massless = cmdline.present("-massless"); 00144 int nev = cmdline.int_val("-nev",1); 00145 bool add_dense_coverage = cmdline.present("-dense"); 00146 00147 bool show_cones = cmdline.present("-cones"); // only works for siscone 00148 00149 // for cone algorithms 00150 // allow -f and -overlap 00151 double overlap_threshold = cmdline.double_val("-overlap",0.5); 00152 overlap_threshold = cmdline.double_val("-f",overlap_threshold); 00153 double seed_threshold = cmdline.double_val("-seed",1.0); 00154 00155 // for printing jets to a file for reading by root 00156 string rootfile = cmdline.value<string>("-root",""); 00157 00158 // The following option causes the Cambridge algo to be used. 00159 // Note that currently the only output that works sensibly here is 00160 // "-incl 0" 00161 fj::JetDefinition jet_def; 00162 if (cmdline.present("-cam")) { 00163 jet_def = fj::JetDefinition(fj::cambridge_algorithm, ktR, strategy); 00164 } else if (cmdline.present("-antikt")) { 00165 jet_def = fj::JetDefinition(fj::antikt_algorithm, ktR, strategy); 00166 } else if (cmdline.present("-midpoint")) { 00167 #ifdef ENABLE_PLUGIN_CDFCONES 00168 typedef fj::CDFMidPointPlugin MPPlug; // for brevity 00169 double cone_area_fraction = 1.0; 00170 int max_pair_size = 2; 00171 int max_iterations = 100; 00172 MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt; 00173 if (cmdline.present("-sm-pttilde")) sm_scale = MPPlug::SM_pttilde; 00174 if (cmdline.present("-sm-pt")) sm_scale = MPPlug::SM_pt; // default 00175 if (cmdline.present("-sm-mt")) sm_scale = MPPlug::SM_mt; 00176 if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et; 00177 jet_def = fj::JetDefinition( new fj::CDFMidPointPlugin ( 00178 seed_threshold, ktR, 00179 cone_area_fraction, max_pair_size, 00180 max_iterations, overlap_threshold, 00181 sm_scale)); 00182 #else // ENABLE_PLUGIN_CDFCONES 00183 cerr << "midpoint requested, but not available for this compilation" << endl; 00184 #endif // ENABLE_PLUGIN_CDFCONES 00185 } else if (cmdline.present("-pxcone")) { 00186 #ifdef ENABLE_PLUGIN_PXCONE 00187 double min_jet_energy = 5.0; 00188 jet_def = fj::JetDefinition( new fj::PxConePlugin ( 00189 ktR, min_jet_energy, 00190 overlap_threshold)); 00191 #else // ENABLE_PLUGIN_PXCONE 00192 cerr << "pxcone requested, but not available for this compilation" << endl; 00193 exit(-1); 00194 #endif // ENABLE_PLUGIN_PXCONE 00195 } else if (cmdline.present("-jetclu")) { 00196 #ifdef ENABLE_PLUGIN_CDFCONES 00197 jet_def = fj::JetDefinition( new fj::CDFJetCluPlugin ( 00198 ktR, overlap_threshold, seed_threshold)); 00199 #else // ENABLE_PLUGIN_CDFCONES 00200 cerr << "jetclu requested, but not available for this compilation" << endl; 00201 #endif // ENABLE_PLUGIN_CDFCONES 00202 } else if (cmdline.present("-siscone")) { 00203 #ifdef ENABLE_PLUGIN_SISCONE 00204 typedef fj::SISConePlugin SISPlug; // for brevity 00205 int npass = cmdline.value("-npass",0); 00206 double sisptmin = cmdline.value("-sisptmin",0.0); 00207 SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass,sisptmin); 00208 if (cmdline.present("-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt); 00209 if (cmdline.present("-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt); 00210 if (cmdline.present("-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et); 00211 if (cmdline.present("-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde); 00212 jet_def = fj::JetDefinition(plugin); 00213 #else // ENABLE_PLUGIN_SISCONE 00214 cerr << "jetclu requested, but not available for this compilation" << endl; 00215 #endif // ENABLE_PLUGIN_SISCONE 00216 } else { 00217 jet_def = fj::JetDefinition(fj::kt_algorithm, ktR, strategy); 00218 } 00219 00220 00221 00222 if (!cmdline.all_options_used()) {cerr << 00223 "Error: some options were not recognized"<<endl; 00224 exit(-1);} 00225 00226 00227 for (int iev = 0; iev < nev; iev++) { 00228 vector<fj::PseudoJet> jets; 00229 string line; 00230 int ndone = 0; 00231 while (getline(cin, line)) { 00232 //cout << line<<endl; 00233 istringstream linestream(line); 00234 if (line == "#END") { 00235 ndone += 1; 00236 if (ndone == combine) {break;} 00237 } 00238 if (line.substr(0,1) == "#") {continue;} 00239 valarray<double> fourvec(4); 00240 if (hydjet) { 00241 // special reading from hydjet.txt event record (though actually 00242 // this is supposed to be a standard pythia event record, so 00243 // being able to read from it is perhaps not so bad an idea...) 00244 int ii, istat,id,m1,m2,d1,d2; 00245 double mass; 00246 linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2 00247 >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass; 00248 // current file contains mass of particle as 4th entry 00249 if (istat == 1) { 00250 fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1]) 00251 +pow2(fourvec[2])+pow2(mass)); 00252 } 00253 } else { 00254 if (massless) { 00255 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2]; 00256 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));} 00257 else { 00258 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3]; 00259 } 00260 } 00261 fj::PseudoJet psjet(fourvec); 00262 if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);} 00263 } 00264 00265 // add a fake underlying event which is very soft, uniformly distributed 00266 // in eta,phi so as to allow one to reconstruct the area that is associated 00267 // with each jet. 00268 if (add_dense_coverage) { 00269 srand(2); 00270 int nphi = 60; 00271 int neta = 100; 00272 double kt = 1e-1; 00273 for (int iphi = 0; iphi<nphi; iphi++) { 00274 for (int ieta = -neta; ieta<neta+1; ieta++) { 00275 double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX; 00276 double eta = ieta * (10.0/neta) + rand()*0.001/RAND_MAX; 00277 kt = 0.0000001*(1+rand()*0.1/RAND_MAX); 00278 double pminus = kt*exp(-eta); 00279 double pplus = kt*exp(+eta); 00280 double px = kt*sin(phi); 00281 double py = kt*cos(phi); 00282 //cout << kt<<" "<<eta<<" "<<phi<<"\n"; 00283 fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus)); 00284 jets.push_back(mom); 00285 } 00286 } 00287 } 00288 00289 for (int irepeat = 0; irepeat < repeat ; irepeat++) { 00290 fj::ClusterSequence clust_seq(jets,jet_def,write); 00291 if (irepeat != 0) {continue;} 00292 cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl; 00293 cout << "strategy used = "<< clust_seq.strategy_string()<< endl; 00294 cout << "Algorithm: " << jet_def.description() << " (" << fj::fastjet_version_string() << ")" << endl; 00295 00296 // now provide some nice output... 00297 if (inclkt >= 0.0) { 00298 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt)); 00299 for (size_t j = 0; j < jets.size(); j++) { 00300 //printf("%5u %15.8f %15.8f %15.8e\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2())); 00301 printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2())); 00302 if (show_constituents) { 00303 vector<fj::PseudoJet> const_jets = clust_seq.constituents(jets[j]); 00304 for (size_t k = 0; k < const_jets.size(); k++) { 00305 printf(" jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(), 00306 const_jets[k].phi(),sqrt(const_jets[k].kt2())); 00307 } 00308 cout << "\n\n"; 00309 } 00310 // family tests 00311 //fj::PseudoJet parent1, parent2; 00312 //bool has_parents = clust_seq.has_parents(jets[j],parent1,parent2); 00313 //if (has_parents) { 00314 // cout << "parent pt's: " << parent1.perp() << " " << parent2.perp() <<endl; 00315 // fj::PseudoJet child_from1, child_from2; 00316 // clust_seq.has_child(parent1, child_from1); 00317 // clust_seq.has_child(parent2, child_from2); 00318 // cout << "parents' children pt: " << child_from1.perp() << " " << child_from2.perp() << endl; 00319 // fj::PseudoJet partner_from1, partner_from2; 00320 // clust_seq.has_partner(parent1, partner_from1); 00321 // clust_seq.has_partner(parent2, partner_from2); 00322 // cout << "parents' partners' pt: " << partner_from1.perp() << " " << partner_from2.perp() << endl; 00323 //} else { 00324 // cout << "has no parents" << endl; 00325 //} 00326 } 00327 if (rootfile != "") { 00328 ofstream ostr(rootfile.c_str()); 00329 clust_seq.print_jets_for_root(jets,ostr); 00330 } 00331 00332 } 00333 00334 if (excln > 0) { 00335 vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln)); 00336 00337 cout << "Printing "<<excln<<" exclusive jets\n"; 00338 for (size_t j = 0; j < jets.size(); j++) { 00339 printf("%5u %15.8f %15.8f %15.8f\n", 00340 //j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2())); 00341 j,jets[j].rap(),jets[j].phi(),jets[j].kt2()); 00342 } 00343 } 00344 00345 if (excld > 0.0) { 00346 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld)); 00347 cout << "Printing exclusive jets for d = "<<excld<<"\n"; 00348 for (size_t j = 0; j < jets.size(); j++) { 00349 printf("%5u %15.8f %15.8f %15.8f\n", 00350 j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2())); 00351 } 00352 } 00353 00354 00355 // useful for testing that recombination sequences are unique 00356 if (unique_write) { 00357 vector<int> unique_history = clust_seq.unique_history_order(); 00358 // construct the inverse of the above mapping 00359 vector<int> inv_unique_history(clust_seq.history().size()); 00360 for (unsigned int i = 0; i < unique_history.size(); i++) { 00361 inv_unique_history[unique_history[i]] = i;} 00362 00363 for (unsigned int i = 0; i < unique_history.size(); i++) { 00364 fj::ClusterSequence::history_element el = 00365 clust_seq.history()[unique_history[i]]; 00366 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1; 00367 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2; 00368 printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2); 00369 } 00370 } 00371 00372 00373 // provide some complementary information for SISCone 00374 if (show_cones) { 00375 const fj::SISConeExtras * extras = 00376 dynamic_cast<const fj::SISConeExtras *>(clust_seq.extras()); 00377 cout << "most ambiguous split (difference in squared dist) = " 00378 << extras->most_ambiguous_split() << endl; 00379 vector<fastjet::PseudoJet> stable_cones(extras->stable_cones()); 00380 stable_cones = sorted_by_rapidity(stable_cones); 00381 for (unsigned int i = 0; i < stable_cones.size(); i++) { 00382 //if (stable_cones[i].phi() < 5.0 && stable_cones[i].phi() > 4.0) { 00383 printf("%5u %15.8f %15.8f %15.8f\n", 00384 i,stable_cones[i].rap(),stable_cones[i].phi(), 00385 stable_cones[i].perp() ); 00386 //} 00387 } 00388 } 00389 } // irepeat 00390 00391 } // iev 00392 }
| double pow2 | ( | const double | x | ) | [inline] |
1.5.2