#include "fastjet/PseudoJet.hh"#include "fastjet/ClusterSequence.hh"#include <iostream>#include <sstream>#include <valarray>#include <vector>#include <cstdlib>#include <cstddef>#include "CmdLine.hh"Include dependency graph for fastjet_timing.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 98 of file fastjet_timing.cc. References CmdLine::all_options_used(), Best, fastjet::cambridge_algorithm, fastjet::ClusterSequence::constituents(), fastjet::ClusterSequence::history_element::dij, CmdLine::double_val(), fastjet::ClusterSequence::exclusive_jets(), fastjet::ClusterSequence::history(), fastjet::ClusterSequence::inclusive_jets(), CmdLine::int_val(), fastjet::kt_algorithm, 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, and fastjet::ClusterSequence::unique_history_order(). 00098 {
00099
00100 CmdLine cmdline(argc,argv);
00101 // allow the use to specify the fj::Strategy either through the
00102 // -clever or the -strategy options (both will take numerical
00103 // values); the latter will override the former.
00104 fj::Strategy strategy = fj::Strategy(cmdline.int_val("-strategy",
00105 cmdline.int_val("-clever", fj::Best)));
00106 int repeat = cmdline.int_val("-repeat",1);
00107 int combine = cmdline.int_val("-combine",1);
00108 bool write = cmdline.present("-write");
00109 bool unique_write = cmdline.present("-unique_write");
00110 bool hydjet = cmdline.present("-hydjet");
00111 double ktR = cmdline.double_val("-r",1.0);
00112 double inclkt = cmdline.double_val("-incl",-1.0);
00113 int excln = cmdline.int_val ("-excln",-1);
00114 double excld = cmdline.double_val("-excld",-1.0);
00115 double etamax = cmdline.double_val("-etamax",1.0e305);
00116 bool show_constituents = cmdline.present("-const");
00117 bool massless = cmdline.present("-massless");
00118 int nev = cmdline.int_val("-nev",1);
00119 bool add_dense_coverage = cmdline.present("-dense");
00120
00121 // The following option causes the Cambridge algo to be used.
00122 // Note that currently the only output that works sensibly here is
00123 // "-incl 0"
00124 fj::JetFinder jet_finder;
00125 if (cmdline.present("-cam")) {
00126 jet_finder = fj::cambridge_algorithm;
00127 } else {
00128 jet_finder = fj::kt_algorithm;
00129 }
00130
00131 if (!cmdline.all_options_used()) {cerr <<
00132 "Error: some options were not recognized"<<endl;
00133 exit(-1);}
00134
00135
00136 for (int iev = 0; iev < nev; iev++) {
00137 vector<fj::PseudoJet> jets;
00138 string line;
00139 int ndone = 0;
00140 while (getline(cin, line)) {
00141 //cout << line<<endl;
00142 istringstream linestream(line);
00143 if (line == "#END") {
00144 ndone += 1;
00145 if (ndone == combine) {break;}
00146 }
00147 if (line.substr(0,1) == "#") {continue;}
00148 valarray<double> fourvec(4);
00149 if (hydjet) {
00150 // special reading from hydjet.txt event record (though actually
00151 // this is supposed to be a standard pythia event record, so
00152 // being able to read from it is perhaps not so bad an idea...)
00153 int ii, istat,id,m1,m2,d1,d2;
00154 double mass;
00155 linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00156 >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00157 // current file contains mass of particle as 4th entry
00158 if (istat == 1) {
00159 fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00160 +pow2(fourvec[2])+pow2(mass));
00161 }
00162 } else {
00163 if (massless) {
00164 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00165 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00166 else {
00167 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00168 }
00169 }
00170 fj::PseudoJet psjet(fourvec);
00171 if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
00172 }
00173
00174 // add a fake underlying event which is very soft, uniformly distributed
00175 // in eta,phi so as to allow one to reconstruct the area that is associated
00176 // with each jet.
00177 if (add_dense_coverage) {
00178 srand(2);
00179 int nphi = 60;
00180 int neta = 100;
00181 double kt = 1e-1;
00182 for (int iphi = 0; iphi<nphi; iphi++) {
00183 for (int ieta = -neta; ieta<neta+1; ieta++) {
00184 double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00185 double eta = ieta * (10.0/neta) + rand()*0.001/RAND_MAX;
00186 kt = 0.0000001*(1+rand()*0.1/RAND_MAX);
00187 double pminus = kt*exp(-eta);
00188 double pplus = kt*exp(+eta);
00189 double px = kt*sin(phi);
00190 double py = kt*cos(phi);
00191 //cout << kt<<" "<<eta<<" "<<phi<<"\n";
00192 fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00193 jets.push_back(mom);
00194 }
00195 }
00196 }
00197
00198 fj::JetDefinition jet_def(jet_finder, ktR, strategy);
00199
00200 for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00201 fj::ClusterSequence clust_seq(jets,jet_def,write);
00202 if (irepeat != 0) {continue;}
00203 cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl;
00204 cout << "strategy used = "<< clust_seq.strategy_string()<< endl;
00205
00206 // now provide some nice output...
00207 if (inclkt >= 0.0) {
00208 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt));
00209 for (size_t j = 0; j < jets.size(); j++) {
00210 printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00211 if (show_constituents) {
00212 vector<fj::PseudoJet> const_jets = clust_seq.constituents(jets[j]);
00213 for (size_t k = 0; k < const_jets.size(); k++) {
00214 printf(" jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
00215 const_jets[k].phi(),sqrt(const_jets[k].kt2()));
00216 }
00217 cout << "\n\n";
00218 }
00219 }
00220 }
00221
00222 if (excln > 0) {
00223 vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
00224
00225 cout << "Printing "<<excln<<" exclusive jets\n";
00226 for (size_t j = 0; j < jets.size(); j++) {
00227 printf("%5u %15.8f %15.8f %15.8f\n",
00228 //j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00229 j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
00230 }
00231 }
00232
00233 if (excld > 0.0) {
00234 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
00235 cout << "Printing exclusive jets for d = "<<excld<<"\n";
00236 for (size_t j = 0; j < jets.size(); j++) {
00237 printf("%5u %15.8f %15.8f %15.8f\n",
00238 j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00239 }
00240 }
00241
00242 // useful for testing that recombination sequences are unique
00243 if (unique_write) {
00244 vector<int> unique_history = clust_seq.unique_history_order();
00245 // construct the inverse of the above mapping
00246 vector<int> inv_unique_history(clust_seq.history().size());
00247 for (unsigned int i = 0; i < unique_history.size(); i++) {
00248 inv_unique_history[unique_history[i]] = i;}
00249
00250 for (unsigned int i = 0; i < unique_history.size(); i++) {
00251 fj::ClusterSequence::history_element el =
00252 clust_seq.history()[unique_history[i]];
00253 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00254 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00255 printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00256 }
00257 }
00258
00259 } // irepeat
00260
00261 } // iev
00262 }
|
|
|
Definition at line 95 of file fastjet_timing.cc. Referenced by main(). 00095 {return x*x;}
|
1.4.2