#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::SISConePlugin::set_split_merge_on_transverse_mass(), 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-mt")) sm_scale = MPPlug::SM_mt;
00159 if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et;
00160 jet_def = fj::JetDefinition( new fj::CDFMidPointPlugin (
00161 seed_threshold, ktR,
00162 cone_area_fraction, max_pair_size,
00163 max_iterations, overlap_threshold,
00164 sm_scale));
00165 } else if (cmdline.present("-pxcone")) {
00166 double min_jet_energy = 5.0;
00167 jet_def = fj::JetDefinition( new fj::PxConePlugin (
00168 ktR, min_jet_energy,
00169 overlap_threshold));
00170 } else if (cmdline.present("-jetclu")) {
00171 double seed_threshold = 1.0;
00172 jet_def = fj::JetDefinition( new fj::CDFJetCluPlugin (
00173 ktR, overlap_threshold, seed_threshold));
00174 } else if (cmdline.present("-siscone")) {
00175 int npass = cmdline.value("-npass",1);
00176 fj::SISConePlugin * plugin = new fj::SISConePlugin (ktR, overlap_threshold,npass);
00177 if (cmdline.present("-nomt")) plugin->set_split_merge_on_transverse_mass(false);
00178 jet_def = fj::JetDefinition(plugin);
00179 } else {
00180 jet_def = fj::JetDefinition(fj::kt_algorithm, ktR, strategy);
00181 }
00182
00183
00184
00185 if (!cmdline.all_options_used()) {cerr <<
00186 "Error: some options were not recognized"<<endl;
00187 exit(-1);}
00188
00189
00190 for (int iev = 0; iev < nev; iev++) {
00191 vector<fj::PseudoJet> jets;
00192 string line;
00193 int ndone = 0;
00194 while (getline(cin, line)) {
00195 //cout << line<<endl;
00196 istringstream linestream(line);
00197 if (line == "#END") {
00198 ndone += 1;
00199 if (ndone == combine) {break;}
00200 }
00201 if (line.substr(0,1) == "#") {continue;}
00202 valarray<double> fourvec(4);
00203 if (hydjet) {
00204 // special reading from hydjet.txt event record (though actually
00205 // this is supposed to be a standard pythia event record, so
00206 // being able to read from it is perhaps not so bad an idea...)
00207 int ii, istat,id,m1,m2,d1,d2;
00208 double mass;
00209 linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00210 >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00211 // current file contains mass of particle as 4th entry
00212 if (istat == 1) {
00213 fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00214 +pow2(fourvec[2])+pow2(mass));
00215 }
00216 } else {
00217 if (massless) {
00218 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00219 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00220 else {
00221 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00222 }
00223 }
00224 fj::PseudoJet psjet(fourvec);
00225 if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
00226 }
00227
00228 // add a fake underlying event which is very soft, uniformly distributed
00229 // in eta,phi so as to allow one to reconstruct the area that is associated
00230 // with each jet.
00231 if (add_dense_coverage) {
00232 srand(2);
00233 int nphi = 60;
00234 int neta = 100;
00235 double kt = 1e-1;
00236 for (int iphi = 0; iphi<nphi; iphi++) {
00237 for (int ieta = -neta; ieta<neta+1; ieta++) {
00238 double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00239 double eta = ieta * (10.0/neta) + rand()*0.001/RAND_MAX;
00240 kt = 0.0000001*(1+rand()*0.1/RAND_MAX);
00241 double pminus = kt*exp(-eta);
00242 double pplus = kt*exp(+eta);
00243 double px = kt*sin(phi);
00244 double py = kt*cos(phi);
00245 //cout << kt<<" "<<eta<<" "<<phi<<"\n";
00246 fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00247 jets.push_back(mom);
00248 }
00249 }
00250 }
00251
00252 for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00253 fj::ClusterSequence clust_seq(jets,jet_def,write);
00254 if (irepeat != 0) {continue;}
00255 cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl;
00256 cout << "strategy used = "<< clust_seq.strategy_string()<< endl;
00257 cout << "Algorithm: " << jet_def.description() << endl;
00258
00259 // now provide some nice output...
00260 if (inclkt >= 0.0) {
00261 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt));
00262 for (size_t j = 0; j < jets.size(); j++) {
00263 //printf("%5u %15.8f %15.8f %15.8e\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00264 printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00265 if (show_constituents) {
00266 vector<fj::PseudoJet> const_jets = clust_seq.constituents(jets[j]);
00267 for (size_t k = 0; k < const_jets.size(); k++) {
00268 printf(" jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
00269 const_jets[k].phi(),sqrt(const_jets[k].kt2()));
00270 }
00271 cout << "\n\n";
00272 }
00273 }
00274 }
00275
00276 if (excln > 0) {
00277 vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
00278
00279 cout << "Printing "<<excln<<" exclusive jets\n";
00280 for (size_t j = 0; j < jets.size(); j++) {
00281 printf("%5u %15.8f %15.8f %15.8f\n",
00282 //j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00283 j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
00284 }
00285 }
00286
00287 if (excld > 0.0) {
00288 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
00289 cout << "Printing exclusive jets for d = "<<excld<<"\n";
00290 for (size_t j = 0; j < jets.size(); j++) {
00291 printf("%5u %15.8f %15.8f %15.8f\n",
00292 j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00293 }
00294 }
00295
00296 // useful for testing that recombination sequences are unique
00297 if (unique_write) {
00298 vector<int> unique_history = clust_seq.unique_history_order();
00299 // construct the inverse of the above mapping
00300 vector<int> inv_unique_history(clust_seq.history().size());
00301 for (unsigned int i = 0; i < unique_history.size(); i++) {
00302 inv_unique_history[unique_history[i]] = i;}
00303
00304 for (unsigned int i = 0; i < unique_history.size(); i++) {
00305 fj::ClusterSequence::history_element el =
00306 clust_seq.history()[unique_history[i]];
00307 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00308 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00309 printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00310 }
00311 }
00312
00313
00314 // provide some complementary information for SISCone
00315 if (show_cones) {
00316 const fj::SISConeExtras * extras =
00317 dynamic_cast<const fj::SISConeExtras *>(clust_seq.extras());
00318 cout << "most ambiguous split (difference in squared dist) = "
00319 << extras->most_ambiguous_split() << endl;
00320 }
00321 } // irepeat
00322
00323 } // iev
00324 }
|
|
|
Definition at line 111 of file fastjet_timing_plugins.cc. 00111 {return x*x;}
|
1.4.2