FastJet  3.1.0-beta.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
fastjet_timing_plugins.cc
1 //STARTHEADER
2 // $Id: fastjet_timing_plugins.cc 3473 2014-07-29 09:45:34Z soyez $
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 /// -compare-strategy STRAT
126 /// compares the output of the default strategy (possibly as specified
127 /// with -strategy) with that from STRAT. Currently compares the history.
128 
129 /// Algorithms
130 /// ----------
131 /// -all-algs runs all algorithms
132 ///
133 /// -kt switch to the longitudinally invariant kt algorithm
134 /// Note: this is the default one.
135 ///
136 /// -cam switch to the inclusive Cambridge/Aachen algorithm --
137 /// note that the option -excld dcut provides a clustering
138 /// up to the dcut which is the minimum squared
139 /// distance between any pair of jets.
140 ///
141 /// -antikt switch to the anti-kt clustering algorithm
142 ///
143 /// -genkt switch to the genkt algorithm
144 /// you can provide the parameter of the alg as an argument to
145 /// -genkt (1 by default)
146 ///
147 /// -eekt switch to the e+e- kt algorithm
148 ///
149 /// -eegenkt switch to the genkt algorithm
150 /// you can provide the parameter of the alg as an argument to
151 /// -ee_genkt (1 by default)
152 ///
153 /// plugins (don't delete this line)
154 ///
155 /// -pxcone switch to the PxCone jet algorithm
156 ///
157 /// -siscone switch to the SISCone jet algorithm (seedless cones)
158 /// -sisconespheri switch to the Spherical SISCone jet algorithm (seedless cones)
159 ///
160 /// -midpoint switch to CDF's midpoint code
161 /// -jetclu switch to CDF's jetclu code
162 ///
163 /// -d0runipre96cone switch to the D0RunIpre96Cone plugin
164 /// -d0runicone switch to the D0RunICone plugin
165 ///
166 /// -d0runiicone switch to D0's run II midpoint cone
167 ///
168 /// -trackjet switch to the TrackJet plugin
169 ///
170 /// -atlascone switch to the ATLASCone plugin
171 ///
172 /// -eecambridge switch to the EECambridge plugin
173 ///
174 /// -jade switch to the Jade plugin
175 ///
176 /// -cmsiterativecone switch to the CMSIterativeCone plugin
177 ///
178 /// -gridjet switch to the GridJet plugin
179 ///
180 /// end of plugins (don't delete this line)
181 ///
182 ///
183 /// Options for running algs
184 /// ------------------------
185 ///
186 /// -r sets the radius of the jet algorithm (default = 1.0)
187 ///
188 /// -overlap | -f sets the overlap fraction in cone algs with split-merge
189 ///
190 /// -seed sets the seed threshold
191 ///
192 /// -strategy N indicate stratgey from the enum fastjet::Strategy (see
193 /// fastjet/JetDefinition.hh).
194 ///
195 
196 
197 #include "fastjet/ClusterSequenceArea.hh"
198 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
199 #include "fastjet/tools/GridMedianBackgroundEstimator.hh"
200 #include "fastjet/Selector.hh"
201 #include<iostream>
202 #include<iomanip>
203 #include<sstream>
204 #include<fstream>
205 #include<valarray>
206 #include<vector>
207 #include <cstdlib>
208 //#include<cstddef> // for size_t
209 #include "CmdLine.hh"
210 
211 // get info on how fastjet was configured
212 #include "fastjet/config.h"
213 
214 // include the installed plugins (don't delete this line)
215 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
216 #include "fastjet/SISConePlugin.hh"
217 #include "fastjet/SISConeSphericalPlugin.hh"
218 #endif
219 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
220 #include "fastjet/CDFMidPointPlugin.hh"
221 #include "fastjet/CDFJetCluPlugin.hh"
222 #endif
223 #ifdef FASTJET_ENABLE_PLUGIN_PXCONE
224 #include "fastjet/PxConePlugin.hh"
225 #endif
226 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
227 #include "fastjet/D0RunIIConePlugin.hh"
228 #endif
229 #ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
230 #include "fastjet/TrackJetPlugin.hh"
231 #endif
232 #ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
233 #include "fastjet/ATLASConePlugin.hh"
234 #endif
235 #ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
236 #include "fastjet/EECambridgePlugin.hh"
237 #endif
238 #ifdef FASTJET_ENABLE_PLUGIN_JADE
239 #include "fastjet/JadePlugin.hh"
240 #endif
241 #ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
242 #include "fastjet/CMSIterativeConePlugin.hh"
243 #endif
244 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
245 #include "fastjet/D0RunIpre96ConePlugin.hh"
246 #include "fastjet/D0RunIConePlugin.hh"
247 #endif
248 #ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
249 #include "fastjet/GridJetPlugin.hh"
250 #endif
251 // end of installed plugins inclusion (don't delete this line)
252 
253 using namespace std;
254 
255 // to avoid excessive typing, use the fastjet namespace
256 using namespace fastjet;
257 
258 inline double pow2(const double x) {return x*x;}
259 
260 // pretty print the jets and their subjets
261 void print_jets_and_sub (const vector<PseudoJet> & jets, double dcut);
262 
263 // have various kinds of subjet finding, to test consistency among them
264 //
265 // this is needed in print_jets_and_sub and declaring it in the
266 // function scope results in errors with older intel compilers (due to
267 // the overloaded == operator in PseudoJet which results in the "a
268 // template argument may not reference a local type" error)
269 enum SubType {subtype_internal, subtype_newclust_dcut, subtype_newclust_R};
270 
271 void do_compare_strategy(int iev,
272  const vector<PseudoJet> & particles,
273  const JetDefinition & jet_def,
274  const ClusterSequence & cs,
275  int compare_strategy);
276 
277 
278 string rootfile;
279 CmdLine * cmdline_p;
280 
281 bool do_areas;
282 
283 /// sort and pretty print jets, with exact behaviour depending on
284 /// whether ee_print is true or not
285 bool ee_print = false;
286 void print_jets(const vector<PseudoJet> & jets, bool show_const = false);
287 
288 bool found_unavailable = false;
289 void is_unavailable(const string & algname) {
290  cerr << algname << " requested, but not available for this compilation" << endl;
291  found_unavailable = true;
292  //exit(0);
293 }
294 
295 
296 /// a program to test and time a range of algorithms as implemented or
297 /// wrapped in fastjet
298 int main (int argc, char ** argv) {
299 
300  ClusterSequence::print_banner();
301 
302  CmdLine cmdline(argc,argv);
303  cmdline_p = &cmdline;
304  // allow the use to specify the Strategy either through the
305  // -clever or the -strategy options (both will take numerical
306  // values); the latter will override the former.
307  Strategy strategy = Strategy(cmdline.int_val("-strategy",
308  cmdline.int_val("-clever", Best)));
309  int repeat = cmdline.int_val("-repeat",1);
310  int combine = cmdline.int_val("-combine",1);
311  bool write = cmdline.present("-write");
312  bool unique_write = cmdline.present("-unique_write");
313  bool hydjet = cmdline.present("-hydjet");
314  double ktR = cmdline.double_val("-r",1.0);
315  ktR = cmdline.double_val("-R",ktR); // allow -r and -R
316  double inclkt = cmdline.double_val("-incl",-1.0);
317  double repeatinclkt = cmdline.double_val("-repeat-incl",-1.0);
318  int excln = cmdline.int_val ("-excln",-1);
319  double excld = cmdline.double_val("-excld",-1.0);
320  double excly = cmdline.double_val("-excly",-1.0);
321  ee_print = cmdline.present("-ee-print");
322  bool get_all_dij = cmdline.present("-get-all-dij");
323  bool get_all_yij = cmdline.present("-get-all-yij");
324  double subdcut = cmdline.double_val("-subdcut",-1.0);
325  double rapmax = cmdline.double_val("-rapmax",1.0e305);
326  if (cmdline.present("-etamax")) {
327  cerr << "WARNING: -etamax options actually sets maximum rapidity (and overrides -rapmax)\n";
328  rapmax = cmdline.double_val("-etamax");
329  }
330  bool show_constituents = cmdline.present("-const");
331  bool massless = cmdline.present("-massless");
332  int nev = cmdline.int_val("-nev",1);
333  int skip = cmdline.int_val("-skip",0);
334  bool add_dense_coverage = cmdline.present("-dense");
335  double ghost_maxrap = cmdline.value("-ghost-maxrap",5.0);
336  bool all_algs = cmdline.present("-all-algs");
337  // have the option of comparing the clustering results to those
338  // obtained with a different clustering strategy; the misuse the
339  // "plugin_strategy" to indicate that no comparison is needed.
340  // Does not currently support areas
341  int compare_strategy = cmdline.value<int>("-compare-strategy", plugin_strategy);
342 
343 
344  Selector particles_sel = (cmdline.present("-nhardest"))
345  ? SelectorNHardest(cmdline.value<unsigned int>("-nhardest"))
346  : SelectorIdentity();
347 
348  do_areas = cmdline.present("-area");
349  AreaDefinition area_def;
350  if (do_areas) {
351  assert(!write); // it's incompatible
352  GhostedAreaSpec ghost_spec(ghost_maxrap,
353  cmdline.value("-area:repeat", 1),
354  cmdline.value("-ghost-area", 0.01));
355  if (cmdline.present("-area:fj2")) ghost_spec.set_fj2_placement(true);
356  if (cmdline.present("-area:explicit")) {
357  area_def = AreaDefinition(active_area_explicit_ghosts, ghost_spec);
358  } else if (cmdline.present("-area:passive")) {
359  area_def = AreaDefinition(passive_area, ghost_spec);
360  } else if (cmdline.present("-area:voronoi")) {
361  double Rfact = cmdline.value<double>("-area:voronoi");
362  area_def = AreaDefinition(voronoi_area,
363  VoronoiAreaSpec(Rfact));
364  } else {
365  cmdline.present("-area:active"); // allow, but do not require, arg
366  area_def = AreaDefinition(active_area, ghost_spec);
367  }
368  }
369  bool do_bkgd = cmdline.present("-bkgd"); // background estimation
370  bool do_bkgd_csab = false, do_bkgd_jetmedian = false, do_bkgd_fj2 = false;
371  bool do_bkgd_gridmedian = false;
372  Selector bkgd_range;
373  if (do_bkgd) {
374  bkgd_range = SelectorAbsRapMax(ghost_maxrap - ktR);
375  if (cmdline.present("-bkgd:csab")) {do_bkgd_csab = true;}
376  else if (cmdline.present("-bkgd:jetmedian")) {do_bkgd_jetmedian = true;
377  do_bkgd_fj2 = cmdline.present("-bkgd:fj2");
378  } else if (cmdline.present("-bkgd:gridmedian")) {do_bkgd_gridmedian = true;
379  } else {
380  throw Error("with the -bkgd option, some particular background must be specified (csab or jetmedian)");
381  }
382  assert(do_areas || do_bkgd_gridmedian);
383  }
384 
385  bool show_cones = cmdline.present("-cones"); // only works for siscone
386 
387  // for cone algorithms
388  // allow -f and -overlap
389  double overlap_threshold = cmdline.double_val("-overlap",0.5);
390  overlap_threshold = cmdline.double_val("-f",overlap_threshold);
391  double seed_threshold = cmdline.double_val("-seed",1.0);
392 
393  // for ee algorithms, allow to specify ycut
394  double ycut = cmdline.double_val("-ycut",0.08);
395 
396  // for printing jets to a file for reading by root
397  rootfile = cmdline.value<string>("-root","");
398 
399  // out default scheme is the E_scheme
400  RecombinationScheme scheme = E_scheme;
401 
402  // The following option causes the Cambridge algo to be used.
403  // Note that currently the only output that works sensibly here is
404  // "-incl 0"
405  vector<JetDefinition> jet_defs;
406  if (all_algs || cmdline.present("-cam") || cmdline.present("-CA")) {
407  jet_defs.push_back( JetDefinition(cambridge_algorithm, ktR, scheme, strategy));
408  }
409  if (all_algs || cmdline.present("-antikt")) {
410  jet_defs.push_back( JetDefinition(antikt_algorithm, ktR, scheme, strategy));
411  }
412  if (all_algs || cmdline.present("-genkt")) {
413  double p;
414  if (cmdline.present("-genkt")) p = cmdline.value<double>("-genkt");
415  else p = -0.5;
416  jet_defs.push_back( JetDefinition(genkt_algorithm, ktR, p, scheme, strategy));
417  }
418  if (all_algs || cmdline.present("-eekt")) {
419  jet_defs.push_back( JetDefinition(ee_kt_algorithm));
420  }
421  if (all_algs || cmdline.present("-eegenkt")) {
422  double p;
423  if (cmdline.present("-eegenkt")) p = cmdline.value<double>("-eegenkt");
424  else p = -0.5;
425  jet_defs.push_back( JetDefinition(ee_genkt_algorithm, ktR, p, scheme, strategy));
426 
427 // checking if one asks to run a plugin (don't delete this line)
428  }
429  if (all_algs || cmdline.present("-midpoint")) {
430 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
431  typedef CDFMidPointPlugin MPPlug; // for brevity
432  double cone_area_fraction = 1.0;
433  int max_pair_size = 2;
434  int max_iterations = 100;
435  MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
436  if (cmdline.present("-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
437  if (cmdline.present("-sm-pt")) sm_scale = MPPlug::SM_pt; // default
438  if (cmdline.present("-sm-mt")) sm_scale = MPPlug::SM_mt;
439  if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et;
440  jet_defs.push_back( JetDefinition( new CDFMidPointPlugin (
441  seed_threshold, ktR,
442  cone_area_fraction, max_pair_size,
443  max_iterations, overlap_threshold,
444  sm_scale)));
445 #else // FASTJET_ENABLE_PLUGIN_CDFCONES
446  is_unavailable("MidPoint");
447 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
448  }
449  if (all_algs || cmdline.present("-pxcone")) {
450 #ifdef FASTJET_ENABLE_PLUGIN_PXCONE
451  double min_jet_energy = 5.0;
452  jet_defs.push_back( JetDefinition( new PxConePlugin (
453  ktR, min_jet_energy,
454  overlap_threshold)));
455 #else // FASTJET_ENABLE_PLUGIN_PXCONE
456  is_unavailable("PxCone");
457 #endif // FASTJET_ENABLE_PLUGIN_PXCONE
458  }
459  if (all_algs || cmdline.present("-jetclu")) {
460 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
461  jet_defs.push_back( JetDefinition( new CDFJetCluPlugin (
462  ktR, overlap_threshold, seed_threshold)));
463 #else // FASTJET_ENABLE_PLUGIN_CDFCONES
464  is_unavailable("JetClu");
465 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
466  }
467  if (all_algs || cmdline.present("-siscone") || cmdline.present("-sisconespheri")) {
468 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
469  typedef SISConePlugin SISPlug; // for brevity
470  int npass = cmdline.value("-npass",0);
471  if (all_algs || cmdline.present("-siscone")) {
472  double sisptmin = cmdline.value("-sisptmin",0.0);
473  SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass,sisptmin);
474  if (cmdline.present("-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
475  if (cmdline.present("-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
476  if (cmdline.present("-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
477  if (cmdline.present("-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
478  // cause it to use the jet-definition's own recombiner
479  plugin->set_use_jet_def_recombiner(true);
480  jet_defs.push_back( JetDefinition(plugin));
481  }
482  if (all_algs || cmdline.present("-sisconespheri")) {
483  double sisEmin = cmdline.value("-sisEmin",0.0);
484  SISConeSphericalPlugin * plugin =
485  new SISConeSphericalPlugin(ktR, overlap_threshold,npass,sisEmin);
486  if (cmdline.present("-ghost-sep")) {
487  plugin->set_ghost_separation_scale(cmdline.value<double>("-ghost-sep"));
488  }
489  jet_defs.push_back( JetDefinition(plugin));
490  }
491 #else // FASTJET_ENABLE_PLUGIN_SISCONE
492  is_unavailable("SISCone");
493 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
494  }
495  if (all_algs || cmdline.present("-d0runiicone")) {
496 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
497  double min_jet_Et = 6.0; // was 8 GeV in earlier work
498  jet_defs.push_back( JetDefinition(new D0RunIIConePlugin(ktR,min_jet_Et)));
499 #else // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
500  is_unavailable("D0RunIICone");
501 #endif // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
502  }
503  if (all_algs || cmdline.present("-trackjet")) {
504 #ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
505  jet_defs.push_back( JetDefinition(new TrackJetPlugin(ktR)));
506 #else // FASTJET_ENABLE_PLUGIN_TRACKJET
507  is_unavailable("TrackJet");
508 #endif // FASTJET_ENABLE_PLUGIN_TRACKJET
509  }
510  if (all_algs || cmdline.present("-atlascone")) {
511 #ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
512  jet_defs.push_back( JetDefinition(new ATLASConePlugin(ktR)));
513 #else // FASTJET_ENABLE_PLUGIN_ATLASCONE
514  is_unavailable("ATLASCone");
515 #endif // FASTJET_ENABLE_PLUGIN_ATLASCONE
516  }
517  if (all_algs || cmdline.present("-eecambridge")) {
518 #ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
519  jet_defs.push_back( JetDefinition(new EECambridgePlugin(ycut)));
520 #else // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
521  is_unavailable("EECambridge");
522 #endif // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
523  }
524  if (all_algs || cmdline.present("-jade")) {
525 #ifdef FASTJET_ENABLE_PLUGIN_JADE
526  jet_defs.push_back( JetDefinition(new JadePlugin()));
527 #else // FASTJET_ENABLE_PLUGIN_JADE
528  is_unavailable("Jade");
529 #endif // FASTJET_ENABLE_PLUGIN_JADE
530  }
531  if (all_algs || cmdline.present("-cmsiterativecone")) {
532 #ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
533  jet_defs.push_back( JetDefinition(new CMSIterativeConePlugin(ktR,seed_threshold)));
534 #else // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
535  is_unavailable("CMSIterativeCone");
536 #endif // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
537  }
538  if (all_algs || cmdline.present("-d0runipre96cone")) {
539 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
540  jet_defs.push_back( JetDefinition(new D0RunIpre96ConePlugin(ktR, seed_threshold, overlap_threshold)));
541 #else // FASTJET_ENABLE_PLUGIN_D0RUNICONE
542  is_unavailable("D0RunIpre96Cone");
543 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
544  }
545  if (all_algs || cmdline.present("-d0runicone")) {
546 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
547  jet_defs.push_back( JetDefinition(new D0RunIConePlugin(ktR, seed_threshold, overlap_threshold)));
548 #else // FASTJET_ENABLE_PLUGIN_D0RUNICONE
549  is_unavailable("D0RunICone");
550 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
551  }
552  if (all_algs || cmdline.present("-gridjet")) {
553 #ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
554  // we want a grid_ymax of 5.0, but when using R=0.4 (i.e. grid
555  // spacing of 0.8), this leads to 12.5 grid cells; depending on
556  // whether this is 12.499999999999 or 12.5000000....1 this gets
557  // converted either to 12 or 13, making the results sensitive to
558  // rounding errors.
559  //
560  // Instead we therefore take 4.9999999999, which avoids this problem.
561  double grid_ymax = 4.9999999999;
562  jet_defs.push_back( JetDefinition(new GridJetPlugin(grid_ymax, ktR*2.0)));
563 #else // FASTJET_ENABLE_PLUGIN_GRIDJET
564  is_unavailable("GridJet");
565 #endif // FASTJET_ENABLE_PLUGIN_GRIDJET
566 // end of checking if one asks to run a plugin (don't delete this line)
567  }
568  if (all_algs ||
569  cmdline.present("-kt") ||
570  (jet_defs.size() == 0 && !found_unavailable)) {
571  jet_defs.push_back( JetDefinition(kt_algorithm, ktR, strategy));
572  }
573 
574  string filename = cmdline.value<string>("-file", "");
575 
576 
577  if (!cmdline.all_options_used()) {cerr <<
578  "Error: some options were not recognized"<<endl;
579  exit(-1);}
580 
581  for (unsigned idef = 0; idef < jet_defs.size(); idef++) {
582  JetDefinition & jet_def = jet_defs[idef];
583  istream * istr;
584  if (filename == "") istr = &cin;
585  else istr = new ifstream(filename.c_str());
586 
587  for (int iev = 0; iev < nev; iev++) {
588  vector<PseudoJet> jets;
589  vector<PseudoJet> particles;
590  string line;
591  int ndone = 0;
592  while (getline(*istr, line)) {
593  //cout << line<<endl;
594  istringstream linestream(line);
595  if (line == "#END") {
596  ndone += 1;
597  if (ndone == combine) {break;}
598  }
599  if (line.substr(0,1) == "#") {continue;}
600  valarray<double> fourvec(4);
601  if (hydjet) {
602  // special reading from hydjet.txt event record (though actually
603  // this is supposed to be a standard pythia event record, so
604  // being able to read from it is perhaps not so bad an idea...)
605  int ii, istat,id,m1,m2,d1,d2;
606  double mass;
607  linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
608  >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
609  // current file contains mass of particle as 4th entry
610  if (istat == 1) {
611  fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
612  +pow2(fourvec[2])+pow2(mass));
613  }
614  } else {
615  if (massless) {
616  linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
617  fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
618  else {
619  linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
620  }
621  }
622  PseudoJet psjet(fourvec);
623  if (abs(psjet.rap()) < rapmax) {particles.push_back(psjet);}
624  }
625 
626  // add a fake underlying event which is very soft, uniformly distributed
627  // in eta,phi so as to allow one to reconstruct the area that is associated
628  // with each jet.
629  if (add_dense_coverage) {
630  GhostedAreaSpec ghosted_area_spec(ghost_maxrap);
631  //GhostedAreaSpec ghosted_area_spec(-2.0,4.0); // asymmetric range
632  // for plots, reduce the scatter default of 1, to avoid "holes"
633  // in the subsequent calorimeter view
634  ghosted_area_spec.set_grid_scatter(0.5);
635  ghosted_area_spec.add_ghosts(particles);
636  //----- old code ------------------
637  // srand(2);
638  // int nphi = 60;
639  // int neta = 100;
640  // double kt = 1e-1;
641  // for (int iphi = 0; iphi<nphi; iphi++) {
642  // for (int ieta = -neta; ieta<neta+1; ieta++) {
643  // double phi = (iphi+0.5) * (twopi/nphi) + rand()*0.001/RAND_MAX;
644  // double eta = ieta * (10.0/neta) + rand()*0.001/RAND_MAX;
645  // kt = 1e-20*(1+rand()*0.1/RAND_MAX);
646  // double pminus = kt*exp(-eta);
647  // double pplus = kt*exp(+eta);
648  // double px = kt*sin(phi);
649  // double py = kt*cos(phi);
650  // //cout << kt<<" "<<eta<<" "<<phi<<"\n";
651  // PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
652  // particles.push_back(mom);
653  // }
654  // }
655  }
656 
657  // select the particles that pass the selection cut
658  particles = particles_sel(particles);
659 
660  // allow user to skip some number of events (e.g. for easier bug-chasing)
661  if (iev < skip) continue;
662 
663  for (int irepeat = 0; irepeat < repeat ; irepeat++) {
664  int nparticles = particles.size();
665  try {
666  auto_ptr<ClusterSequence> clust_seq;
667  if (do_areas) {
668  clust_seq.reset(new ClusterSequenceArea(particles,jet_def,area_def));
669  } else {
670  clust_seq.reset(new ClusterSequence(particles,jet_def,write));
671  }
672  if (compare_strategy != plugin_strategy) {
673  do_compare_strategy(iev, particles, jet_def, *clust_seq, compare_strategy);
674  }
675 
676  // repetitive output
677  if (repeatinclkt >= 0.0) {
678  vector<PseudoJet> jets_local = sorted_by_pt(clust_seq->inclusive_jets(repeatinclkt));
679  }
680 
681  if (irepeat != 0) {continue;}
682  cout << "iev "<<iev<< ": number of particles = "<< nparticles << endl;
683  cout << "strategy used = "<< clust_seq->strategy_string()<< endl;
684  if (iev == 0) cout << "Jet Definition: " << jet_def.description() << " (" << fastjet_version_string() << ")" << endl;
685  if (do_areas && iev == 0) cout << "Area definition: " << area_def.description() << endl;
686 
687  // now provide some nice output...
688  if (inclkt >= 0.0) {
689  vector<PseudoJet> jets_local = sorted_by_pt(clust_seq->inclusive_jets(inclkt));
690  print_jets(jets_local, show_constituents);
691 
692  }
693 
694  if (excln > 0) {
695  cout << "Printing "<<excln<<" exclusive jets\n";
696  print_jets(clust_seq->exclusive_jets(excln), show_constituents);
697  }
698 
699  if (excld > 0.0) {
700  cout << "Printing exclusive jets for d = "<<excld<<"\n";
701  print_jets(clust_seq->exclusive_jets(excld), show_constituents);
702  }
703 
704  if (excly > 0.0) {
705  cout << "Printing exclusive jets for ycut = "<<excly<<"\n";
706  print_jets(clust_seq->exclusive_jets_ycut(excly), show_constituents);
707  }
708 
709  if (get_all_dij) {
710  for (int i = nparticles-1; i >= 0; i--) {
711  printf("d for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_dmerge(i));
712  }
713  }
714  if (get_all_yij) {
715  for (int i = nparticles-1; i >= 0; i--) {
716  printf("y for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_ymerge(i));
717  }
718  }
719 
720  // have the option of printing out the subjets (at scale dcut) of
721  // each inclusive jet
722  if (subdcut >= 0.0) {
723  print_jets_and_sub(clust_seq->inclusive_jets(), subdcut);
724  }
725 
726  // useful for testing that recombination sequences are unique
727  if (unique_write) {
728  vector<int> unique_history = clust_seq->unique_history_order();
729  // construct the inverse of the above mapping
730  vector<int> inv_unique_history(clust_seq->history().size());
731  for (unsigned int i = 0; i < unique_history.size(); i++) {
732  inv_unique_history[unique_history[i]] = i;}
733 
734  for (unsigned int i = 0; i < unique_history.size(); i++) {
735  ClusterSequence::history_element el =
736  clust_seq->history()[unique_history[i]];
737  int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
738  int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
739  printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
740  }
741  }
742 
743 
744 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
745  // provide some complementary information for SISCone
746  if (show_cones) {
747  const SISConeExtras * extras =
748  dynamic_cast<const SISConeExtras *>(clust_seq->extras());
749  if (extras == 0)
750  throw fastjet::Error("extras object for SISCone was null (this can happen with certain area types)");
751  cout << "most ambiguous split (difference in squared dist) = "
752  << extras->most_ambiguous_split() << endl;
753  vector<fastjet::PseudoJet> stable_cones(extras->stable_cones());
754  stable_cones = sorted_by_rapidity(stable_cones);
755  for (unsigned int i = 0; i < stable_cones.size(); i++) {
756  //if (stable_cones[i].phi() < 5.0 && stable_cones[i].phi() > 4.0) {
757  printf("%5u %15.8f %15.8f %15.8f\n",
758  i,stable_cones[i].rap(),stable_cones[i].phi(),
759  stable_cones[i].perp() );
760  //}
761  }
762 
763  // also show passes for jets
764  vector<PseudoJet> sisjets = clust_seq->inclusive_jets();
765  printf("\n%15s %15s %15s %12s %8s %8s\n","rap","phi","pt","user-index","pass","nconst");
766  for (unsigned i = 0; i < sisjets.size(); i++) {
767  printf("%15.8f %15.8f %15.8f %12d %8d %8u\n",
768  sisjets[i].rap(), sisjets[i].phi(), sisjets[i].perp(),
769  sisjets[i].user_index(), extras->pass(sisjets[i]),
770  (unsigned int) clust_seq->constituents(sisjets[i]).size()
771  );
772 
773  }
774  }
775 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
776 
777  if (do_bkgd) {
778  double rho, sigma, mean_area, empty_area, n_empty_jets;
779  ClusterSequenceAreaBase * csab =
780  dynamic_cast<ClusterSequenceAreaBase *>(clust_seq.get());
781  if (do_bkgd_csab) {
782  csab->get_median_rho_and_sigma(bkgd_range, true, rho, sigma, mean_area);
783  empty_area = csab->empty_area(bkgd_range);
784  n_empty_jets = csab->n_empty_jets(bkgd_range);
785  } else if (do_bkgd_jetmedian) {
786  JetMedianBackgroundEstimator bge(bkgd_range);
787  bge.set_provide_fj2_sigma(do_bkgd_fj2);
788  bge.set_cluster_sequence(*csab);
789  rho = bge.rho();
790  sigma = bge.sigma();
791  mean_area = bge.mean_area();
792  empty_area = bge.empty_area();
793  n_empty_jets = bge.n_empty_jets();
794  } else {
795  assert(do_bkgd_gridmedian);
796  double grid_rapmin, grid_rapmax;
797  bkgd_range.get_rapidity_extent(grid_rapmin, grid_rapmax);
798  GridMedianBackgroundEstimator bge(grid_rapmax, 2*ktR);
799  bge.set_particles(particles);
800  rho = bge.rho();
801  sigma = bge.sigma();
802  mean_area = bge.mean_area();
803  empty_area = 0;
804  n_empty_jets = 0;
805  }
806  cout << " rho = " << rho
807  << ", sigma = " << sigma
808  << ", mean_area = " << mean_area
809  << ", empty_area = " << empty_area
810  << ", n_empty_jets = " << n_empty_jets
811  << endl;
812  }
813  } // try
814  catch (fastjet::Error fjerr) {
815  cout << "Caught fastjet error, exiting gracefully" << endl;
816  exit(0);
817  }
818 
819  } // irepeat
820  } // iev
821  // if we've instantiated a plugin, delete it
822  if (jet_def.strategy()==plugin_strategy){
823  delete jet_def.plugin();
824  }
825  // close any file that we've opened
826  if (istr != &cin) delete istr;
827  } // jet_defs
828 
829 }
830 
831 
832 
833 
834 //------ HELPER ROUTINES -----------------------------------------------
835 /// print a single jet
836 void print_jet (const PseudoJet & jet) {
837  unsigned int n_constituents = jet.constituents().size();
838  printf("%15.8f %15.8f %15.8f %8u\n",
839  jet.rap(), jet.phi(), jet.perp(), n_constituents);
840 }
841 
842 
843 //----------------------------------------------------------------------
844 void print_jets(const vector<PseudoJet> & jets_in, bool show_constituents) {
845  vector<PseudoJet> jets;
846  if (ee_print) {
847  jets = sorted_by_E(jets_in);
848  for (unsigned int j = 0; j < jets.size(); j++) {
849  printf("%5u %15.8f %15.8f %15.8f %15.8f\n",
850  j,jets[j].px(),jets[j].py(),jets[j].pz(),jets[j].E());
851  if (show_constituents) {
852  vector<PseudoJet> const_jets = jets[j].constituents();
853  for (unsigned int k = 0; k < const_jets.size(); k++) {
854  printf(" jet%03u %15.8f %15.8f %15.8f %15.8f\n",j,const_jets[k].px(),
855  const_jets[k].py(),const_jets[k].pz(),const_jets[k].E());
856  }
857  cout << "\n\n";
858  }
859 
860  }
861  } else {
862  jets = sorted_by_pt(jets_in);
863  for (unsigned int j = 0; j < jets.size(); j++) {
864  printf("%5u %15.8f %15.8f %15.8f",
865  j,jets[j].rap(),jets[j].phi(),jets[j].perp());
866  // also print out the scalar area and the perp component of the
867  // 4-vector (just enough to check a reasonable 4-vector?)
868  if (do_areas) printf(" %15.8f %15.8f", jets[j].area(),
869  jets[j].area_4vector().perp());
870  cout << "\n";
871 
872  if (show_constituents) {
873  vector<PseudoJet> const_jets = jets[j].constituents();
874  for (unsigned int k = 0; k < const_jets.size(); k++) {
875  printf(" jet%03u %15.8f %15.8f %15.8f %5d\n",j,const_jets[k].rap(),
876  const_jets[k].phi(),sqrt(const_jets[k].kt2()), const_jets[k].cluster_hist_index());
877  }
878  cout << "\n\n";
879  }
880  }
881  }
882 
883  if (rootfile != "") {
884  ofstream ostr(rootfile.c_str());
885  ostr << "# " << cmdline_p->command_line() << endl;
886  ostr << "# output for root" << endl;
887  assert(jets.size() > 0);
888  jets[0].validated_cs()->print_jets_for_root(jets,ostr);
889  }
890 
891 }
892 
893 
894 //----- SUBJETS --------------------------------------------------------
895 /// a function that pretty prints a list of jets and the subjets for each
896 /// one
897 void print_jets_and_sub (const vector<PseudoJet> & jets, double dcut) {
898 
899  // sort jets into increasing pt
900  vector<PseudoJet> sorted_jets = sorted_by_pt(jets);
901 
902  // label the columns
903  printf("Printing jets and their subjets with subdcut = %10.5f\n",dcut);
904  printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity",
905  "phi", "pt", "n constituents");
906 
907  // the kind of subjet finding used to test consistency among them
908  SubType sub_type = subtype_internal;
909  //SubType sub_type = subtype_newclust_dcut;
910  //SubType sub_type = subtype_newclust_R;
911 
912  // print out the details for each jet
913  //for (unsigned int i = 0; i < sorted_jets.size(); i++) {
914  for (vector<PseudoJet>::const_iterator jet = sorted_jets.begin();
915  jet != sorted_jets.end(); jet++) {
916  const JetDefinition & jet_def = jet->validated_cs()->jet_def();
917 
918  // if jet pt^2 < dcut with kt alg, then some methods of
919  // getting subjets will return nothing -- so skip the jet
920  if (jet_def.jet_algorithm() == kt_algorithm
921  && jet->perp2() < dcut) continue;
922 
923 
924  printf("%5u ",(unsigned int) (jet - sorted_jets.begin()));
925  print_jet(*jet);
926  vector<PseudoJet> subjets;
927  ClusterSequence * cspoint;
928  if (sub_type == subtype_internal) {
929  cspoint = 0;
930  subjets = jet->exclusive_subjets(dcut);
931  double ddnp1 = jet->exclusive_subdmerge_max(subjets.size());
932  double ddn = jet->exclusive_subdmerge_max(subjets.size()-1);
933  cout << " for " << ddnp1 << " < d < " << ddn << " one has " << endl;
934  } else if (sub_type == subtype_newclust_dcut) {
935  cspoint = new ClusterSequence(jet->constituents(), jet_def);
936  subjets = cspoint->exclusive_jets(dcut);
937  } else if (sub_type == subtype_newclust_R) {
938  assert(jet_def.jet_algorithm() == cambridge_algorithm);
939  JetDefinition subjd(jet_def.jet_algorithm(), jet_def.R()*sqrt(dcut));
940  cspoint = new ClusterSequence(jet->constituents(), subjd);
941  subjets = cspoint->inclusive_jets();
942  } else {
943  cerr << "unrecognized subtype for subjet finding" << endl;
944  exit(-1);
945  }
946 
947  subjets = sorted_by_pt(subjets);
948  for (unsigned int j = 0; j < subjets.size(); j++) {
949  printf(" -sub-%02u ",j);
950  print_jet(subjets[j]);
951  }
952 
953  if (cspoint != 0) delete cspoint;
954 
955  //ClusterSequence subseq(clust_seq->constituents(sorted_jets[i]),
956  // JetDefinition(cambridge_algorithm, 0.4));
957  //vector<PseudoJet> subjets = sorted_by_pt(subseq.inclusive_jets());
958  //for (unsigned int j = 0; j < subjets.size(); j++) {
959  // printf(" -sub-%02u ",j);
960  // print_jet(subseq, subjets[j]);
961  //}
962  }
963 
964 }
965 
966 
967 //----------------------------------------------------------------------
968 void signal_failed_comparison(int iev,
969  const string & message,
970  const vector<PseudoJet> & particles) {
971  cout << "# failed comparison, reason is " << message << endl;
972  cout << "# iev = " << iev << endl;
973  cout << setprecision(16);
974  for (unsigned i = 0; i < particles.size(); i++) {
975  const PseudoJet & p = particles[i];
976  cout << p.px() << " "
977  << p.py() << " "
978  << p.pz() << " "
979  << p.E () << endl;
980  }
981  cout << "#END" << endl;
982 }
983 
984 //----------------------------------------------------------------------
985 void do_compare_strategy(int iev,
986  const vector<PseudoJet> & particles,
987  const JetDefinition & jet_def,
988  const ClusterSequence & cs,
989  int compare_strategy) {
990 
991  // create a jet def with the reference comparison strategy
992  JetDefinition jet_def_ref(jet_def.jet_algorithm(),
993  jet_def.R(),
994  jet_def.recombination_scheme(),
995  Strategy(compare_strategy));
996  // do the clustering
997  ClusterSequence cs_ref(particles, jet_def_ref);
998 
999  // now compare the outputs. At this stage just based on the clustering
1000  // sequence - get more sophisticated later...
1001  const vector<ClusterSequence::history_element> & history_in = cs.history();
1002  const vector<ClusterSequence::history_element> & history_ref = cs_ref.history();
1003 
1004  if (history_in.size() != history_ref.size()) {
1005  signal_failed_comparison(iev, "history sizes do not match", particles);
1006  }
1007 
1008  // now run over each clustering step to do the comparison
1009  for (unsigned i = cs.n_particles(); i < history_in.size(); i++) {
1010  bool fail_parents = (history_in[i].parent1 != history_ref[i].parent1 ||
1011  history_in[i].parent2 != history_ref[i].parent2);
1012  bool fail_dij = (history_in[i].dij != history_ref[i].dij);
1013  if (fail_parents || fail_dij) {
1014  ostringstream ostr;
1015  ostr << "at step " << i << ", ";
1016  if (fail_parents) ostr << "parents do not match, ";
1017  if (fail_dij) ostr << "dij does not match, ";
1018  ostr << "history in (p1, p2, dij) = "
1019  << history_in[i].parent1 << " " << history_in[i].parent2 << " " << history_in[i].dij;
1020  ostr << ", history ref (p1, p2, dij) = "
1021  << history_ref[i].parent1 << " " << history_ref[i].parent2 << " " << history_ref[i].dij;
1022  signal_failed_comparison(iev, ostr.str(), particles);
1023  break;
1024  }
1025  }
1026 }
Generated on Fri Aug 15 2014 16:12:48 for FastJet by   doxygen 1.8.4