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<valarray>
00094 #include<vector>
00095 #include <cstdlib>
00096 #include<cstddef>
00097 #include "CmdLine.hh"
00098
00099
00100 #include "PxConePlugin.hh"
00101 #include "SISConePlugin.hh"
00102 #include "CDFMidPointPlugin.hh"
00103 #include "CDFJetCluPlugin.hh"
00104
00105 using namespace std;
00106
00107
00108
00109 namespace fj = fastjet;
00110
00111 inline double pow2(const double x) {return x*x;}
00112
00114 int main (int argc, char ** argv) {
00115
00116 CmdLine cmdline(argc,argv);
00117
00118
00119
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);
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");
00139
00140
00141
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
00147
00148
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;
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-pttilde")) sm_scale = MPPlug::SM_pttilde;
00159 if (cmdline.present("-sm-pt")) sm_scale = MPPlug::SM_pt;
00160 if (cmdline.present("-sm-mt")) sm_scale = MPPlug::SM_mt;
00161 if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et;
00162 jet_def = fj::JetDefinition( new fj::CDFMidPointPlugin (
00163 seed_threshold, ktR,
00164 cone_area_fraction, max_pair_size,
00165 max_iterations, overlap_threshold,
00166 sm_scale));
00167 } else if (cmdline.present("-pxcone")) {
00168 double min_jet_energy = 5.0;
00169 jet_def = fj::JetDefinition( new fj::PxConePlugin (
00170 ktR, min_jet_energy,
00171 overlap_threshold));
00172 } else if (cmdline.present("-jetclu")) {
00173 double seed_threshold = 1.0;
00174 jet_def = fj::JetDefinition( new fj::CDFJetCluPlugin (
00175 ktR, overlap_threshold, seed_thshold));
00176 } else if (cmdline.present("-siscone")) {
00177 typedef fj::SISConePlugin SISPlug;
00178 int npass = cmdline.value("-npass",1);
00179 SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass);
00180 if (cmdline.present("-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
00181 if (cmdline.present("-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
00182 if (cmdline.present("-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
00183 if (cmdline.present("-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
00184 jet_def = fj::JetDefinition(plugin);
00185 } else {
00186 jet_def = fj::JetDefinition(fj::kt_algorithm, ktR, strategy);
00187 }
00188
00189
00190
00191 if (!cmdline.all_options_used()) {cerr <<
00192 "Error: some options were not recognized"<<endl;
00193 exit(-1);}
00194
00195
00196 for (int iev = 0; iev < nev; iev++) {
00197 vector<fj::PseudoJet> jets;
00198 string line;
00199 int ndone = 0;
00200 while (getline(cin, line)) {
00201
00202 istringstream linestream(line);
00203 if (line == "#END") {
00204 ndone += 1;
00205 if (ndone == combine) {break;}
00206 }
00207 if (line.substr(0,1) == "#") {continue;}
00208 valarray<double> fourvec(4);
00209 if (hydjet) {
00210
00211
00212
00213 int ii, istat,id,m1,m2,d1,d2;
00214 double mass;
00215 linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00216 >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00217
00218 if (istat == 1) {
00219 fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00220 +pow2(fourvec[2])+pow2(mass));
00221 }
00222 } else {
00223 if (massless) {
00224 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00225 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00226 else {
00227 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00228 }
00229 }
00230 fj::PseudoJet psjet(fourvec);
00231 if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
00232 }
00233
00234
00235
00236
00237 if (add_dense_coverage) {
00238 srand(2);
00239 int nphi = 60;
00240 int neta = 100;
00241 double kt = 1e-1;
00242 for (int iphi = 0; iphi<nphi; iphi++) {
00243 for (int ieta = -neta; ieta<neta+1; ieta++) {
00244 double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00245 double eta = ieta * (10.0/neta) + rand()*0.001/RAND_MAX;
00246 kt = 0.0000001*(1+rand()*0.1/RAND_MAX);
00247 double pminus = kt*exp(-eta);
00248 double pplus = kt*exp(+eta);
00249 double px = kt*sin(phi);
00250 double py = kt*cos(phi);
00251
00252 fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00253 jets.push_back(mom);
00254 }
00255 }
00256 }
00257
00258 for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00259 fj::ClusterSequence clust_seq(jets,jet_def,write);
00260 if (irepeat != 0) {continue;}
00261 cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl;
00262 cout << "strategy used = "<< clust_seq.strategy_string()<< endl;
00263 cout << "Algorithm: " << jet_def.description() << endl;
00264
00265
00266 if (inclkt >= 0.0) {
00267 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt));
00268 for (size_t j = 0; j < jets.size(); j++) {
00269
00270 printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00271 if (show_constituents) {
00272 vector<fj::PseudoJet> const_jets = clust_seq.constituents(jets[j]);
00273 for (size_t k = 0; k < const_jets.size(); k++) {
00274 printf(" jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
00275 const_jets[k].phi(),sqrt(const_jets[k].kt2()));
00276 }
00277 cout << "\n\n";
00278 }
00279 }
00280 }
00281
00282 if (excln > 0) {
00283 vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
00284
00285 cout << "Printing "<<excln<<" exclusive jets\n";
00286 for (size_t j = 0; j < jets.size(); j++) {
00287 printf("%5u %15.8f %15.8f %15.8f\n",
00288
00289 j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
00290 }
00291 }
00292
00293 if (excld > 0.0) {
00294 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
00295 cout << "Printing exclusive jets for d = "<<excld<<"\n";
00296 for (size_t j = 0; j < jets.size(); j++) {
00297 printf("%5u %15.8f %15.8f %15.8f\n",
00298 j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00299 }
00300 }
00301
00302
00303 if (unique_write) {
00304 vector<int> unique_history = clust_seq.unique_history_order();
00305
00306 vector<int> inv_unique_history(clust_seq.history().size());
00307 for (unsigned int i = 0; i < unique_history.size(); i++) {
00308 inv_unique_history[unique_history[i]] = i;}
00309
00310 for (unsigned int i = 0; i < unique_history.size(); i++) {
00311 fj::ClusterSequence::history_element el =
00312 clust_seq.history()[unique_history[i]];
00313 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00314 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00315 printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00316 }
00317 }
00318
00319
00320
00321 if (show_cones) {
00322 const fj::SISConeExtras * extras =
00323 dynamic_cast<const fj::SISConeExtras *>(clust_seq.extras());
00324 cout << "most ambiguous split (difference in squared dist) = "
00325 << extras->most_ambiguous_split() << endl;
00326 }
00327 }
00328
00329 }
00330 }