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