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