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