#include <iostream>
#include <sstream>
#include <valarray>
#include <vector>
#include <cstddef>
#include "CmdLine.hh"
#include "fastjet/internal/numconsts.hh"
#include "KtJet/KtEvent.h"
#include "KtJet/KtLorentzVector.h"
Go to the source code of this file.
Functions | |
double | pow2 (const double x) |
int | main (int argc, char **argv) |
a program to test and time the kt algorithm as implemented in ktjet |
int main | ( | int | argc, | |
char ** | argv | |||
) |
a program to test and time the kt algorithm as implemented in ktjet
Retrieve the final state jets from KtEvent sorted by Pt
Print out jets 4-momentum and Pt
Definition at line 27 of file ktjet_timing.cc.
References CmdLine::double_val(), CmdLine::int_val(), fastjet::d0::inline_maths::phi(), pow2(), CmdLine::present(), and fastjet::twopi.
00027 { 00028 00029 CmdLine cmdline(argc,argv); 00030 //bool clever = !cmdline.present("-dumb"); 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 //cout << line<<endl; 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 // set KtEvent flags 00068 int type, angle, recom; 00069 ostringstream info; 00070 if (cmdline.present("-eekt")) { 00071 type = 1; // e+e- 00072 angle = 1; // angular 00073 recom = 1; // E 00074 info << "Algorithm: KtJet e+e- kt algorithm" ; 00075 } else { 00076 type = 4; // PP 00077 angle = 2; // delta R 00078 recom = 1; // E 00079 info << "Algorithm: KtJet (long.inv.) with R = " << ktR ; 00080 } 00081 //double rparameter = 1.0; 00082 00083 for (int i = 0; i < repeat ; i++) { 00084 // Construct the KtEvent object 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 // Print out the number of final state jets 00093 //std::cout << "Number of final state jets: " << ev.getNJets() << std::endl; 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 // Retrieve the final state jets from KtEvent sorted by Pt 00107 std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt(); 00108 00109 // Print out index, rap, phi, pt 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 }
double pow2 | ( | const double | x | ) | [inline] |
Definition at line 24 of file ktjet_timing.cc.