#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 "fastjet/SISConePlugin.hh"#include "fastjet/CDFMidPointPlugin.hh"#include "fastjet/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 123 of file fastjet_timing_plugins.cc. References CmdLine::all_options_used(), fastjet::antikt_algorithm, Best, fastjet::cambridge_algorithm, CmdLine::command_line(), fastjet::ClusterSequence::constituents(), fastjet::JetDefinition::description(), fastjet::ClusterSequence::history_element::dij, CmdLine::double_val(), fastjet::ClusterSequence::exclusive_jets(), fastjet::ClusterSequence::extras(), fastjet::fastjet_version_string(), 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::ClusterSequence::print_jets_for_root(), fastjet::PseudoJet::rap(), fastjet::sorted_by_E(), fastjet::sorted_by_pt(), fastjet::sorted_by_rapidity(), fastjet::SISConeExtras::stable_cones(), fastjet::ClusterSequence::strategy_string(), twopi, fastjet::ClusterSequence::unique_history_order(), 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 << "siscone 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 ostr << "# " << cmdline.command_line() << endl;
00330 ostr << "# output for root" << endl;
00331 clust_seq.print_jets_for_root(jets,ostr);
00332 }
00333
00334 }
00335
00336 if (excln > 0) {
00337 vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
00338
00339 cout << "Printing "<<excln<<" exclusive jets\n";
00340 for (size_t j = 0; j < jets.size(); j++) {
00341 printf("%5u %15.8f %15.8f %15.8f\n",
00342 //j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00343 j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
00344 }
00345 }
00346
00347 if (excld > 0.0) {
00348 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
00349 cout << "Printing exclusive jets for d = "<<excld<<"\n";
00350 for (size_t j = 0; j < jets.size(); j++) {
00351 printf("%5u %15.8f %15.8f %15.8f\n",
00352 j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00353 }
00354 }
00355
00356
00357 // useful for testing that recombination sequences are unique
00358 if (unique_write) {
00359 vector<int> unique_history = clust_seq.unique_history_order();
00360 // construct the inverse of the above mapping
00361 vector<int> inv_unique_history(clust_seq.history().size());
00362 for (unsigned int i = 0; i < unique_history.size(); i++) {
00363 inv_unique_history[unique_history[i]] = i;}
00364
00365 for (unsigned int i = 0; i < unique_history.size(); i++) {
00366 fj::ClusterSequence::history_element el =
00367 clust_seq.history()[unique_history[i]];
00368 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00369 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00370 printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00371 }
00372 }
00373
00374
00375 #ifdef ENABLE_PLUGIN_SISCONE
00376 // provide some complementary information for SISCone
00377 if (show_cones) {
00378 const fj::SISConeExtras * extras =
00379 dynamic_cast<const fj::SISConeExtras *>(clust_seq.extras());
00380 cout << "most ambiguous split (difference in squared dist) = "
00381 << extras->most_ambiguous_split() << endl;
00382 vector<fastjet::PseudoJet> stable_cones(extras->stable_cones());
00383 stable_cones = sorted_by_rapidity(stable_cones);
00384 for (unsigned int i = 0; i < stable_cones.size(); i++) {
00385 //if (stable_cones[i].phi() < 5.0 && stable_cones[i].phi() > 4.0) {
00386 printf("%5u %15.8f %15.8f %15.8f\n",
00387 i,stable_cones[i].rap(),stable_cones[i].phi(),
00388 stable_cones[i].perp() );
00389 //}
00390 }
00391 }
00392 #endif // ENABLE_PLUGIN_SISCONE
00393 } // irepeat
00394 } // iev
00395 }
|
|
|
Definition at line 120 of file fastjet_timing_plugins.cc. 00120 {return x*x;}
|
1.4.2