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
00079 #include "fastjet/PseudoJet.hh"
00080 #include "fastjet/ClusterSequence.hh"
00081 #include<iostream>
00082 #include<sstream>
00083 #include<valarray>
00084 #include<vector>
00085 #include <cstdio>
00086 #include <cstdlib>
00087 #include<cstddef>
00088 #include "CmdLine.hh"
00089
00090 using namespace std;
00091
00092
00093
00094 namespace fj = fastjet;
00095
00096 inline double pow2(const double x) {return x*x;}
00097
00099 int main (int argc, char ** argv) {
00100
00101 CmdLine cmdline(argc,argv);
00102
00103
00104
00105 fj::Strategy strategy = fj::Strategy(cmdline.int_val("-strategy",
00106 cmdline.int_val("-clever", fj::Best)));
00107 int repeat = cmdline.int_val("-repeat",1);
00108 int combine = cmdline.int_val("-combine",1);
00109 bool write = cmdline.present("-write");
00110 bool unique_write = cmdline.present("-unique_write");
00111 bool hydjet = cmdline.present("-hydjet");
00112 double ktR = cmdline.double_val("-r",1.0);
00113 double inclkt = cmdline.double_val("-incl",-1.0);
00114 int excln = cmdline.int_val ("-excln",-1);
00115 double excld = cmdline.double_val("-excld",-1.0);
00116 double etamax = cmdline.double_val("-etamax",1.0e305);
00117 bool show_constituents = cmdline.present("-const");
00118 bool massless = cmdline.present("-massless");
00119 int nev = cmdline.int_val("-nev",1);
00120 bool add_dense_coverage = cmdline.present("-dense");
00121
00122
00123
00124
00125 fj::JetAlgorithm jet_algorithm;
00126 if (cmdline.present("-cam")) {
00127 jet_algorithm = fj::cambridge_algorithm;
00128 } else {
00129 jet_algorithm = fj::kt_algorithm;
00130 }
00131
00132 if (!cmdline.all_options_used()) {cerr <<
00133 "Error: some options were not recognized"<<endl;
00134 exit(-1);}
00135
00136
00137 for (int iev = 0; iev < nev; iev++) {
00138 vector<fj::PseudoJet> jets;
00139 string line;
00140 int ndone = 0;
00141 while (getline(cin, line)) {
00142
00143 istringstream linestream(line);
00144 if (line == "#END") {
00145 ndone += 1;
00146 if (ndone == combine) {break;}
00147 }
00148 if (line.substr(0,1) == "#") {continue;}
00149 valarray<double> fourvec(4);
00150 if (hydjet) {
00151
00152
00153
00154 int ii, istat,id,m1,m2,d1,d2;
00155 double mass;
00156 linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00157 >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00158
00159 if (istat == 1) {
00160 fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00161 +pow2(fourvec[2])+pow2(mass));
00162 }
00163 } else {
00164 if (massless) {
00165 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00166 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00167 else {
00168 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00169 }
00170 }
00171 fj::PseudoJet psjet(fourvec);
00172 if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
00173 }
00174
00175
00176
00177
00178 if (add_dense_coverage) {
00179 srand(2);
00180 int nphi = 60;
00181 int neta = 100;
00182 double kt = 1e-1;
00183 for (int iphi = 0; iphi<nphi; iphi++) {
00184 for (int ieta = -neta; ieta<neta+1; ieta++) {
00185 double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00186 double eta = ieta * (10.0/neta) + rand()*0.001/RAND_MAX;
00187 kt = 0.0000001*(1+rand()*0.1/RAND_MAX);
00188 double pminus = kt*exp(-eta);
00189 double pplus = kt*exp(+eta);
00190 double px = kt*sin(phi);
00191 double py = kt*cos(phi);
00192
00193 fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00194 jets.push_back(mom);
00195 }
00196 }
00197 }
00198
00199 fj::JetDefinition jet_def(jet_algorithm, ktR, strategy);
00200
00201 for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00202 fj::ClusterSequence clust_seq(jets,jet_def,write);
00203 if (irepeat != 0) {continue;}
00204 cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl;
00205 cout << "strategy used = "<< clust_seq.strategy_string()<< endl;
00206
00207
00208 if (inclkt >= 0.0) {
00209 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt));
00210 for (size_t j = 0; j < jets.size(); j++) {
00211 printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00212 if (show_constituents) {
00213 vector<fj::PseudoJet> const_jets = clust_seq.constituents(jets[j]);
00214 for (size_t k = 0; k < const_jets.size(); k++) {
00215 printf(" jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
00216 const_jets[k].phi(),sqrt(const_jets[k].kt2()));
00217 }
00218 cout << "\n\n";
00219 }
00220 }
00221 }
00222
00223 if (excln > 0) {
00224 vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
00225
00226 cout << "Printing "<<excln<<" exclusive jets\n";
00227 for (size_t j = 0; j < jets.size(); j++) {
00228 printf("%5u %15.8f %15.8f %15.8f\n",
00229
00230 j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
00231 }
00232 }
00233
00234 if (excld > 0.0) {
00235 vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
00236 cout << "Printing exclusive jets for d = "<<excld<<"\n";
00237 for (size_t j = 0; j < jets.size(); j++) {
00238 printf("%5u %15.8f %15.8f %15.8f\n",
00239 j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
00240 }
00241 }
00242
00243
00244 if (unique_write) {
00245 vector<int> unique_history = clust_seq.unique_history_order();
00246
00247 vector<int> inv_unique_history(clust_seq.history().size());
00248 for (unsigned int i = 0; i < unique_history.size(); i++) {
00249 inv_unique_history[unique_history[i]] = i;}
00250
00251 for (unsigned int i = 0; i < unique_history.size(); i++) {
00252 fj::ClusterSequence::history_element el =
00253 clust_seq.history()[unique_history[i]];
00254 int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00255 int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00256 printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00257 }
00258 }
00259
00260 }
00261
00262 }
00263 }