#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 |
int main | ( | int | argc, | |
char ** | argv | |||
) |
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, CmdLine::double_val(), CmdLine::int_val(), fastjet::kt_algorithm, pow2(), CmdLine::present(), fastjet::PseudoJet::rap(), fastjet::sorted_by_E(), fastjet::sorted_by_pt(), fastjet::ClusterSequence::strategy_string(), and twopi.
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::JetAlgorithm jet_algorithm; 00125 if (cmdline.present("-cam")) { 00126 jet_algorithm = fj::cambridge_algorithm; 00127 } else { 00128 jet_algorithm = 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_algorithm, 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 }
double pow2 | ( | const double | x | ) | [inline] |