FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: fastjet_timing.cc 2577 2011-09-13 15:11:38Z salam $ 00003 // 00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of FastJet. 00008 // 00009 // FastJet is free software; you can redistribute it and/or modify 00010 // it under the terms of the GNU General Public License as published by 00011 // the Free Software Foundation; either version 2 of the License, or 00012 // (at your option) any later version. 00013 // 00014 // The algorithms that underlie FastJet have required considerable 00015 // development and are described in hep-ph/0512210. If you use 00016 // FastJet as part of work towards a scientific publication, please 00017 // include a citation to the FastJet paper. 00018 // 00019 // FastJet is distributed in the hope that it will be useful, 00020 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00021 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00022 // GNU General Public License for more details. 00023 // 00024 // You should have received a copy of the GNU General Public License 00025 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00026 //---------------------------------------------------------------------- 00027 //ENDHEADER 00028 00029 00030 //---------------------------------------------------------------------- 00031 /// fastjet_timing.cc: Program to help time and test the fastjet package 00032 /// 00033 /// It reads files containing multiple events in the format 00034 /// p1x p1y p1z E1 00035 /// p2x p2y p2z E2 00036 /// ... 00037 /// #END 00038 /// 00039 /// An example input file containing 10 events is included as 00040 /// data/Pythia-PtMin1000-LHC-10ev.dat 00041 /// 00042 /// Usage: 00043 /// fastjet_timing [-strategy NUMBER] [-repeat nrepeats] [-massive] \ 00044 /// [-combine nevents] [-r Rparameter] [-incl ptmin] [...] \ 00045 /// < data_file 00046 /// 00047 /// where the clustering can be repeated to aid timing and multiple 00048 /// events can be combined to get to larger multiplicities. Some options: 00049 /// 00050 /// -strategy N indicate stratgey from the enum FjStrategy (see 00051 /// FjClusterSequence.hh). 00052 /// 00053 /// -combine nev for combining multiple events from the data file in order 00054 /// to get to large multiplicities. 00055 /// 00056 /// -incl ptmin output of all inclusive jets with pt > ptmin is obtained 00057 /// with the -incl option. 00058 /// 00059 /// -excld dcut output of all exclusive jets as obtained in a clustering 00060 /// with dcut 00061 /// 00062 /// -massless read in only the 3-momenta and deduce energies assuming 00063 /// that particles are massless 00064 /// 00065 /// -write for writing out detailed clustering sequence (valuable 00066 /// for testing purposes) 00067 /// 00068 /// -unique_write writes out the sequence of dij's according to the 00069 /// "unique_history_order" (useful for verifying consistency 00070 /// between different clustering strategies). 00071 /// 00072 /// -cam switch to the inclusive Cambridge/Aachen algorithm -- 00073 /// note that the option -excld dcut provides a clustering 00074 /// up to the dcut which is the minimum squared 00075 /// distance between any pair of jets. 00076 /// 00077 #include "fastjet/PseudoJet.hh" 00078 #include "fastjet/ClusterSequence.hh" 00079 #include<iostream> 00080 #include<sstream> 00081 #include<valarray> 00082 #include<vector> 00083 #include <cstdio> 00084 #include <cstdlib> 00085 #include<cstddef> // for size_t 00086 #include "CmdLine.hh" 00087 00088 using namespace std; 00089 00090 // to avoid excessive typing, define an abbreviation for the 00091 // fastjet namespace 00092 namespace fj = fastjet; 00093 00094 inline double pow2(const double x) {return x*x;} 00095 00096 /// a program to test and time the kt algorithm as implemented in fastjet 00097 int main (int argc, char ** argv) { 00098 00099 CmdLine cmdline(argc,argv); 00100 // allow the use to specify the fj::Strategy either through the 00101 // -clever or the -strategy options (both will take numerical 00102 // values); the latter will override the former. 00103 fj::Strategy strategy = fj::Strategy(cmdline.int_val("-strategy", 00104 cmdline.int_val("-clever", fj::Best))); 00105 int repeat = cmdline.int_val("-repeat",1); 00106 int combine = cmdline.int_val("-combine",1); 00107 bool write = cmdline.present("-write"); 00108 bool unique_write = cmdline.present("-unique_write"); 00109 bool hydjet = cmdline.present("-hydjet"); 00110 double ktR = cmdline.double_val("-r",1.0); 00111 double inclkt = cmdline.double_val("-incl",-1.0); 00112 int excln = cmdline.int_val ("-excln",-1); 00113 double excld = cmdline.double_val("-excld",-1.0); 00114 double etamax = cmdline.double_val("-etamax",1.0e305); 00115 bool show_constituents = cmdline.present("-const"); 00116 bool massless = cmdline.present("-massless"); 00117 int nev = cmdline.int_val("-nev",1); 00118 bool add_dense_coverage = cmdline.present("-dense"); 00119 00120 // The following option causes the Cambridge algo to be used. 00121 // Note that currently the only output that works sensibly here is 00122 // "-incl 0" 00123 fj::JetAlgorithm jet_algorithm; 00124 if (cmdline.present("-cam")) { 00125 jet_algorithm = fj::cambridge_algorithm; 00126 } else { 00127 jet_algorithm = fj::kt_algorithm; 00128 } 00129 00130 if (!cmdline.all_options_used()) {cerr << 00131 "Error: some options were not recognized"<<endl; 00132 exit(-1);} 00133 00134 00135 for (int iev = 0; iev < nev; iev++) { 00136 vector<fj::PseudoJet> jets; 00137 string line; 00138 int ndone = 0; 00139 while (getline(cin, line)) { 00140 //cout << line<<endl; 00141 istringstream linestream(line); 00142 if (line == "#END") { 00143 ndone += 1; 00144 if (ndone == combine) {break;} 00145 } 00146 if (line.substr(0,1) == "#") {continue;} 00147 valarray<double> fourvec(4); 00148 if (hydjet) { 00149 // special reading from hydjet.txt event record (though actually 00150 // this is supposed to be a standard pythia event record, so 00151 // being able to read from it is perhaps not so bad an idea...) 00152 int ii, istat,id,m1,m2,d1,d2; 00153 double mass; 00154 linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2 00155 >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass; 00156 // current file contains mass of particle as 4th entry 00157 if (istat == 1) { 00158 fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1]) 00159 +pow2(fourvec[2])+pow2(mass)); 00160 } 00161 } else { 00162 if (massless) { 00163 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2]; 00164 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));} 00165 else { 00166 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3]; 00167 } 00168 } 00169 fj::PseudoJet psjet(fourvec); 00170 if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);} 00171 } 00172 00173 // add a fake underlying event which is very soft, uniformly distributed 00174 // in eta,phi so as to allow one to reconstruct the area that is associated 00175 // with each jet. 00176 if (add_dense_coverage) { 00177 srand(2); 00178 int nphi = 60; 00179 int neta = 100; 00180 double kt = 1e-1; 00181 for (int iphi = 0; iphi<nphi; iphi++) { 00182 for (int ieta = -neta; ieta<neta+1; ieta++) { 00183 double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX; 00184 double eta = ieta * (10.0/neta) + rand()*0.001/RAND_MAX; 00185 kt = 0.0000001*(1+rand()*0.1/RAND_MAX); 00186 double pminus = kt*exp(-eta); 00187 double pplus = kt*exp(+eta); 00188 double px = kt*sin(phi); 00189 double py = kt*cos(phi); 00190 //cout << kt<<" "<<eta<<" "<<phi<<"\n"; 00191 fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus)); 00192 jets.push_back(mom); 00193 } 00194 } 00195 } 00196 00197 fj::JetDefinition jet_def(jet_algorithm, ktR, strategy); 00198 00199 for (int irepeat = 0; irepeat < repeat ; irepeat++) { 00200 fj::ClusterSequence clust_seq(jets,jet_def,write); 00201 if (irepeat != 0) {continue;} 00202 cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl; 00203 cout << "strategy used = "<< clust_seq.strategy_string()<< endl; 00204 00205 // now provide some nice output... 00206 if (inclkt >= 0.0) { 00207 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt)); 00208 for (size_t j = 0; j < jets.size(); j++) { 00209 printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2())); 00210 if (show_constituents) { 00211 vector<fj::PseudoJet> const_jets = jets[j].constituents(); 00212 for (size_t k = 0; k < const_jets.size(); k++) { 00213 printf(" jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(), 00214 const_jets[k].phi(),sqrt(const_jets[k].kt2())); 00215 } 00216 cout << "\n\n"; 00217 } 00218 } 00219 } 00220 00221 if (excln > 0) { 00222 vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln)); 00223 00224 cout << "Printing "<<excln<<" exclusive jets\n"; 00225 for (size_t j = 0; j < jets.size(); j++) { 00226 printf("%5u %15.8f %15.8f %15.8f\n", 00227 //j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2())); 00228 j,jets[j].rap(),jets[j].phi(),jets[j].kt2()); 00229 } 00230 } 00231 00232 if (excld > 0.0) { 00233 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld)); 00234 cout << "Printing exclusive jets for d = "<<excld<<"\n"; 00235 for (size_t j = 0; j < jets.size(); j++) { 00236 printf("%5u %15.8f %15.8f %15.8f\n", 00237 j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2())); 00238 } 00239 } 00240 00241 // useful for testing that recombination sequences are unique 00242 if (unique_write) { 00243 vector<int> unique_history = clust_seq.unique_history_order(); 00244 // construct the inverse of the above mapping 00245 vector<int> inv_unique_history(clust_seq.history().size()); 00246 for (unsigned int i = 0; i < unique_history.size(); i++) { 00247 inv_unique_history[unique_history[i]] = i;} 00248 00249 for (unsigned int i = 0; i < unique_history.size(); i++) { 00250 fj::ClusterSequence::history_element el = 00251 clust_seq.history()[unique_history[i]]; 00252 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1; 00253 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2; 00254 printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2); 00255 } 00256 } 00257 00258 } // irepeat 00259 00260 } // iev 00261 }