FastJet  3.0.6
fastjet_timing_plugins.cc
1 //STARTHEADER
2 // $Id: fastjet_timing_plugins.cc 3209 2013-09-15 20:03:30Z 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 /// Options for reading
51 /// -------------------
52 ///
53 /// -nev n number of events to run
54 ///
55 /// -combine n for combining multiple events from the data file in order
56 /// to get a single high-multipicity event to run.
57 ///
58 /// -massless read in only the 3-momenta and deduce energies assuming
59 /// that particles are massless
60 ///
61 /// -dense adds dense ghost coverage
62 ///
63 /// -repeat n repeats each event n times
64 ///
65 /// -nhardest n keep only the n hardest particles in the event
66 ///
67 /// -file name read from the corresponding file rather than stdin.
68 /// (The file will be reopened for each new jet alg.; in
69 /// constrast, if you use stdin, each new alg will take a
70 /// new event).
71 ///
72 /// Output Options
73 /// --------------
74 ///
75 /// -incl ptmin output of all inclusive jets with pt > ptmin is obtained
76 /// with the -incl option.
77 ///
78 /// -repeat-incl ptmin
79 /// same as -incl ptmin but do it for each repetition
80 /// of the clustering
81 ///
82 /// -excld dcut output of all exclusive jets as obtained in a clustering
83 /// with dcut
84 ///
85 /// -excly ycut output of all exclusive jets as obtained in a clustering
86 /// with ycut
87 ///
88 /// -excln n output of clustering to n exclusive jets
89 ///
90 /// -ee-print print things as px,py,pz,E
91 ///
92 /// -get-all-dij print out all dij values
93 /// -get-all-yij print out all yij values
94 ///
95 /// -const show jet constituents (works with excl jets)
96 ///
97 /// -write for writing out detailed clustering sequence (valuable
98 /// for testing purposes)
99 ///
100 /// -unique_write writes out the sequence of dij's according to the
101 /// "unique_history_order" (useful for verifying consistency
102 /// between different clustering strategies).
103 ///
104 /// -root file sends output to file that can be read in with the script in
105 /// root/ so as to show a lego-plot of the event
106 ///
107 /// -cones show extra info about internal steps for SISCone
108 ///
109 /// -area calculate areas. Additional options include
110 /// -area:active
111 /// -area:passive
112 /// -area:explicit
113 /// -area:voronoi Rfact
114 /// -area:repeat nrepeat
115 /// -ghost-area area
116 /// -ghost-maxrap maxrap
117 /// -area:fj2 place ghosts as in fj2
118 ///
119 /// -bkgd calculate the background density. Additional options include
120 /// -bkgd:csab use the old ClusterSequenceAreaBase methods
121 /// -bkgd:jetmedian use the new JetMedianBackgroundEstimator class
122 /// -bkgd:fj2 force jetmedian to calculate sigma as in fj2
123 /// -bkgd:gridmedian use GridMedianBackgroundEstimator with grid up to ghost_maxrap-ktR and grid spacing of 2ktR
124 ///
125 /// Algorithms
126 /// ----------
127 /// -all-algs runs all algorithms
128 ///
129 /// -kt switch to the longitudinally invariant kt algorithm
130 /// Note: this is the default one.
131 ///
132 /// -cam switch to the inclusive Cambridge/Aachen algorithm --
133 /// note that the option -excld dcut provides a clustering
134 /// up to the dcut which is the minimum squared
135 /// distance between any pair of jets.
136 ///
137 /// -antikt switch to the anti-kt clustering algorithm
138 ///
139 /// -genkt switch to the genkt algorithm
140 /// you can provide the parameter of the alg as an argument to
141 /// -genkt (1 by default)
142 ///
143 /// -eekt switch to the e+e- kt algorithm
144 ///
145 /// -eegenkt switch to the genkt algorithm
146 /// you can provide the parameter of the alg as an argument to
147 /// -ee_genkt (1 by default)
148 ///
149 /// plugins (don't delete this line)
150 ///
151 /// -pxcone switch to the PxCone jet algorithm
152 ///
153 /// -siscone switch to the SISCone jet algorithm (seedless cones)
154 /// -sisconespheri switch to the Spherical SISCone jet algorithm (seedless cones)
155 ///
156 /// -midpoint switch to CDF's midpoint code
157 /// -jetclu switch to CDF's jetclu code
158 ///
159 /// -d0runipre96cone switch to the D0RunIpre96Cone plugin
160 /// -d0runicone switch to the D0RunICone plugin
161 ///
162 /// -d0runiicone switch to D0's run II midpoint cone
163 ///
164 /// -trackjet switch to the TrackJet plugin
165 ///
166 /// -atlascone switch to the ATLASCone plugin
167 ///
168 /// -eecambridge switch to the EECambridge plugin
169 ///
170 /// -jade switch to the Jade plugin
171 ///
172 /// -cmsiterativecone switch to the CMSIterativeCone plugin
173 ///
174 /// -gridjet switch to the GridJet plugin
175 ///
176 /// end of plugins (don't delete this line)
177 ///
178 ///
179 /// Options for running algs
180 /// ------------------------
181 ///
182 /// -r sets the radius of the jet algorithm (default = 1.0)
183 ///
184 /// -overlap | -f sets the overlap fraction in cone algs with split-merge
185 ///
186 /// -seed sets the seed threshold
187 ///
188 /// -strategy N indicate stratgey from the enum fastjet::Strategy (see
189 /// fastjet/JetDefinition.hh).
190 ///
191 
192 
193 #include "fastjet/ClusterSequenceArea.hh"
194 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
195 #include "fastjet/tools/GridMedianBackgroundEstimator.hh"
196 #include "fastjet/Selector.hh"
197 #include<iostream>
198 #include<sstream>
199 #include<fstream>
200 #include<valarray>
201 #include<vector>
202 #include <cstdlib>
203 //#include<cstddef> // for size_t
204 #include "CmdLine.hh"
205 
206 // get info on how fastjet was configured
207 #include "fastjet/config.h"
208 
209 // include the installed plugins (don't delete this line)
210 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
211 #include "fastjet/SISConePlugin.hh"
212 #include "fastjet/SISConeSphericalPlugin.hh"
213 #endif
214 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
215 #include "fastjet/CDFMidPointPlugin.hh"
216 #include "fastjet/CDFJetCluPlugin.hh"
217 #endif
218 #ifdef FASTJET_ENABLE_PLUGIN_PXCONE
219 #include "fastjet/PxConePlugin.hh"
220 #endif
221 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
222 #include "fastjet/D0RunIIConePlugin.hh"
223 #endif
224 #ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
225 #include "fastjet/TrackJetPlugin.hh"
226 #endif
227 #ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
228 #include "fastjet/ATLASConePlugin.hh"
229 #endif
230 #ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
231 #include "fastjet/EECambridgePlugin.hh"
232 #endif
233 #ifdef FASTJET_ENABLE_PLUGIN_JADE
234 #include "fastjet/JadePlugin.hh"
235 #endif
236 #ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
237 #include "fastjet/CMSIterativeConePlugin.hh"
238 #endif
239 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
240 #include "fastjet/D0RunIpre96ConePlugin.hh"
241 #include "fastjet/D0RunIConePlugin.hh"
242 #endif
243 #ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
244 #include "fastjet/GridJetPlugin.hh"
245 #endif
246 // end of installed plugins inclusion (don't delete this line)
247 
248 using namespace std;
249 
250 // to avoid excessive typing, define an abbreviation for the
251 // fastjet namespace
252 namespace fj = fastjet;
253 
254 inline double pow2(const double x) {return x*x;}
255 
256 // pretty print the jets and their subjets
257 void print_jets_and_sub (const vector<fj::PseudoJet> & jets, double dcut);
258 
259 string rootfile;
260 CmdLine * cmdline_p;
261 
262 bool do_areas;
263 
264 /// sort and pretty print jets, with exact behaviour depending on
265 /// whether ee_print is true or not
266 bool ee_print = false;
267 void print_jets(const vector<fj::PseudoJet> & jets, bool show_const = false);
268 
269 bool found_unavailable = false;
270 void is_unavailable(const string & algname) {
271  cerr << algname << " requested, but not available for this compilation" << endl;
272  found_unavailable = true;
273  //exit(0);
274 }
275 
276 
277 /// a program to test and time a range of algorithms as implemented or
278 /// wrapped in fastjet
279 int main (int argc, char ** argv) {
280 
282 
283  CmdLine cmdline(argc,argv);
284  cmdline_p = &cmdline;
285  // allow the use to specify the fj::Strategy either through the
286  // -clever or the -strategy options (both will take numerical
287  // values); the latter will override the former.
288  fj::Strategy strategy = fj::Strategy(cmdline.int_val("-strategy",
289  cmdline.int_val("-clever", fj::Best)));
290  int repeat = cmdline.int_val("-repeat",1);
291  int combine = cmdline.int_val("-combine",1);
292  bool write = cmdline.present("-write");
293  bool unique_write = cmdline.present("-unique_write");
294  bool hydjet = cmdline.present("-hydjet");
295  double ktR = cmdline.double_val("-r",1.0);
296  ktR = cmdline.double_val("-R",ktR); // allow -r and -R
297  double inclkt = cmdline.double_val("-incl",-1.0);
298  double repeatinclkt = cmdline.double_val("-repeat-incl",-1.0);
299  int excln = cmdline.int_val ("-excln",-1);
300  double excld = cmdline.double_val("-excld",-1.0);
301  double excly = cmdline.double_val("-excly",-1.0);
302  ee_print = cmdline.present("-ee-print");
303  bool get_all_dij = cmdline.present("-get-all-dij");
304  bool get_all_yij = cmdline.present("-get-all-yij");
305  double subdcut = cmdline.double_val("-subdcut",-1.0);
306  double rapmax = cmdline.double_val("-rapmax",1.0e305);
307  if (cmdline.present("-etamax")) {
308  cerr << "WARNING: -etamax options actually sets maximum rapidity (and overrides -rapmax)\n";
309  rapmax = cmdline.double_val("-etamax");
310  }
311  bool show_constituents = cmdline.present("-const");
312  bool massless = cmdline.present("-massless");
313  int nev = cmdline.int_val("-nev",1);
314  bool add_dense_coverage = cmdline.present("-dense");
315  double ghost_maxrap = cmdline.value("-ghost-maxrap",5.0);
316  bool all_algs = cmdline.present("-all-algs");
317 
318  fj::Selector particles_sel = (cmdline.present("-nhardest"))
319  ? fj::SelectorNHardest(cmdline.value<unsigned int>("-nhardest"))
320  : fj::SelectorIdentity();
321 
322  do_areas = cmdline.present("-area");
323  fj::AreaDefinition area_def;
324  if (do_areas) {
325  assert(!write); // it's incompatible
326  fj::GhostedAreaSpec ghost_spec(ghost_maxrap,
327  cmdline.value("-area:repeat", 1),
328  cmdline.value("-ghost-area", 0.01));
329  if (cmdline.present("-area:fj2")) ghost_spec.set_fj2_placement(true);
330  if (cmdline.present("-area:explicit")) {
331  area_def = fj::AreaDefinition(fj::active_area_explicit_ghosts, ghost_spec);
332  } else if (cmdline.present("-area:passive")) {
333  area_def = fj::AreaDefinition(fj::passive_area, ghost_spec);
334  } else if (cmdline.present("-area:voronoi")) {
335  double Rfact = cmdline.value<double>("-area:voronoi");
336  area_def = fj::AreaDefinition(fj::voronoi_area,
337  fj::VoronoiAreaSpec(Rfact));
338  } else {
339  cmdline.present("-area:active"); // allow, but do not require, arg
340  area_def = fj::AreaDefinition(fj::active_area, ghost_spec);
341  }
342  }
343  bool do_bkgd = cmdline.present("-bkgd"); // background estimation
344  bool do_bkgd_csab = false, do_bkgd_jetmedian = false, do_bkgd_fj2 = false;
345  bool do_bkgd_gridmedian = false;
346  fj::Selector bkgd_range;
347  if (do_bkgd) {
348  bkgd_range = fj::SelectorAbsRapMax(ghost_maxrap - ktR);
349  if (cmdline.present("-bkgd:csab")) {do_bkgd_csab = true;}
350  else if (cmdline.present("-bkgd:jetmedian")) {do_bkgd_jetmedian = true;
351  do_bkgd_fj2 = cmdline.present("-bkgd:fj2");
352  } else if (cmdline.present("-bkgd:gridmedian")) {do_bkgd_gridmedian = true;
353  } else {
354  throw fj::Error("with the -bkgd option, some particular background must be specified (csab or jetmedian)");
355  }
356  assert(do_areas || do_bkgd_gridmedian);
357  }
358 
359  bool show_cones = cmdline.present("-cones"); // only works for siscone
360 
361  // for cone algorithms
362  // allow -f and -overlap
363  double overlap_threshold = cmdline.double_val("-overlap",0.5);
364  overlap_threshold = cmdline.double_val("-f",overlap_threshold);
365  double seed_threshold = cmdline.double_val("-seed",1.0);
366 
367  // for ee algorithms, allow to specify ycut
368  double ycut = cmdline.double_val("-ycut",0.08);
369 
370  // for printing jets to a file for reading by root
371  rootfile = cmdline.value<string>("-root","");
372 
373  // out default scheme is the E_scheme
375 
376  // The following option causes the Cambridge algo to be used.
377  // Note that currently the only output that works sensibly here is
378  // "-incl 0"
379  vector<fj::JetDefinition> jet_defs;
380  if (all_algs || cmdline.present("-cam") || cmdline.present("-CA")) {
381  jet_defs.push_back( fj::JetDefinition(fj::cambridge_algorithm, ktR, scheme, strategy));
382  }
383  if (all_algs || cmdline.present("-antikt")) {
384  jet_defs.push_back( fj::JetDefinition(fj::antikt_algorithm, ktR, scheme, strategy));
385  }
386  if (all_algs || cmdline.present("-genkt")) {
387  double p;
388  if (cmdline.present("-genkt")) p = cmdline.value<double>("-genkt");
389  else p = -0.5;
390  jet_defs.push_back( fj::JetDefinition(fj::genkt_algorithm, ktR, p, scheme, strategy));
391  }
392  if (all_algs || cmdline.present("-eekt")) {
393  jet_defs.push_back( fj::JetDefinition(fj::ee_kt_algorithm));
394  }
395  if (all_algs || cmdline.present("-eegenkt")) {
396  double p;
397  if (cmdline.present("-eegenkt")) p = cmdline.value<double>("-eegenkt");
398  else p = -0.5;
399  jet_defs.push_back( fj::JetDefinition(fj::ee_genkt_algorithm, ktR, p, scheme, strategy));
400 
401 // checking if one asks to run a plugin (don't delete this line)
402  }
403  if (all_algs || cmdline.present("-midpoint")) {
404 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
405  typedef fj::CDFMidPointPlugin MPPlug; // for brevity
406  double cone_area_fraction = 1.0;
407  int max_pair_size = 2;
408  int max_iterations = 100;
409  MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
410  if (cmdline.present("-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
411  if (cmdline.present("-sm-pt")) sm_scale = MPPlug::SM_pt; // default
412  if (cmdline.present("-sm-mt")) sm_scale = MPPlug::SM_mt;
413  if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et;
414  jet_defs.push_back( fj::JetDefinition( new fj::CDFMidPointPlugin (
415  seed_threshold, ktR,
416  cone_area_fraction, max_pair_size,
417  max_iterations, overlap_threshold,
418  sm_scale)));
419 #else // FASTJET_ENABLE_PLUGIN_CDFCONES
420  is_unavailable("MidPoint");
421 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
422  }
423  if (all_algs || cmdline.present("-pxcone")) {
424 #ifdef FASTJET_ENABLE_PLUGIN_PXCONE
425  double min_jet_energy = 5.0;
426  jet_defs.push_back( fj::JetDefinition( new fj::PxConePlugin (
427  ktR, min_jet_energy,
428  overlap_threshold)));
429 #else // FASTJET_ENABLE_PLUGIN_PXCONE
430  is_unavailable("PxCone");
431 #endif // FASTJET_ENABLE_PLUGIN_PXCONE
432  }
433  if (all_algs || cmdline.present("-jetclu")) {
434 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
435  jet_defs.push_back( fj::JetDefinition( new fj::CDFJetCluPlugin (
436  ktR, overlap_threshold, seed_threshold)));
437 #else // FASTJET_ENABLE_PLUGIN_CDFCONES
438  is_unavailable("JetClu");
439 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
440  }
441  if (all_algs || cmdline.present("-siscone") || cmdline.present("-sisconespheri")) {
442 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
443  typedef fj::SISConePlugin SISPlug; // for brevity
444  int npass = cmdline.value("-npass",0);
445  if (all_algs || cmdline.present("-siscone")) {
446  double sisptmin = cmdline.value("-sisptmin",0.0);
447  SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass,sisptmin);
448  if (cmdline.present("-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
449  if (cmdline.present("-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
450  if (cmdline.present("-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
451  if (cmdline.present("-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
452  // cause it to use the jet-definition's own recombiner
453  plugin->set_use_jet_def_recombiner(true);
454  jet_defs.push_back( fj::JetDefinition(plugin));
455  }
456  if (all_algs || cmdline.present("-sisconespheri")) {
457  double sisEmin = cmdline.value("-sisEmin",0.0);
458  fj::SISConeSphericalPlugin * plugin =
459  new fj::SISConeSphericalPlugin(ktR, overlap_threshold,npass,sisEmin);
460  if (cmdline.present("-ghost-sep")) {
461  plugin->set_ghost_separation_scale(cmdline.value<double>("-ghost-sep"));
462  }
463  jet_defs.push_back( fj::JetDefinition(plugin));
464  }
465 #else // FASTJET_ENABLE_PLUGIN_SISCONE
466  is_unavailable("SISCone");
467 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
468  }
469  if (all_algs || cmdline.present("-d0runiicone")) {
470 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
471  double min_jet_Et = 6.0; // was 8 GeV in earlier work
472  jet_defs.push_back( fj::JetDefinition(new fj::D0RunIIConePlugin(ktR,min_jet_Et)));
473 #else // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
474  is_unavailable("D0RunIICone");
475 #endif // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
476  }
477  if (all_algs || cmdline.present("-trackjet")) {
478 #ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
479  jet_defs.push_back( fj::JetDefinition(new fj::TrackJetPlugin(ktR)));
480 #else // FASTJET_ENABLE_PLUGIN_TRACKJET
481  is_unavailable("TrackJet");
482 #endif // FASTJET_ENABLE_PLUGIN_TRACKJET
483  }
484  if (all_algs || cmdline.present("-atlascone")) {
485 #ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
486  jet_defs.push_back( fj::JetDefinition(new fj::ATLASConePlugin(ktR)));
487 #else // FASTJET_ENABLE_PLUGIN_ATLASCONE
488  is_unavailable("ATLASCone");
489 #endif // FASTJET_ENABLE_PLUGIN_ATLASCONE
490  }
491  if (all_algs || cmdline.present("-eecambridge")) {
492 #ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
493  jet_defs.push_back( fj::JetDefinition(new fj::EECambridgePlugin(ycut)));
494 #else // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
495  is_unavailable("EECambridge");
496 #endif // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
497  }
498  if (all_algs || cmdline.present("-jade")) {
499 #ifdef FASTJET_ENABLE_PLUGIN_JADE
500  jet_defs.push_back( fj::JetDefinition(new fj::JadePlugin()));
501 #else // FASTJET_ENABLE_PLUGIN_JADE
502  is_unavailable("Jade");
503 #endif // FASTJET_ENABLE_PLUGIN_JADE
504  }
505  if (all_algs || cmdline.present("-cmsiterativecone")) {
506 #ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
507  jet_defs.push_back( fj::JetDefinition(new fj::CMSIterativeConePlugin(ktR,seed_threshold)));
508 #else // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
509  is_unavailable("CMSIterativeCone");
510 #endif // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
511  }
512  if (all_algs || cmdline.present("-d0runipre96cone")) {
513 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
514  jet_defs.push_back( fj::JetDefinition(new fj::D0RunIpre96ConePlugin(ktR, seed_threshold, overlap_threshold)));
515 #else // FASTJET_ENABLE_PLUGIN_D0RUNICONE
516  is_unavailable("D0RunIpre96Cone");
517 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
518  }
519  if (all_algs || cmdline.present("-d0runicone")) {
520 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
521  jet_defs.push_back( fj::JetDefinition(new fj::D0RunIConePlugin(ktR, seed_threshold, overlap_threshold)));
522 #else // FASTJET_ENABLE_PLUGIN_D0RUNICONE
523  is_unavailable("D0RunICone");
524 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
525  }
526  if (all_algs || cmdline.present("-gridjet")) {
527 #ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
528  // we want a grid_ymax of 5.0, but when using R=0.4 (i.e. grid
529  // spacing of 0.8), this leads to 12.5 grid cells; depending on
530  // whether this is 12.499999999999 or 12.5000000....1 this gets
531  // converted either to 12 or 13, making the results sensitive to
532  // rounding errors.
533  //
534  // Instead we therefore take 4.9999999999, which avoids this problem.
535  double grid_ymax = 4.9999999999;
536  jet_defs.push_back( fj::JetDefinition(new fj::GridJetPlugin(grid_ymax, ktR*2.0)));
537 #else // FASTJET_ENABLE_PLUGIN_GRIDJET
538  is_unavailable("GridJet");
539 #endif // FASTJET_ENABLE_PLUGIN_GRIDJET
540 // end of checking if one asks to run a plugin (don't delete this line)
541  }
542  if (all_algs ||
543  cmdline.present("-kt") ||
544  (jet_defs.size() == 0 && !found_unavailable)) {
545  jet_defs.push_back( fj::JetDefinition(fj::kt_algorithm, ktR, strategy));
546  }
547 
548  string filename = cmdline.value<string>("-file", "");
549 
550 
551  if (!cmdline.all_options_used()) {cerr <<
552  "Error: some options were not recognized"<<endl;
553  exit(-1);}
554 
555  for (unsigned idef = 0; idef < jet_defs.size(); idef++) {
556  fj::JetDefinition & jet_def = jet_defs[idef];
557  istream * istr;
558  if (filename == "") istr = &cin;
559  else istr = new ifstream(filename.c_str());
560 
561  for (int iev = 0; iev < nev; iev++) {
562  vector<fj::PseudoJet> jets;
563  vector<fj::PseudoJet> particles;
564  string line;
565  int ndone = 0;
566  while (getline(*istr, line)) {
567  //cout << line<<endl;
568  istringstream linestream(line);
569  if (line == "#END") {
570  ndone += 1;
571  if (ndone == combine) {break;}
572  }
573  if (line.substr(0,1) == "#") {continue;}
574  valarray<double> fourvec(4);
575  if (hydjet) {
576  // special reading from hydjet.txt event record (though actually
577  // this is supposed to be a standard pythia event record, so
578  // being able to read from it is perhaps not so bad an idea...)
579  int ii, istat,id,m1,m2,d1,d2;
580  double mass;
581  linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
582  >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
583  // current file contains mass of particle as 4th entry
584  if (istat == 1) {
585  fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
586  +pow2(fourvec[2])+pow2(mass));
587  }
588  } else {
589  if (massless) {
590  linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
591  fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
592  else {
593  linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
594  }
595  }
596  fj::PseudoJet psjet(fourvec);
597  if (abs(psjet.rap()) < rapmax) {particles.push_back(psjet);}
598  }
599 
600  // add a fake underlying event which is very soft, uniformly distributed
601  // in eta,phi so as to allow one to reconstruct the area that is associated
602  // with each jet.
603  if (add_dense_coverage) {
604  fj::GhostedAreaSpec ghosted_area_spec(ghost_maxrap);
605  //fj::GhostedAreaSpec ghosted_area_spec(-2.0,4.0); // asymmetric range
606  // for plots, reduce the scatter default of 1, to avoid "holes"
607  // in the subsequent calorimeter view
608  ghosted_area_spec.set_grid_scatter(0.5);
609  ghosted_area_spec.add_ghosts(particles);
610  //----- old code ------------------
611  // srand(2);
612  // int nphi = 60;
613  // int neta = 100;
614  // double kt = 1e-1;
615  // for (int iphi = 0; iphi<nphi; iphi++) {
616  // for (int ieta = -neta; ieta<neta+1; ieta++) {
617  // double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
618  // double eta = ieta * (10.0/neta) + rand()*0.001/RAND_MAX;
619  // kt = 1e-20*(1+rand()*0.1/RAND_MAX);
620  // double pminus = kt*exp(-eta);
621  // double pplus = kt*exp(+eta);
622  // double px = kt*sin(phi);
623  // double py = kt*cos(phi);
624  // //cout << kt<<" "<<eta<<" "<<phi<<"\n";
625  // fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
626  // particles.push_back(mom);
627  // }
628  // }
629  }
630 
631  // select the particles that pass the selection cut
632  particles = particles_sel(particles);
633 
634  for (int irepeat = 0; irepeat < repeat ; irepeat++) {
635  int nparticles = particles.size();
636  try {
637  auto_ptr<fj::ClusterSequence> clust_seq;
638  if (do_areas) {
639  clust_seq.reset(new fj::ClusterSequenceArea(particles,jet_def,area_def));
640  } else {
641  clust_seq.reset(new fj::ClusterSequence(particles,jet_def,write));
642  }
643 
644  // repetitive output
645  if (repeatinclkt >= 0.0) {
646  vector<fj::PseudoJet> jets_local = sorted_by_pt(clust_seq->inclusive_jets(repeatinclkt));
647  }
648 
649  if (irepeat != 0) {continue;}
650  cout << "iev "<<iev<< ": number of particles = "<< nparticles << endl;
651  cout << "strategy used = "<< clust_seq->strategy_string()<< endl;
652  if (iev == 0) cout << "Jet Definition: " << jet_def.description() << " (" << fj::fastjet_version_string() << ")" << endl;
653  if (do_areas && iev == 0) cout << "Area definition: " << area_def.description() << endl;
654 
655  // now provide some nice output...
656  if (inclkt >= 0.0) {
657  vector<fj::PseudoJet> jets_local = sorted_by_pt(clust_seq->inclusive_jets(inclkt));
658  print_jets(jets_local, show_constituents);
659 
660  }
661 
662  if (excln > 0) {
663  cout << "Printing "<<excln<<" exclusive jets\n";
664  print_jets(clust_seq->exclusive_jets(excln), show_constituents);
665  }
666 
667  if (excld > 0.0) {
668  cout << "Printing exclusive jets for d = "<<excld<<"\n";
669  print_jets(clust_seq->exclusive_jets(excld), show_constituents);
670  }
671 
672  if (excly > 0.0) {
673  cout << "Printing exclusive jets for ycut = "<<excly<<"\n";
674  print_jets(clust_seq->exclusive_jets_ycut(excly), show_constituents);
675  }
676 
677  if (get_all_dij) {
678  for (int i = nparticles-1; i >= 0; i--) {
679  printf("d for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_dmerge(i));
680  }
681  }
682  if (get_all_yij) {
683  for (int i = nparticles-1; i >= 0; i--) {
684  printf("y for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_ymerge(i));
685  }
686  }
687 
688  // have the option of printing out the subjets (at scale dcut) of
689  // each inclusive jet
690  if (subdcut >= 0.0) {
691  print_jets_and_sub(clust_seq->inclusive_jets(), subdcut);
692  }
693 
694  // useful for testing that recombination sequences are unique
695  if (unique_write) {
696  vector<int> unique_history = clust_seq->unique_history_order();
697  // construct the inverse of the above mapping
698  vector<int> inv_unique_history(clust_seq->history().size());
699  for (unsigned int i = 0; i < unique_history.size(); i++) {
700  inv_unique_history[unique_history[i]] = i;}
701 
702  for (unsigned int i = 0; i < unique_history.size(); i++) {
703  fj::ClusterSequence::history_element el =
704  clust_seq->history()[unique_history[i]];
705  int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
706  int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
707  printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
708  }
709  }
710 
711 
712 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
713  // provide some complementary information for SISCone
714  if (show_cones) {
715  const fj::SISConeExtras * extras =
716  dynamic_cast<const fj::SISConeExtras *>(clust_seq->extras());
717  if (extras == 0)
718  throw fastjet::Error("extras object for SISCone was null (this can happen with certain area types)");
719  cout << "most ambiguous split (difference in squared dist) = "
720  << extras->most_ambiguous_split() << endl;
721  vector<fastjet::PseudoJet> stable_cones(extras->stable_cones());
722  stable_cones = sorted_by_rapidity(stable_cones);
723  for (unsigned int i = 0; i < stable_cones.size(); i++) {
724  //if (stable_cones[i].phi() < 5.0 && stable_cones[i].phi() > 4.0) {
725  printf("%5u %15.8f %15.8f %15.8f\n",
726  i,stable_cones[i].rap(),stable_cones[i].phi(),
727  stable_cones[i].perp() );
728  //}
729  }
730 
731  // also show passes for jets
732  vector<fj::PseudoJet> sisjets = clust_seq->inclusive_jets();
733  printf("\n%15s %15s %15s %12s %8s %8s\n","rap","phi","pt","user-index","pass","nconst");
734  for (unsigned i = 0; i < sisjets.size(); i++) {
735  printf("%15.8f %15.8f %15.8f %12d %8d %8u\n",
736  sisjets[i].rap(), sisjets[i].phi(), sisjets[i].perp(),
737  sisjets[i].user_index(), extras->pass(sisjets[i]),
738  (unsigned int) clust_seq->constituents(sisjets[i]).size()
739  );
740 
741  }
742  }
743 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
744 
745  if (do_bkgd) {
746  double rho, sigma, mean_area, empty_area, n_empty_jets;
748  dynamic_cast<fj::ClusterSequenceAreaBase *>(clust_seq.get());
749  if (do_bkgd_csab) {
750  csab->get_median_rho_and_sigma(bkgd_range, true, rho, sigma, mean_area);
751  empty_area = csab->empty_area(bkgd_range);
752  n_empty_jets = csab->n_empty_jets(bkgd_range);
753  } else if (do_bkgd_jetmedian) {
754  fj::JetMedianBackgroundEstimator bge(bkgd_range);
755  bge.set_provide_fj2_sigma(do_bkgd_fj2);
756  bge.set_cluster_sequence(*csab);
757  rho = bge.rho();
758  sigma = bge.sigma();
759  mean_area = bge.mean_area();
760  empty_area = bge.empty_area();
761  n_empty_jets = bge.n_empty_jets();
762  } else {
763  assert(do_bkgd_gridmedian);
764  double rapmin, rapmax;
765  bkgd_range.get_rapidity_extent(rapmin, rapmax);
766  fj::GridMedianBackgroundEstimator bge(rapmax, 2*ktR);
767  bge.set_particles(particles);
768  rho = bge.rho();
769  sigma = bge.sigma();
770  mean_area = bge.mean_area();
771  empty_area = 0;
772  n_empty_jets = 0;
773  }
774  cout << " rho = " << rho
775  << ", sigma = " << sigma
776  << ", mean_area = " << mean_area
777  << ", empty_area = " << empty_area
778  << ", n_empty_jets = " << n_empty_jets
779  << endl;
780  }
781  } // try
782  catch (fastjet::Error fjerr) {
783  cout << "Caught fastjet error, exiting gracefully" << endl;
784  exit(0);
785  }
786 
787  } // irepeat
788  } // iev
789  // if we've instantiated a plugin, delete it
790  if (jet_def.strategy()==fj::plugin_strategy){
791  delete jet_def.plugin();
792  }
793  // close any file that we've opened
794  if (istr != &cin) delete istr;
795  } // jet_defs
796 
797 }
798 
799 
800 
801 
802 //------ HELPER ROUTINES -----------------------------------------------
803 /// print a single jet
804 void print_jet (const fj::PseudoJet & jet) {
805  unsigned int n_constituents = jet.constituents().size();
806  printf("%15.8f %15.8f %15.8f %8u\n",
807  jet.rap(), jet.phi(), jet.perp(), n_constituents);
808 }
809 
810 
811 //----------------------------------------------------------------------
812 void print_jets(const vector<fj::PseudoJet> & jets_in, bool show_constituents) {
813  vector<fj::PseudoJet> jets;
814  if (ee_print) {
815  jets = sorted_by_E(jets_in);
816  for (unsigned int j = 0; j < jets.size(); j++) {
817  printf("%5u %15.8f %15.8f %15.8f %15.8f\n",
818  j,jets[j].px(),jets[j].py(),jets[j].pz(),jets[j].E());
819  if (show_constituents) {
820  vector<fj::PseudoJet> const_jets = jets[j].constituents();
821  for (unsigned int k = 0; k < const_jets.size(); k++) {
822  printf(" jet%03u %15.8f %15.8f %15.8f %15.8f\n",j,const_jets[k].px(),
823  const_jets[k].py(),const_jets[k].pz(),const_jets[k].E());
824  }
825  cout << "\n\n";
826  }
827 
828  }
829  } else {
830  jets = sorted_by_pt(jets_in);
831  for (unsigned int j = 0; j < jets.size(); j++) {
832  printf("%5u %15.8f %15.8f %15.8f",
833  j,jets[j].rap(),jets[j].phi(),jets[j].perp());
834  // also print out the scalar area and the perp component of the
835  // 4-vector (just enough to check a reasonable 4-vector?)
836  if (do_areas) printf(" %15.8f %15.8f", jets[j].area(),
837  jets[j].area_4vector().perp());
838  cout << "\n";
839 
840  if (show_constituents) {
841  vector<fj::PseudoJet> const_jets = jets[j].constituents();
842  for (unsigned int k = 0; k < const_jets.size(); k++) {
843  printf(" jet%03u %15.8f %15.8f %15.8f %5d\n",j,const_jets[k].rap(),
844  const_jets[k].phi(),sqrt(const_jets[k].kt2()), const_jets[k].cluster_hist_index());
845  }
846  cout << "\n\n";
847  }
848  }
849  }
850 
851  if (rootfile != "") {
852  ofstream ostr(rootfile.c_str());
853  ostr << "# " << cmdline_p->command_line() << endl;
854  ostr << "# output for root" << endl;
855  assert(jets.size() > 0);
856  jets[0].validated_cs()->print_jets_for_root(jets,ostr);
857  }
858 
859 }
860 
861 
862 //----- SUBJETS --------------------------------------------------------
863 /// a function that pretty prints a list of jets and the subjets for each
864 /// one
865 void print_jets_and_sub (const vector<fj::PseudoJet> & jets, double dcut) {
866 
867  // sort jets into increasing pt
868  vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);
869 
870  // label the columns
871  printf("Printing jets and their subjets with subdcut = %10.5f\n",dcut);
872  printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity",
873  "phi", "pt", "n constituents");
874 
875  // have various kinds of subjet finding, to test consistency among them
876  enum SubType {internal, newclust_dcut, newclust_R};
877  SubType subtype = internal;
878  //SubType subtype = newclust_dcut;
879  //SubType subtype = newclust_R;
880 
881  // print out the details for each jet
882  //for (unsigned int i = 0; i < sorted_jets.size(); i++) {
883  for (vector<fj::PseudoJet>::const_iterator jet = sorted_jets.begin();
884  jet != sorted_jets.end(); jet++) {
885  const fj::JetDefinition & jet_def = jet->validated_cs()->jet_def();
886 
887  // if jet pt^2 < dcut with kt alg, then some methods of
888  // getting subjets will return nothing -- so skip the jet
889  if (jet_def.jet_algorithm() == fj::kt_algorithm
890  && jet->perp2() < dcut) continue;
891 
892 
893  printf("%5u ",(unsigned int) (jet - sorted_jets.begin()));
894  print_jet(*jet);
895  vector<fj::PseudoJet> subjets;
896  fj::ClusterSequence * cspoint;
897  if (subtype == internal) {
898  cspoint = 0;
899  subjets = jet->exclusive_subjets(dcut);
900  double ddnp1 = jet->exclusive_subdmerge_max(subjets.size());
901  double ddn = jet->exclusive_subdmerge_max(subjets.size()-1);
902  cout << " for " << ddnp1 << " < d < " << ddn << " one has " << endl;
903  } else if (subtype == newclust_dcut) {
904  cspoint = new fj::ClusterSequence(jet->constituents(), jet_def);
905  subjets = cspoint->exclusive_jets(dcut);
906  } else if (subtype == newclust_R) {
907  assert(jet_def.jet_algorithm() == fj::cambridge_algorithm);
908  fj::JetDefinition subjd(jet_def.jet_algorithm(), jet_def.R()*sqrt(dcut));
909  cspoint = new fj::ClusterSequence(jet->constituents(), subjd);
910  subjets = cspoint->inclusive_jets();
911  } else {
912  cerr << "unrecognized subtype for subjet finding" << endl;
913  exit(-1);
914  }
915 
916  subjets = sorted_by_pt(subjets);
917  for (unsigned int j = 0; j < subjets.size(); j++) {
918  printf(" -sub-%02u ",j);
919  print_jet(subjets[j]);
920  }
921 
922  if (cspoint != 0) delete cspoint;
923 
924  //fj::ClusterSequence subseq(clust_seq->constituents(sorted_jets[i]),
925  // fj::JetDefinition(fj::cambridge_algorithm, 0.4));
926  //vector<fj::PseudoJet> subjets = sorted_by_pt(subseq.inclusive_jets());
927  //for (unsigned int j = 0; j < subjets.size(); j++) {
928  // printf(" -sub-%02u ",j);
929  // print_jet(subseq, subjets[j]);
930  //}
931  }
932 
933 }
934 
const JetDefinition & jet_def() const
return a reference to the jet definition
summing the 4-momenta
const Plugin * plugin() const
return a pointer to the plugin
JetAlgorithm jet_algorithm() const
return information about the definition...
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:785
deals with clustering
Class to estimate the pt density of the background per unit area, using the median of the distributio...
General class for user to obtain ClusterSequence with additional area information.
Selector SelectorNHardest(unsigned int n)
select the n hardest objects
Definition: Selector.cc:1022
const ClusterSequence * validated_cs() const
shorthand for validated_cluster_sequence()
Definition: PseudoJet.cc:443
const std::vector< PseudoJet > & stable_cones() const
returns a reference to the vector of stable cones (aka protocones)
std::vector< PseudoJet > exclusive_subjets(const double &dcut) const
return a vector of all subjets of the current jet (in the sense of the exclusive algorithm) that woul...
Definition: PseudoJet.cc:589
static void print_banner()
This is the function that is automatically called during clustering to print the FastJet banner...
Selector SelectorAbsRapMax(double absrapmax)
select objects with |rap| <= absrapmax
Definition: Selector.cc:828
Implementation of the D0 Run II Cone (plugin for fastjet v2.1 upwards)
the longitudinally invariant kt algorithm
vector< PseudoJet > sorted_by_rapidity(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing rapidity
Definition: PseudoJet.cc:793
class that holds a generic area definition
virtual double n_empty_jets(const Selector &selector) const
return something similar to the number of pure ghost jets in the given selector&#39;s range in an active ...
Implementation of the spherical version of the SISCone algorithm (plugin for fastjet v2...
void get_rapidity_extent(double &rapmin, double &rapmax) const
returns the rapidity range for which it may return "true"
Definition: Selector.hh:213
double most_ambiguous_split() const
return the smallest difference in squared distance encountered during splitting between a particle an...
the plugin has been used...
Implementation of the TrackJet algorithm (plugin for fastjet v2.4 upwards)
std::vector< PseudoJet > exclusive_jets(const double &dcut) const
return a vector of all jets (in the sense of the exclusive algorithm) that would be obtained when run...
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition: PseudoJet.cc:566
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:141
int main()
an example program showing how to use fastjet
Definition: 01-basic.cc:50
Specification for the computation of the Voronoi jet area.
the FastJet namespace
base class that sets interface for extensions of ClusterSequence that provide information about the a...
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
Definition: PseudoJet.cc:801
Implementation of the e+e- Cambridge algorithm (plugin for fastjet v2.4 upwards)
void print_jets(const vector< fastjet::PseudoJet > &)
a function that pretty prints a list of jets
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:41
RecombinationScheme
the various recombination schemes
string fastjet_version_string()
return a string containing information about the release
double perp2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:139
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:106
double exclusive_subdmerge_max(int nsub) const
return the maximum dij that occurred in the whole event at the stage that the nsub+1 -> nsub merge of...
Definition: PseudoJet.cc:648
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:121
Implementation of the PxCone algorithm (plugin for fastjet v2.1 upwards)
Definition: PxConePlugin.hh:66
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:141
virtual void get_median_rho_and_sigma(const Selector &selector, bool use_area_4vector, double &median, double &sigma, double &mean_area) const
using jets withing the selector range (and with 4-vector areas if use_area_4vector), calculate the median pt/area, as well as an "error" (uncertainty), which is defined as the 1-sigma half-width of the distribution of pt/A, obtained by looking for the point below which we have (1-0.6827)/2 of the jets (including empty jets).
virtual double empty_area(const Selector &selector) const
return the total area, corresponding to the given Selector, that is free of jets, in general based on...
Implementation of the CMS Iterative Cone (plugin for fastjet v2.4 upwards)
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
std::vector< PseudoJet > inclusive_jets(const double &ptmin=0.0) const
return a vector of all jets (in the sense of the inclusive algorithm) with pt >= ptmin.
void set_fj2_placement(bool val)
if val is true, set ghost placement as it was in FastJet 2.X.
Implementation of the e+e- Jade algorithm (plugin for fastjet v2.4 upwards)
Definition: JadePlugin.hh:75
std::string description() const
return a textual description of the current jet definition
Class that provides extra information about a SISCone clustering.
Parameters to configure the computation of jet areas using ghosts.
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
Implementation of the SISCone algorithm (plugin for fastjet v2.1 upwards)
plugin for fastjet (v3.0 upwards) that clusters particles such that all particles in a given cell of ...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:65
class that is intended to hold a full definition of the jet clusterer
int pass(const PseudoJet &jet) const
return the # of the pass at which a given jet was found; will return -1 if the pass is invalid ...