00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00089 #include "fastjet/PseudoJet.hh"
00090 #include "fastjet/ClusterSequence.hh"
00091 #include<iostream>
00092 #include<sstream>
00093 #include<fstream>
00094 #include<valarray>
00095 #include<vector>
00096 #include <cstdlib>
00097 #include<cstddef>
00098 #include "CmdLine.hh"
00099
00100
00101 #include "fastjet/config.h"
00102
00103 #ifdef ENABLE_PLUGIN_SISCONE
00104 #include "fastjet/SISConePlugin.hh"
00105 #endif
00106 #ifdef ENABLE_PLUGIN_CDFCONES
00107 #include "fastjet/CDFMidPointPlugin.hh"
00108 #include "fastjet/CDFJetCluPlugin.hh"
00109 #endif
00110 #ifdef ENABLE_PLUGIN_PXCONE
00111 #include "fastjet/PxConePlugin.hh"
00112 #endif
00113
00114 using namespace std;
00115
00116
00117
00118 namespace fj = fastjet;
00119
00120 inline double pow2(const double x) {return x*x;}
00121
00123 int main (int argc, char ** argv) {
00124
00125 CmdLine cmdline(argc,argv);
00126
00127
00128
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);
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");
00148
00149
00150
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
00156 string rootfile = cmdline.value<string>("-root","");
00157
00158
00159
00160
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;
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;
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;
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 << "jetclu 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
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
00242
00243
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
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
00266
00267
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
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
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
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
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326 }
00327 if (rootfile != "") {
00328 ofstream ostr(rootfile.c_str());
00329 clust_seq.print_jets_for_root(jets,ostr);
00330 }
00331
00332 }
00333
00334 if (excln > 0) {
00335 vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
00336
00337 cout << "Printing "<<excln<<" exclusive jets\n";
00338 for (size_t j = 0; j < jets.size(); j++) {
00339 printf("%5u %15.8f %15.8f %15.8f\n",
00340
00341 j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
00342 }
00343 }
00344
00345 if (excld > 0.0) {
00346 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
00347 cout << "Printing exclusive jets for d = "<<excld<<"\n";
00348 for (size_t j = 0; j < jets.size(); j++) {
00349 printf("%5u %15.8f %15.8f %15.8f\n",
00350 j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00351 }
00352 }
00353
00354
00355
00356 if (unique_write) {
00357 vector<int> unique_history = clust_seq.unique_history_order();
00358
00359 vector<int> inv_unique_history(clust_seq.history().size());
00360 for (unsigned int i = 0; i < unique_history.size(); i++) {
00361 inv_unique_history[unique_history[i]] = i;}
00362
00363 for (unsigned int i = 0; i < unique_history.size(); i++) {
00364 fj::ClusterSequence::history_element el =
00365 clust_seq.history()[unique_history[i]];
00366 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00367 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00368 printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00369 }
00370 }
00371
00372
00373
00374 if (show_cones) {
00375 const fj::SISConeExtras * extras =
00376 dynamic_cast<const fj::SISConeExtras *>(clust_seq.extras());
00377 cout << "most ambiguous split (difference in squared dist) = "
00378 << extras->most_ambiguous_split() << endl;
00379 vector<fastjet::PseudoJet> stable_cones(extras->stable_cones());
00380 stable_cones = sorted_by_rapidity(stable_cones);
00381 for (unsigned int i = 0; i < stable_cones.size(); i++) {
00382
00383 printf("%5u %15.8f %15.8f %15.8f\n",
00384 i,stable_cones[i].rap(),stable_cones[i].phi(),
00385 stable_cones[i].perp() );
00386
00387 }
00388 }
00389 }
00390
00391 }
00392 }