FastJet  3.1.0-beta.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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"
22 using namespace std;
23 using namespace KtJet;
24 
25 inline double pow2(const double x) {return x*x;}
26 
27 /// a program to test and time the kt algorithm as implemented in ktjet
28 int 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 }