1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
//----------------------------------------------------------------------
// ktjet_timing.cc: Program similar to fastjet_timing.cc, to be used
// for timing and result comparisons. For usage, see
// explanations given at the start of fastjet_timing.cc
//
// Note: - some options present in fastjet_timing are missing
// - one or two behave differently, notably -write
//
#include<iostream>
#include<sstream>
#include<valarray>
#include<vector>
#include<cstddef> // for size_t
#include "CmdLine.hh"
#include "fastjet/internal/numconsts.hh"
/** Need to include these KtJet Headers */
#include "KtJet/KtEvent.h"
#include "KtJet/KtLorentzVector.h"
using namespace std;
using namespace KtJet;
inline double pow2(const double x) {return x*x;}
/// a program to test and time the kt algorithm as implemented in ktjet
int main (int argc, char ** argv) {
CmdLine cmdline(argc,argv);
//bool clever = !cmdline.present("-dumb");
int repeat = cmdline.int_val("-repeat",1);
int combine = cmdline.int_val("-combine",1);
bool write = cmdline.present("-write");
double ktR = cmdline.double_val("-r",1.0);
double inclkt = cmdline.double_val("-incl",-1.0);
int excln = cmdline.int_val ("-excln",-1);
double excld = cmdline.double_val("-excld",-1.0);
int nev = cmdline.int_val("-nev",1);
bool massless = cmdline.present("-massless");
bool get_all_dij = cmdline.present("-get-all-dij");
for (int iev = 0; iev < nev; iev++) {
vector<KtJet::KtLorentzVector> jets;
string line;
int ndone = 0;
while (getline(cin, line)) {
//cout << line<<endl;
istringstream linestream(line);
if (line == "#END") {
ndone += 1;
if (ndone == combine) {break;}
}
if (line.substr(0,1) == "#") {continue;}
valarray<double> fourvec(4);
linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
if (massless) {
linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
else {
linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
}
KtJet::KtLorentzVector p(fourvec[0],fourvec[1],fourvec[2],fourvec[3]);
jets.push_back(p);
}
// set KtEvent flags
int type, angle, recom;
ostringstream info;
if (cmdline.present("-eekt")) {
type = 1; // e+e-
angle = 1; // angular
recom = 1; // E
info << "Algorithm: KtJet e+e- kt algorithm" ;
} else {
type = 4; // PP
angle = 2; // delta R
recom = 1; // E
info << "Algorithm: KtJet (long.inv.) with R = " << ktR ;
}
//double rparameter = 1.0;
for (int i = 0; i < repeat ; i++) {
// Construct the KtEvent object
KtJet::KtEvent ev(jets,type,angle,recom,ktR);
if (i!=0) {continue;}
int nparticles = jets.size();
cout << "Number of particles = "<< nparticles << endl;
cout << info.str() << endl;
// Print out the number of final state jets
//std::cout << "Number of final state jets: " << ev.getNJets() << std::endl;
if (write) {
/** Retrieve the final state jets from KtEvent sorted by Pt*/
std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
/** Print out jets 4-momentum and Pt */
std::vector<KtJet::KtLorentzVector>::const_iterator itr = jets.begin();
for( ; itr != jets.end() ; ++itr) {
std::cout << "Jets Pt2: " << pow2((*itr).perp()) << std::endl;
}
}
if (inclkt >= 0.0) {
// Retrieve the final state jets from KtEvent sorted by Pt
std::vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
// Print out index, rap, phi, pt
for (size_t j = 0; j < jets.size(); j++) {
if (jets[j].perp() < inclkt) {break;}
double phi = jets[j].phi();
if (phi < 0.0) {phi += fastjet::twopi;}
printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rapidity(),phi,jets[j].perp());
}
}
if (excln > 0) {
ev.findJetsN(excln);
vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
cout << "Printing "<<excln<<" exclusive jets\n";
for (size_t j = 0; j < jets.size(); j++) {
double phi = jets[j].phi();
if (phi < 0) phi += fastjet::twopi;
printf("%5u %15.8f %15.8f %15.8f\n",j,
jets[j].rapidity(),phi,jets[j].perp());
}
}
if (excld > 0.0) {
ev.findJetsD(excld);
vector<KtJet::KtLorentzVector> jets = ev.getJetsPt();
cout << "Printing exclusive jets for d = "<<excld<<"\n";
for (size_t j = 0; j < jets.size(); j++) {
double phi = jets[j].phi();
if (phi < 0) phi += fastjet::twopi;
printf("%5u %15.8f %15.8f %15.8f\n",j,
jets[j].rapidity(),phi,jets[j].perp());
}
}
if (get_all_dij) {
for (int i = nparticles-1; i > 0; i--) {
printf("d for n = %4d -> %4d is %14.5e\n", i+1, i, ev.getDMerge(i));
}
}
}
}
}
|
