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 <cstdlib>
00086 #include<cstddef>
00087 #include "CmdLine.hh"
00088
00089 using namespace std;
00090
00091
00092
00093 namespace fj = fastjet;
00094
00095 inline double pow2(const double x) {return x*x;}
00096
00098 int main (int argc, char ** argv) {
00099
00100 CmdLine cmdline(argc,argv);
00101
00102
00103
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
00122
00123
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
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
00151
00152
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
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
00175
00176
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
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
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
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
00243 if (unique_write) {
00244 vector<int> unique_history = clust_seq.unique_history_order();
00245
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 }
00260
00261 }
00262 }