ktjet_timing.cc
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009 #include<iostream>
00010 #include<sstream>
00011 #include<valarray>
00012 #include<vector>
00013 #include<cstddef>
00014 #include "CmdLine.hh"
00015 #include "fastjet/internal/numconsts.hh"
00016
00017
00019 #include "KtJet/KtEvent.h"
00020 #include "KtJet/KtLorentzVector.h"
00021 using namespace std;
00022 using namespace KtJet;
00023
00024 inline double pow2(const double x) {return x*x;}
00025
00027 int main (int argc, char ** argv) {
00028
00029 CmdLine cmdline(argc,argv);
00030
00031 int repeat = cmdline.int_val("-repeat",1);
00032 int combine = cmdline.int_val("-combine",1);
00033 bool write = cmdline.present("-write");
00034 double ktR = cmdline.double_val("-r",1.0);
00035 double inclkt = cmdline.double_val("-incl",-1.0);
00036 int excln = cmdline.int_val ("-excln",-1);
00037 double excld = cmdline.double_val("-excld",-1.0);
00038 int nev = cmdline.int_val("-nev",1);
00039 bool massless = cmdline.present("-massless");
00040 bool get_all_dij = cmdline.present("-get-all-dij");
00041
00042
00043 for (int iev = 0; iev < nev; iev++) {
00044 vector<KtJet::KtLorentzVector> jets;
00045 string line;
00046 int ndone = 0;
00047 while (getline(cin, line)) {
00048
00049 istringstream linestream(line);
00050 if (line == "#END") {
00051 ndone += 1;
00052 if (ndone == combine) {break;}
00053 }
00054 if (line.substr(0,1) == "#") {continue;}
00055 valarray<double> fourvec(4);
00056 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00057 if (massless) {
00058 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00059 fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00060 else {
00061 linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00062 }
00063 KtJet::KtLorentzVector p(fourvec[0],fourvec[1],fourvec[2],fourvec[3]);
00064 jets.push_back(p);
00065 }
00066
00067
00068 int type, angle, recom;
00069 ostringstream info;
00070 if (cmdline.present("-eekt")) {
00071 type = 1;
00072 angle = 1;
00073 recom = 1;
00074 info << "Algorithm: KtJet e+e- kt algorithm" ;
00075 } else {
00076 type = 4;
00077 angle = 2;
00078 recom = 1;
00079 info << "Algorithm: KtJet (long.inv.) with R = " << ktR ;
00080 }
00081
00082
00083 for (int i = 0; i < repeat ; i++) {
00084
00085 KtJet::KtEvent ev(jets,type,angle,recom,ktR);
00086
00087 if (i!=0) {continue;}
00088 int nparticles = jets.size();
00089 cout << "Number of particles = "<< nparticles << endl;
00090 cout << info.str() << endl;
00091
00092
00093
00094 if (write) {
00096 std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00097
00099 std::vector<KtJet::KtLorentzVector>::const_iterator itr = jets.begin();
00100 for( ; itr != jets.end() ; ++itr) {
00101 std::cout << "Jets Pt2: " << pow2((*itr).perp()) << std::endl;
00102 }
00103 }
00104
00105 if (inclkt >= 0.0) {
00106
00107 std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00108
00109
00110 for (size_t j = 0; j < jets.size(); j++) {
00111 if (jets[j].perp() < inclkt) {break;}
00112 double phi = jets[j].phi();
00113 if (phi < 0.0) {phi += fastjet::twopi;}
00114 printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rapidity(),phi,jets[j].perp());
00115 }
00116 }
00117
00118 if (excln > 0) {
00119 ev.findJetsN(excln);
00120 vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00121 cout << "Printing "<<excln<<" exclusive jets\n";
00122 for (size_t j = 0; j < jets.size(); j++) {
00123 double phi = jets[j].phi();
00124 if (phi < 0) phi += fastjet::twopi;
00125 printf("%5u %15.8f %15.8f %15.8f\n",j,
00126 jets[j].rapidity(),phi,jets[j].perp());
00127 }
00128 }
00129
00130 if (excld > 0.0) {
00131 ev.findJetsD(excld);
00132 vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
00133 cout << "Printing exclusive jets for d = "<<excld<<"\n";
00134 for (size_t j = 0; j < jets.size(); j++) {
00135 double phi = jets[j].phi();
00136 if (phi < 0) phi += fastjet::twopi;
00137 printf("%5u %15.8f %15.8f %15.8f\n",j,
00138 jets[j].rapidity(),phi,jets[j].perp());
00139 }
00140 }
00141
00142 if (get_all_dij) {
00143 for (int i = nparticles-1; i > 0; i--) {
00144 printf("d for n = %4d -> %4d is %14.5e\n", i+1, i, ev.getDMerge(i));
00145 }
00146 }
00147
00148 }
00149
00150 }
00151 }