FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
fastjet_timing.cc
1 //STARTHEADER
2 // $Id: fastjet_timing.cc 2577 2011-09-13 15:11:38Z salam $
3 //
4 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 
30 //----------------------------------------------------------------------
31 /// fastjet_timing.cc: Program to help time and test the fastjet package
32 ///
33 /// It reads files containing multiple events in the format
34 /// p1x p1y p1z E1
35 /// p2x p2y p2z E2
36 /// ...
37 /// #END
38 ///
39 /// An example input file containing 10 events is included as
40 /// data/Pythia-PtMin1000-LHC-10ev.dat
41 ///
42 /// Usage:
43 /// fastjet_timing [-strategy NUMBER] [-repeat nrepeats] [-massive] \
44 /// [-combine nevents] [-r Rparameter] [-incl ptmin] [...] \
45 /// < data_file
46 ///
47 /// where the clustering can be repeated to aid timing and multiple
48 /// events can be combined to get to larger multiplicities. Some options:
49 ///
50 /// -strategy N indicate stratgey from the enum FjStrategy (see
51 /// FjClusterSequence.hh).
52 ///
53 /// -combine nev for combining multiple events from the data file in order
54 /// to get to large multiplicities.
55 ///
56 /// -incl ptmin output of all inclusive jets with pt > ptmin is obtained
57 /// with the -incl option.
58 ///
59 /// -excld dcut output of all exclusive jets as obtained in a clustering
60 /// with dcut
61 ///
62 /// -massless read in only the 3-momenta and deduce energies assuming
63 /// that particles are massless
64 ///
65 /// -write for writing out detailed clustering sequence (valuable
66 /// for testing purposes)
67 ///
68 /// -unique_write writes out the sequence of dij's according to the
69 /// "unique_history_order" (useful for verifying consistency
70 /// between different clustering strategies).
71 ///
72 /// -cam switch to the inclusive Cambridge/Aachen algorithm --
73 /// note that the option -excld dcut provides a clustering
74 /// up to the dcut which is the minimum squared
75 /// distance between any pair of jets.
76 ///
77 #include "fastjet/PseudoJet.hh"
78 #include "fastjet/ClusterSequence.hh"
79 #include<iostream>
80 #include<sstream>
81 #include<valarray>
82 #include<vector>
83 #include <cstdio>
84 #include <cstdlib>
85 #include<cstddef> // for size_t
86 #include "CmdLine.hh"
87 
88 using namespace std;
89 
90 // to avoid excessive typing, define an abbreviation for the
91 // fastjet namespace
92 namespace fj = fastjet;
93 
94 inline double pow2(const double x) {return x*x;}
95 
96 /// a program to test and time the kt algorithm as implemented in fastjet
97 int main (int argc, char ** argv) {
98 
99  CmdLine cmdline(argc,argv);
100  // allow the use to specify the fj::Strategy either through the
101  // -clever or the -strategy options (both will take numerical
102  // values); the latter will override the former.
103  fj::Strategy strategy = fj::Strategy(cmdline.int_val("-strategy",
104  cmdline.int_val("-clever", fj::Best)));
105  int repeat = cmdline.int_val("-repeat",1);
106  int combine = cmdline.int_val("-combine",1);
107  bool write = cmdline.present("-write");
108  bool unique_write = cmdline.present("-unique_write");
109  bool hydjet = cmdline.present("-hydjet");
110  double ktR = cmdline.double_val("-r",1.0);
111  double inclkt = cmdline.double_val("-incl",-1.0);
112  int excln = cmdline.int_val ("-excln",-1);
113  double excld = cmdline.double_val("-excld",-1.0);
114  double etamax = cmdline.double_val("-etamax",1.0e305);
115  bool show_constituents = cmdline.present("-const");
116  bool massless = cmdline.present("-massless");
117  int nev = cmdline.int_val("-nev",1);
118  bool add_dense_coverage = cmdline.present("-dense");
119 
120  // The following option causes the Cambridge algo to be used.
121  // Note that currently the only output that works sensibly here is
122  // "-incl 0"
123  fj::JetAlgorithm jet_algorithm;
124  if (cmdline.present("-cam")) {
125  jet_algorithm = fj::cambridge_algorithm;
126  } else {
127  jet_algorithm = fj::kt_algorithm;
128  }
129 
130  if (!cmdline.all_options_used()) {cerr <<
131  "Error: some options were not recognized"<<endl;
132  exit(-1);}
133 
134 
135  for (int iev = 0; iev < nev; iev++) {
136  vector<fj::PseudoJet> jets;
137  string line;
138  int ndone = 0;
139  while (getline(cin, line)) {
140  //cout << line<<endl;
141  istringstream linestream(line);
142  if (line == "#END") {
143  ndone += 1;
144  if (ndone == combine) {break;}
145  }
146  if (line.substr(0,1) == "#") {continue;}
147  valarray<double> fourvec(4);
148  if (hydjet) {
149  // special reading from hydjet.txt event record (though actually
150  // this is supposed to be a standard pythia event record, so
151  // being able to read from it is perhaps not so bad an idea...)
152  int ii, istat,id,m1,m2,d1,d2;
153  double mass;
154  linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
155  >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
156  // current file contains mass of particle as 4th entry
157  if (istat == 1) {
158  fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
159  +pow2(fourvec[2])+pow2(mass));
160  }
161  } else {
162  if (massless) {
163  linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
164  fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
165  else {
166  linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
167  }
168  }
169  fj::PseudoJet psjet(fourvec);
170  if (abs(psjet.rap() < etamax)) {jets.push_back(psjet);}
171  }
172 
173  // add a fake underlying event which is very soft, uniformly distributed
174  // in eta,phi so as to allow one to reconstruct the area that is associated
175  // with each jet.
176  if (add_dense_coverage) {
177  srand(2);
178  int nphi = 60;
179  int neta = 100;
180  double kt = 1e-1;
181  for (int iphi = 0; iphi<nphi; iphi++) {
182  for (int ieta = -neta; ieta<neta+1; ieta++) {
183  double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
184  double eta = ieta * (10.0/neta) + rand()*0.001/RAND_MAX;
185  kt = 0.0000001*(1+rand()*0.1/RAND_MAX);
186  double pminus = kt*exp(-eta);
187  double pplus = kt*exp(+eta);
188  double px = kt*sin(phi);
189  double py = kt*cos(phi);
190  //cout << kt<<" "<<eta<<" "<<phi<<"\n";
191  fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
192  jets.push_back(mom);
193  }
194  }
195  }
196 
197  fj::JetDefinition jet_def(jet_algorithm, ktR, strategy);
198 
199  for (int irepeat = 0; irepeat < repeat ; irepeat++) {
200  fj::ClusterSequence clust_seq(jets,jet_def,write);
201  if (irepeat != 0) {continue;}
202  cout << "iev "<<iev<< ": number of particles = "<< jets.size() << endl;
203  cout << "strategy used = "<< clust_seq.strategy_string()<< endl;
204 
205  // now provide some nice output...
206  if (inclkt >= 0.0) {
207  vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.inclusive_jets(inclkt));
208  for (size_t j = 0; j < jets.size(); j++) {
209  printf("%5u %15.8f %15.8f %15.8f\n",j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
210  if (show_constituents) {
211  vector<fj::PseudoJet> const_jets = jets[j].constituents();
212  for (size_t k = 0; k < const_jets.size(); k++) {
213  printf(" jet%03u %15.8f %15.8f %15.8f\n",j,const_jets[k].rap(),
214  const_jets[k].phi(),sqrt(const_jets[k].kt2()));
215  }
216  cout << "\n\n";
217  }
218  }
219  }
220 
221  if (excln > 0) {
222  vector<fj::PseudoJet> jets = sorted_by_E(clust_seq.exclusive_jets(excln));
223 
224  cout << "Printing "<<excln<<" exclusive jets\n";
225  for (size_t j = 0; j < jets.size(); j++) {
226  printf("%5u %15.8f %15.8f %15.8f\n",
227  //j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
228  j,jets[j].rap(),jets[j].phi(),jets[j].kt2());
229  }
230  }
231 
232  if (excld > 0.0) {
233  vector<fj::PseudoJet> jets = sorted_by_pt(clust_seq.exclusive_jets(excld));
234  cout << "Printing exclusive jets for d = "<<excld<<"\n";
235  for (size_t j = 0; j < jets.size(); j++) {
236  printf("%5u %15.8f %15.8f %15.8f\n",
237  j,jets[j].rap(),jets[j].phi(),sqrt(jets[j].kt2()));
238  }
239  }
240 
241  // useful for testing that recombination sequences are unique
242  if (unique_write) {
243  vector<int> unique_history = clust_seq.unique_history_order();
244  // construct the inverse of the above mapping
245  vector<int> inv_unique_history(clust_seq.history().size());
246  for (unsigned int i = 0; i < unique_history.size(); i++) {
247  inv_unique_history[unique_history[i]] = i;}
248 
249  for (unsigned int i = 0; i < unique_history.size(); i++) {
250  fj::ClusterSequence::history_element el =
251  clust_seq.history()[unique_history[i]];
252  int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
253  int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
254  printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
255  }
256  }
257 
258  } // irepeat
259 
260  } // iev
261 }
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:799
deals with clustering
the longitudinally invariant kt algorithm
int main()
an example program showing how to use fastjet
Definition: 01-basic.cc:50
the FastJet namespace
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
Definition: PseudoJet.cc:815
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
JetAlgorithm
the various families of jet-clustering algorithm
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
class that is intended to hold a full definition of the jet clusterer