FastJet 3.0.6
fastjet_timing_plugins.cc
00001 //STARTHEADER
00002 // $Id: fastjet_timing_plugins.cc 3209 2013-09-15 20:03:30Z salam $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 
00030 //----------------------------------------------------------------------
00031 /// fastjet_timing.cc: Program to help time and test the fastjet package
00032 /// 
00033 /// It reads files containing multiple events in the format 
00034 /// p1x p1y p1z E1
00035 /// p2x p2y p2z E2
00036 /// ...
00037 /// #END
00038 /// 
00039 /// An example input file containing 10 events is included as 
00040 /// data/Pythia-PtMin1000-LHC-10ev.dat
00041 ///
00042 /// Usage:
00043 ///   fastjet_timing [-strategy NUMBER] [-repeat nrepeats] [-massive] \
00044 ///                  [-combine nevents] [-r Rparameter] [-incl ptmin] [...] \
00045 ///                  < data_file
00046 ///
00047 /// where the clustering can be repeated to aid timing and multiple
00048 /// events can be combined to get to larger multiplicities. Some options:
00049 ///
00050 /// Options for reading
00051 /// -------------------
00052 ///
00053 ///   -nev     n    number of events to run
00054 ///
00055 ///   -combine n    for combining multiple events from the data file in order
00056 ///                 to get a single high-multipicity event to run.
00057 ///
00058 ///   -massless     read in only the 3-momenta and deduce energies assuming
00059 ///                 that particles are massless
00060 ///
00061 ///   -dense        adds dense ghost coverage
00062 ///
00063 ///   -repeat n     repeats each event n times
00064 ///
00065 ///   -nhardest n   keep only the n hardest particles in the event
00066 ///
00067 ///   -file name    read from the corresponding file rather than stdin.
00068 ///                 (The file will be reopened for each new jet alg.; in
00069 ///                 constrast, if you use stdin, each new alg will take a
00070 ///                 new event).
00071 /// 
00072 /// Output Options
00073 /// --------------
00074 ///
00075 ///   -incl ptmin   output of all inclusive jets with pt > ptmin is obtained
00076 ///                 with the -incl option.
00077 ///
00078 ///   -repeat-incl ptmin
00079 ///                 same as -incl ptmin but do it for each repetition
00080 ///                 of the clustering
00081 ///
00082 ///   -excld dcut   output of all exclusive jets as obtained in a clustering
00083 ///                 with dcut
00084 ///
00085 ///   -excly ycut   output of all exclusive jets as obtained in a clustering
00086 ///                 with ycut
00087 ///
00088 ///   -excln n      output of clustering to n exclusive jets
00089 ///
00090 ///   -ee-print     print things as px,py,pz,E
00091 ///
00092 ///   -get-all-dij  print out all dij values
00093 ///   -get-all-yij  print out all yij values
00094 ///
00095 ///   -const        show jet constituents (works with excl jets)
00096 ///
00097 ///   -write        for writing out detailed clustering sequence (valuable
00098 ///                 for testing purposes)
00099 ///
00100 ///   -unique_write writes out the sequence of dij's according to the
00101 ///                 "unique_history_order" (useful for verifying consistency
00102 ///                 between different clustering strategies).
00103 ///
00104 ///   -root file    sends output to file that can be read in with the script in
00105 ///                 root/ so as to show a lego-plot of the event
00106 ///
00107 ///   -cones        show extra info about internal steps for SISCone
00108 ///
00109 ///   -area         calculate areas. Additional options include
00110 ///                   -area:active
00111 ///                   -area:passive
00112 ///                   -area:explicit
00113 ///                   -area:voronoi Rfact
00114 ///                   -area:repeat nrepeat
00115 ///                   -ghost-area area
00116 ///                   -ghost-maxrap maxrap
00117 ///                   -area:fj2               place ghosts as in fj2
00118 ///
00119 ///   -bkgd         calculate the background density. Additional options include
00120 ///                   -bkgd:csab       use the old ClusterSequenceAreaBase methods
00121 ///                   -bkgd:jetmedian  use the new JetMedianBackgroundEstimator class
00122 ///                   -bkgd:fj2        force jetmedian to calculate sigma as in fj2
00123 ///                   -bkgd:gridmedian use GridMedianBackgroundEstimator with grid up to ghost_maxrap-ktR and grid spacing of 2ktR
00124 ///
00125 /// Algorithms
00126 /// ----------
00127 ///   -all-algs     runs all algorithms
00128 ///
00129 ///   -kt           switch to the longitudinally invariant kt algorithm
00130 ///                 Note: this is the default one.
00131 ///
00132 ///   -cam          switch to the inclusive Cambridge/Aachen algorithm --
00133 ///                 note that the option -excld dcut provides a clustering
00134 ///                 up to the dcut which is the minimum squared
00135 ///                 distance between any pair of jets.
00136 ///
00137 ///   -antikt       switch to the anti-kt clustering algorithm
00138 ///
00139 ///   -genkt        switch to the genkt algorithm
00140 ///                 you can provide the parameter of the alg as an argument to 
00141 ///                 -genkt (1 by default)
00142 ///                 
00143 ///   -eekt         switch to the e+e- kt algorithm
00144 ///
00145 ///   -eegenkt      switch to the genkt algorithm
00146 ///                 you can provide the parameter of the alg as an argument to 
00147 ///                 -ee_genkt (1 by default)
00148 ///                 
00149 /// plugins (don't delete this line)
00150 ///
00151 ///   -pxcone             switch to the PxCone jet algorithm
00152 /// 
00153 ///   -siscone            switch to the SISCone jet algorithm (seedless cones)
00154 ///   -sisconespheri      switch to the Spherical SISCone jet algorithm (seedless cones)
00155 ///
00156 ///   -midpoint           switch to CDF's midpoint code
00157 ///   -jetclu             switch to CDF's jetclu code
00158 ///
00159 ///   -d0runipre96cone    switch to the D0RunIpre96Cone plugin
00160 ///   -d0runicone         switch to the D0RunICone plugin
00161 ///
00162 ///   -d0runiicone        switch to D0's run II midpoint cone
00163 ///
00164 ///   -trackjet           switch to the TrackJet plugin
00165 ///
00166 ///   -atlascone          switch to the ATLASCone plugin
00167 ///
00168 ///   -eecambridge        switch to the EECambridge plugin
00169 ///
00170 ///   -jade               switch to the Jade plugin
00171 ///
00172 ///   -cmsiterativecone   switch to the CMSIterativeCone plugin
00173 ///
00174 ///   -gridjet     switch to the GridJet plugin
00175 ///
00176 ///  end of plugins (don't delete this line)
00177 ///
00178 ///
00179 /// Options for running algs
00180 /// ------------------------
00181 ///
00182 ///   -r            sets the radius of the jet algorithm (default = 1.0)
00183 ///
00184 ///   -overlap | -f sets the overlap fraction in cone algs with split-merge
00185 ///
00186 ///   -seed         sets the seed threshold
00187 ///
00188 ///   -strategy N   indicate stratgey from the enum fastjet::Strategy (see
00189 ///                 fastjet/JetDefinition.hh).
00190 ///
00191 
00192 
00193 #include "fastjet/ClusterSequenceArea.hh"
00194 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
00195 #include "fastjet/tools/GridMedianBackgroundEstimator.hh"
00196 #include "fastjet/Selector.hh"
00197 #include<iostream>
00198 #include<sstream>
00199 #include<fstream>
00200 #include<valarray>
00201 #include<vector>
00202 #include <cstdlib>
00203 //#include<cstddef> // for size_t
00204 #include "CmdLine.hh"
00205 
00206 // get info on how fastjet was configured
00207 #include "fastjet/config.h"
00208 
00209 // include the installed plugins (don't delete this line)
00210 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
00211 #include "fastjet/SISConePlugin.hh"
00212 #include "fastjet/SISConeSphericalPlugin.hh"
00213 #endif
00214 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
00215 #include "fastjet/CDFMidPointPlugin.hh"
00216 #include "fastjet/CDFJetCluPlugin.hh"
00217 #endif
00218 #ifdef FASTJET_ENABLE_PLUGIN_PXCONE
00219 #include "fastjet/PxConePlugin.hh"
00220 #endif
00221 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
00222 #include "fastjet/D0RunIIConePlugin.hh"
00223 #endif 
00224 #ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
00225 #include "fastjet/TrackJetPlugin.hh"
00226 #endif
00227 #ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
00228 #include "fastjet/ATLASConePlugin.hh"
00229 #endif
00230 #ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
00231 #include "fastjet/EECambridgePlugin.hh"
00232 #endif
00233 #ifdef FASTJET_ENABLE_PLUGIN_JADE
00234 #include "fastjet/JadePlugin.hh"
00235 #endif
00236 #ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
00237 #include "fastjet/CMSIterativeConePlugin.hh"
00238 #endif
00239 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
00240 #include "fastjet/D0RunIpre96ConePlugin.hh"
00241 #include "fastjet/D0RunIConePlugin.hh"
00242 #endif
00243 #ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
00244 #include "fastjet/GridJetPlugin.hh"
00245 #endif
00246 // end of installed plugins inclusion (don't delete this line)
00247 
00248 using namespace std;
00249 
00250 // to avoid excessive typing, define an abbreviation for the 
00251 // fastjet namespace
00252 namespace fj = fastjet;
00253 
00254 inline double pow2(const double x) {return x*x;}
00255 
00256 // pretty print the jets and their subjets
00257 void print_jets_and_sub (const vector<fj::PseudoJet> & jets, double dcut);
00258 
00259 string rootfile;
00260 CmdLine * cmdline_p;
00261 
00262 bool do_areas;
00263 
00264 /// sort and pretty print jets, with exact behaviour depending on 
00265 /// whether ee_print is true or not
00266 bool ee_print = false;
00267 void print_jets(const vector<fj::PseudoJet> & jets, bool show_const = false);
00268 
00269 bool found_unavailable = false;
00270 void is_unavailable(const string & algname) {
00271   cerr << algname << " requested, but not available for this compilation" << endl;
00272   found_unavailable = true;
00273   //exit(0);
00274 }
00275 
00276 
00277 /// a program to test and time a range of algorithms as implemented or
00278 /// wrapped in fastjet
00279 int main (int argc, char ** argv) {
00280 
00281   fj::ClusterSequence::print_banner();
00282 
00283   CmdLine cmdline(argc,argv);
00284   cmdline_p = &cmdline;
00285   // allow the use to specify the fj::Strategy either through the
00286   // -clever or the -strategy options (both will take numerical
00287   // values); the latter will override the former.
00288   fj::Strategy  strategy  = fj::Strategy(cmdline.int_val("-strategy",
00289                                         cmdline.int_val("-clever", fj::Best)));
00290   int  repeat  = cmdline.int_val("-repeat",1);
00291   int  combine = cmdline.int_val("-combine",1);
00292   bool write   = cmdline.present("-write");
00293   bool unique_write = cmdline.present("-unique_write");
00294   bool hydjet  = cmdline.present("-hydjet");
00295   double ktR   = cmdline.double_val("-r",1.0);
00296   ktR   = cmdline.double_val("-R",ktR); // allow -r and -R
00297   double inclkt = cmdline.double_val("-incl",-1.0);
00298   double repeatinclkt = cmdline.double_val("-repeat-incl",-1.0);
00299   int    excln  = cmdline.int_val   ("-excln",-1);
00300   double excld  = cmdline.double_val("-excld",-1.0);
00301   double excly  = cmdline.double_val("-excly",-1.0);
00302   ee_print = cmdline.present("-ee-print");
00303   bool   get_all_dij   = cmdline.present("-get-all-dij");
00304   bool   get_all_yij   = cmdline.present("-get-all-yij");
00305   double subdcut = cmdline.double_val("-subdcut",-1.0);
00306   double rapmax = cmdline.double_val("-rapmax",1.0e305);
00307   if (cmdline.present("-etamax")) {
00308     cerr << "WARNING: -etamax options actually sets maximum rapidity (and overrides -rapmax)\n";
00309     rapmax = cmdline.double_val("-etamax");
00310   }
00311   bool   show_constituents = cmdline.present("-const");
00312   bool   massless = cmdline.present("-massless");
00313   int    nev     = cmdline.int_val("-nev",1);
00314   bool   add_dense_coverage = cmdline.present("-dense");
00315   double ghost_maxrap = cmdline.value("-ghost-maxrap",5.0);
00316   bool   all_algs = cmdline.present("-all-algs");
00317 
00318   fj::Selector particles_sel = (cmdline.present("-nhardest"))
00319     ? fj::SelectorNHardest(cmdline.value<unsigned int>("-nhardest"))
00320     : fj::SelectorIdentity();
00321 
00322   do_areas = cmdline.present("-area");
00323   fj::AreaDefinition area_def;
00324   if (do_areas) {
00325     assert(!write); // it's incompatible
00326     fj::GhostedAreaSpec ghost_spec(ghost_maxrap, 
00327                                    cmdline.value("-area:repeat", 1),
00328                                    cmdline.value("-ghost-area", 0.01));
00329     if (cmdline.present("-area:fj2")) ghost_spec.set_fj2_placement(true);
00330     if (cmdline.present("-area:explicit")) {
00331       area_def = fj::AreaDefinition(fj::active_area_explicit_ghosts, ghost_spec);
00332     } else if (cmdline.present("-area:passive")) {
00333       area_def = fj::AreaDefinition(fj::passive_area, ghost_spec);
00334     } else if (cmdline.present("-area:voronoi")) {
00335       double Rfact = cmdline.value<double>("-area:voronoi");
00336       area_def = fj::AreaDefinition(fj::voronoi_area, 
00337                                     fj::VoronoiAreaSpec(Rfact));
00338     } else {
00339       cmdline.present("-area:active"); // allow, but do not require, arg
00340       area_def = fj::AreaDefinition(fj::active_area, ghost_spec);
00341     }
00342   }
00343   bool do_bkgd = cmdline.present("-bkgd"); // background estimation
00344   bool do_bkgd_csab = false, do_bkgd_jetmedian = false, do_bkgd_fj2 = false;
00345   bool do_bkgd_gridmedian = false;
00346   fj::Selector bkgd_range;
00347   if (do_bkgd) {
00348     bkgd_range = fj::SelectorAbsRapMax(ghost_maxrap - ktR); 
00349     if      (cmdline.present("-bkgd:csab"))      {do_bkgd_csab = true;}
00350     else if (cmdline.present("-bkgd:jetmedian")) {do_bkgd_jetmedian = true;
00351       do_bkgd_fj2 = cmdline.present("-bkgd:fj2");
00352     } else if (cmdline.present("-bkgd:gridmedian")) {do_bkgd_gridmedian = true;
00353     } else {
00354       throw fj::Error("with the -bkgd option, some particular background must be specified (csab or jetmedian)");
00355     }
00356     assert(do_areas || do_bkgd_gridmedian);
00357   }
00358 
00359   bool show_cones = cmdline.present("-cones"); // only works for siscone
00360 
00361   // for cone algorithms
00362   // allow -f and -overlap
00363   double overlap_threshold = cmdline.double_val("-overlap",0.5);
00364   overlap_threshold = cmdline.double_val("-f",overlap_threshold); 
00365   double seed_threshold = cmdline.double_val("-seed",1.0);
00366 
00367   // for ee algorithms, allow to specify ycut
00368   double ycut = cmdline.double_val("-ycut",0.08);
00369 
00370   // for printing jets to a file for reading by root
00371   rootfile = cmdline.value<string>("-root","");
00372 
00373   // out default scheme is the E_scheme
00374   fj::RecombinationScheme scheme = fj::E_scheme;
00375 
00376   // The following option causes the Cambridge algo to be used.
00377   // Note that currently the only output that works sensibly here is
00378   // "-incl 0"
00379   vector<fj::JetDefinition> jet_defs;
00380   if (all_algs || cmdline.present("-cam") || cmdline.present("-CA")) {
00381     jet_defs.push_back( fj::JetDefinition(fj::cambridge_algorithm, ktR, scheme, strategy));
00382   } 
00383   if (all_algs || cmdline.present("-antikt")) {
00384     jet_defs.push_back( fj::JetDefinition(fj::antikt_algorithm, ktR, scheme, strategy));
00385   } 
00386   if (all_algs || cmdline.present("-genkt")) {
00387     double p;
00388     if (cmdline.present("-genkt")) p = cmdline.value<double>("-genkt");
00389     else                           p = -0.5;
00390     jet_defs.push_back( fj::JetDefinition(fj::genkt_algorithm, ktR, p, scheme, strategy));
00391   } 
00392   if (all_algs || cmdline.present("-eekt")) {
00393     jet_defs.push_back( fj::JetDefinition(fj::ee_kt_algorithm));
00394   } 
00395   if (all_algs || cmdline.present("-eegenkt")) {
00396     double p;
00397     if (cmdline.present("-eegenkt")) p = cmdline.value<double>("-eegenkt");
00398     else                             p = -0.5;
00399     jet_defs.push_back( fj::JetDefinition(fj::ee_genkt_algorithm, ktR, p, scheme, strategy));
00400 
00401 // checking if one asks to run a plugin (don't delete this line)
00402   } 
00403   if (all_algs || cmdline.present("-midpoint")) {
00404 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
00405     typedef fj::CDFMidPointPlugin MPPlug; // for brevity
00406     double cone_area_fraction = 1.0;
00407     int    max_pair_size = 2;
00408     int    max_iterations = 100;
00409     MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
00410     if (cmdline.present("-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
00411     if (cmdline.present("-sm-pt")) sm_scale = MPPlug::SM_pt; // default
00412     if (cmdline.present("-sm-mt")) sm_scale = MPPlug::SM_mt;
00413     if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et;
00414     jet_defs.push_back( fj::JetDefinition( new fj::CDFMidPointPlugin (
00415                                       seed_threshold, ktR, 
00416                                       cone_area_fraction, max_pair_size,
00417                                       max_iterations, overlap_threshold,
00418                                       sm_scale)));
00419 #else  // FASTJET_ENABLE_PLUGIN_CDFCONES
00420     is_unavailable("MidPoint");
00421 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
00422   } 
00423   if (all_algs || cmdline.present("-pxcone")) {
00424 #ifdef FASTJET_ENABLE_PLUGIN_PXCONE
00425     double min_jet_energy = 5.0;
00426     jet_defs.push_back( fj::JetDefinition( new fj::PxConePlugin (
00427                                       ktR, min_jet_energy,
00428                                       overlap_threshold)));
00429 #else  // FASTJET_ENABLE_PLUGIN_PXCONE
00430     is_unavailable("PxCone");
00431 #endif // FASTJET_ENABLE_PLUGIN_PXCONE
00432   } 
00433   if (all_algs || cmdline.present("-jetclu")) {
00434 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
00435     jet_defs.push_back( fj::JetDefinition( new fj::CDFJetCluPlugin (
00436                                                                     ktR, overlap_threshold, seed_threshold)));
00437 #else  // FASTJET_ENABLE_PLUGIN_CDFCONES
00438     is_unavailable("JetClu");
00439 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
00440   } 
00441   if (all_algs || cmdline.present("-siscone") || cmdline.present("-sisconespheri")) {
00442 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
00443     typedef fj::SISConePlugin SISPlug; // for brevity
00444     int npass = cmdline.value("-npass",0);
00445     if (all_algs || cmdline.present("-siscone")) {
00446       double sisptmin = cmdline.value("-sisptmin",0.0);
00447       SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass,sisptmin);
00448       if (cmdline.present("-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
00449       if (cmdline.present("-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
00450       if (cmdline.present("-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
00451       if (cmdline.present("-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
00452       // cause it to use the jet-definition's own recombiner
00453       plugin->set_use_jet_def_recombiner(true);
00454       jet_defs.push_back( fj::JetDefinition(plugin));
00455     } 
00456     if (all_algs || cmdline.present("-sisconespheri")) {
00457       double sisEmin = cmdline.value("-sisEmin",0.0);
00458       fj::SISConeSphericalPlugin * plugin = 
00459         new fj::SISConeSphericalPlugin(ktR, overlap_threshold,npass,sisEmin);
00460       if (cmdline.present("-ghost-sep")) {
00461         plugin->set_ghost_separation_scale(cmdline.value<double>("-ghost-sep"));
00462       }
00463       jet_defs.push_back( fj::JetDefinition(plugin));
00464     }
00465 #else  // FASTJET_ENABLE_PLUGIN_SISCONE
00466     is_unavailable("SISCone");
00467 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
00468   } 
00469   if (all_algs || cmdline.present("-d0runiicone")) {
00470 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
00471     double min_jet_Et = 6.0; // was 8 GeV in earlier work
00472     jet_defs.push_back( fj::JetDefinition(new fj::D0RunIIConePlugin(ktR,min_jet_Et)));
00473 #else  // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
00474     is_unavailable("D0RunIICone");
00475 #endif // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
00476   } 
00477   if (all_algs || cmdline.present("-trackjet")) {
00478 #ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
00479     jet_defs.push_back( fj::JetDefinition(new fj::TrackJetPlugin(ktR)));
00480 #else  // FASTJET_ENABLE_PLUGIN_TRACKJET
00481     is_unavailable("TrackJet");
00482 #endif // FASTJET_ENABLE_PLUGIN_TRACKJET
00483   } 
00484   if (all_algs || cmdline.present("-atlascone")) {
00485 #ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
00486     jet_defs.push_back( fj::JetDefinition(new fj::ATLASConePlugin(ktR)));
00487 #else  // FASTJET_ENABLE_PLUGIN_ATLASCONE
00488     is_unavailable("ATLASCone");
00489 #endif // FASTJET_ENABLE_PLUGIN_ATLASCONE
00490   } 
00491   if (all_algs || cmdline.present("-eecambridge")) {
00492 #ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
00493     jet_defs.push_back( fj::JetDefinition(new fj::EECambridgePlugin(ycut)));
00494 #else  // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
00495     is_unavailable("EECambridge");
00496 #endif // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
00497   } 
00498   if (all_algs || cmdline.present("-jade")) {
00499 #ifdef FASTJET_ENABLE_PLUGIN_JADE
00500     jet_defs.push_back( fj::JetDefinition(new fj::JadePlugin()));
00501 #else  // FASTJET_ENABLE_PLUGIN_JADE
00502     is_unavailable("Jade");
00503 #endif // FASTJET_ENABLE_PLUGIN_JADE
00504   } 
00505   if (all_algs || cmdline.present("-cmsiterativecone")) {
00506 #ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
00507     jet_defs.push_back( fj::JetDefinition(new fj::CMSIterativeConePlugin(ktR,seed_threshold)));
00508 #else  // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
00509     is_unavailable("CMSIterativeCone");
00510 #endif // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
00511   } 
00512   if (all_algs || cmdline.present("-d0runipre96cone")) {
00513 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
00514     jet_defs.push_back( fj::JetDefinition(new fj::D0RunIpre96ConePlugin(ktR, seed_threshold, overlap_threshold)));
00515 #else  // FASTJET_ENABLE_PLUGIN_D0RUNICONE
00516     is_unavailable("D0RunIpre96Cone");
00517 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
00518   } 
00519   if (all_algs || cmdline.present("-d0runicone")) {
00520 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
00521     jet_defs.push_back( fj::JetDefinition(new fj::D0RunIConePlugin(ktR, seed_threshold, overlap_threshold)));
00522 #else  // FASTJET_ENABLE_PLUGIN_D0RUNICONE
00523     is_unavailable("D0RunICone");
00524 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
00525   } 
00526   if (all_algs || cmdline.present("-gridjet")) {
00527 #ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
00528     // we want a grid_ymax of 5.0, but when using R=0.4 (i.e. grid
00529     // spacing of 0.8), this leads to 12.5 grid cells; depending on
00530     // whether this is 12.499999999999 or 12.5000000....1 this gets
00531     // converted either to 12 or 13, making the results sensitive to
00532     // rounding errors.
00533     //
00534     // Instead we therefore take 4.9999999999, which avoids this problem.
00535     double grid_ymax = 4.9999999999;
00536     jet_defs.push_back( fj::JetDefinition(new fj::GridJetPlugin(grid_ymax, ktR*2.0)));
00537 #else  // FASTJET_ENABLE_PLUGIN_GRIDJET
00538     is_unavailable("GridJet");
00539 #endif // FASTJET_ENABLE_PLUGIN_GRIDJET
00540 // end of checking if one asks to run a plugin (don't delete this line)
00541   } 
00542   if (all_algs || 
00543       cmdline.present("-kt") || 
00544       (jet_defs.size() == 0 && !found_unavailable))  {
00545     jet_defs.push_back( fj::JetDefinition(fj::kt_algorithm, ktR, strategy));
00546   }
00547 
00548   string filename = cmdline.value<string>("-file", "");
00549 
00550 
00551   if (!cmdline.all_options_used()) {cerr << 
00552       "Error: some options were not recognized"<<endl; 
00553     exit(-1);}
00554 
00555   for (unsigned idef = 0; idef < jet_defs.size(); idef++) {
00556   fj::JetDefinition & jet_def = jet_defs[idef];
00557   istream * istr;
00558   if (filename == "") istr = &cin;
00559   else                istr = new ifstream(filename.c_str());
00560 
00561   for (int iev = 0; iev < nev; iev++) {
00562   vector<fj::PseudoJet> jets;
00563   vector<fj::PseudoJet> particles;
00564   string line;
00565   int  ndone = 0;
00566   while (getline(*istr, line)) {
00567       //cout << line<<endl;
00568     istringstream linestream(line);
00569     if (line == "#END") {
00570       ndone += 1;
00571       if (ndone == combine) {break;}
00572     }
00573     if (line.substr(0,1) == "#") {continue;}
00574     valarray<double> fourvec(4);
00575     if (hydjet) {
00576       // special reading from hydjet.txt event record (though actually
00577       // this is supposed to be a standard pythia event record, so
00578       // being able to read from it is perhaps not so bad an idea...)
00579       int ii, istat,id,m1,m2,d1,d2;
00580       double mass;
00581       linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00582                  >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00583       // current file contains mass of particle as 4th entry
00584       if (istat == 1) {
00585         fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00586                           +pow2(fourvec[2])+pow2(mass));
00587       }
00588     } else {
00589       if (massless) {
00590         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00591         fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00592       else {
00593         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00594       }
00595     }
00596     fj::PseudoJet psjet(fourvec);
00597     if (abs(psjet.rap()) < rapmax) {particles.push_back(psjet);}
00598   }
00599 
00600   // add a fake underlying event which is very soft, uniformly distributed
00601   // in eta,phi so as to allow one to reconstruct the area that is associated
00602   // with each jet.
00603   if (add_dense_coverage) {
00604     fj::GhostedAreaSpec ghosted_area_spec(ghost_maxrap);
00605     //fj::GhostedAreaSpec ghosted_area_spec(-2.0,4.0); // asymmetric range
00606     // for plots, reduce the scatter default of 1, to avoid "holes"
00607     // in the subsequent calorimeter view
00608     ghosted_area_spec.set_grid_scatter(0.5); 
00609     ghosted_area_spec.add_ghosts(particles);
00610     //----- old code ------------------
00611     // srand(2);
00612     // int nphi = 60;
00613     // int neta = 100;
00614     // double kt = 1e-1;
00615     // for (int iphi = 0; iphi<nphi; iphi++) {
00616     //   for (int ieta = -neta; ieta<neta+1; ieta++) {
00617     //  double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00618     //  double eta = ieta * (10.0/neta)  + rand()*0.001/RAND_MAX;
00619     //  kt = 1e-20*(1+rand()*0.1/RAND_MAX);
00620     //  double pminus = kt*exp(-eta);
00621     //  double pplus  = kt*exp(+eta);
00622     //  double px = kt*sin(phi);
00623     //  double py = kt*cos(phi);
00624     //  //cout << kt<<" "<<eta<<" "<<phi<<"\n";
00625     //  fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00626     //  particles.push_back(mom);
00627     //   }
00628     // }
00629   }
00630 
00631   // select the particles that pass the selection cut
00632   particles = particles_sel(particles);
00633   
00634   for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00635     int nparticles = particles.size();
00636     try {
00637     auto_ptr<fj::ClusterSequence> clust_seq;
00638     if (do_areas) {
00639       clust_seq.reset(new fj::ClusterSequenceArea(particles,jet_def,area_def));
00640     } else {
00641       clust_seq.reset(new fj::ClusterSequence(particles,jet_def,write));
00642     }
00643 
00644     // repetitive output
00645     if (repeatinclkt >= 0.0) {
00646       vector<fj::PseudoJet> jets_local = sorted_by_pt(clust_seq->inclusive_jets(repeatinclkt));
00647     }
00648 
00649     if (irepeat != 0) {continue;}
00650     cout << "iev "<<iev<< ": number of particles = "<< nparticles << endl;
00651     cout << "strategy used =  "<< clust_seq->strategy_string()<< endl;
00652     if (iev == 0) cout << "Jet Definition: " << jet_def.description() << " (" << fj::fastjet_version_string() << ")" << endl;
00653     if (do_areas && iev == 0) cout << "Area definition: " << area_def.description() << endl;
00654 
00655     // now provide some nice output...
00656     if (inclkt >= 0.0) {
00657       vector<fj::PseudoJet> jets_local = sorted_by_pt(clust_seq->inclusive_jets(inclkt));
00658       print_jets(jets_local, show_constituents);
00659 
00660     }
00661 
00662     if (excln > 0) {
00663       cout << "Printing "<<excln<<" exclusive jets\n";
00664       print_jets(clust_seq->exclusive_jets(excln), show_constituents);
00665     }
00666 
00667     if (excld > 0.0) {
00668       cout << "Printing exclusive jets for d = "<<excld<<"\n";
00669       print_jets(clust_seq->exclusive_jets(excld), show_constituents);
00670     }
00671 
00672     if (excly > 0.0) {
00673       cout << "Printing exclusive jets for ycut = "<<excly<<"\n";
00674       print_jets(clust_seq->exclusive_jets_ycut(excly), show_constituents);
00675     }
00676 
00677     if (get_all_dij) {
00678       for (int i = nparticles-1; i >= 0; i--) {
00679         printf("d for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_dmerge(i));
00680       }
00681     }
00682     if (get_all_yij) {
00683       for (int i = nparticles-1; i >= 0; i--) {
00684         printf("y for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_ymerge(i));
00685       }
00686     }
00687 
00688     // have the option of printing out the subjets (at scale dcut) of
00689     // each inclusive jet
00690     if (subdcut >= 0.0) {
00691       print_jets_and_sub(clust_seq->inclusive_jets(), subdcut);
00692     }
00693     
00694     // useful for testing that recombination sequences are unique
00695     if (unique_write) {
00696       vector<int> unique_history = clust_seq->unique_history_order();
00697       // construct the inverse of the above mapping
00698       vector<int> inv_unique_history(clust_seq->history().size());
00699       for (unsigned int i = 0; i < unique_history.size(); i++) {
00700         inv_unique_history[unique_history[i]] = i;}
00701 
00702       for (unsigned int i = 0; i < unique_history.size(); i++) {
00703         fj::ClusterSequence::history_element el = 
00704           clust_seq->history()[unique_history[i]];
00705         int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00706         int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00707         printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00708       }
00709     }
00710 
00711 
00712 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
00713     // provide some complementary information for SISCone 
00714     if (show_cones) {
00715       const fj::SISConeExtras * extras = 
00716         dynamic_cast<const fj::SISConeExtras *>(clust_seq->extras());
00717       if (extras == 0) 
00718         throw fastjet::Error("extras object for SISCone was null (this can happen with certain area types)");
00719       cout << "most ambiguous split (difference in squared dist) = "
00720            << extras->most_ambiguous_split() << endl;
00721       vector<fastjet::PseudoJet> stable_cones(extras->stable_cones()); 
00722       stable_cones = sorted_by_rapidity(stable_cones);
00723       for (unsigned int i = 0; i < stable_cones.size(); i++) {
00724       //if (stable_cones[i].phi() < 5.0 && stable_cones[i].phi() > 4.0) {
00725         printf("%5u %15.8f %15.8f %15.8f\n",
00726                i,stable_cones[i].rap(),stable_cones[i].phi(),
00727                stable_cones[i].perp() );
00728       //}
00729       }
00730       
00731       // also show passes for jets
00732       vector<fj::PseudoJet> sisjets = clust_seq->inclusive_jets();
00733       printf("\n%15s %15s %15s %12s %8s %8s\n","rap","phi","pt","user-index","pass","nconst");
00734       for (unsigned i = 0; i < sisjets.size(); i++) {
00735         printf("%15.8f %15.8f %15.8f %12d %8d %8u\n",
00736                sisjets[i].rap(), sisjets[i].phi(), sisjets[i].perp(), 
00737                sisjets[i].user_index(), extras->pass(sisjets[i]),
00738                (unsigned int) clust_seq->constituents(sisjets[i]).size()
00739                );
00740         
00741       }
00742     }
00743 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
00744 
00745     if (do_bkgd) {
00746       double rho, sigma, mean_area, empty_area, n_empty_jets;
00747       fj::ClusterSequenceAreaBase * csab = 
00748         dynamic_cast<fj::ClusterSequenceAreaBase *>(clust_seq.get());
00749       if (do_bkgd_csab) {
00750         csab->get_median_rho_and_sigma(bkgd_range, true, rho, sigma, mean_area);
00751         empty_area = csab->empty_area(bkgd_range);
00752         n_empty_jets = csab->n_empty_jets(bkgd_range);
00753       } else if (do_bkgd_jetmedian) {
00754         fj::JetMedianBackgroundEstimator bge(bkgd_range);
00755         bge.set_provide_fj2_sigma(do_bkgd_fj2);
00756         bge.set_cluster_sequence(*csab);
00757         rho = bge.rho();
00758         sigma = bge.sigma();
00759         mean_area = bge.mean_area();
00760         empty_area = bge.empty_area();
00761         n_empty_jets = bge.n_empty_jets();
00762       } else {
00763         assert(do_bkgd_gridmedian);
00764         double rapmin, rapmax;
00765         bkgd_range.get_rapidity_extent(rapmin, rapmax);
00766         fj::GridMedianBackgroundEstimator bge(rapmax, 2*ktR);
00767         bge.set_particles(particles);
00768         rho = bge.rho();
00769         sigma = bge.sigma();
00770         mean_area = bge.mean_area();
00771         empty_area = 0;
00772         n_empty_jets = 0;
00773       }
00774       cout << "  rho = " << rho 
00775            << ", sigma = " << sigma 
00776            << ", mean_area = " << mean_area
00777            << ", empty_area = " << empty_area
00778            << ", n_empty_jets = " << n_empty_jets
00779            << endl;
00780     }
00781   } // try
00782   catch (fastjet::Error fjerr) {
00783     cout << "Caught fastjet error, exiting gracefully" << endl;
00784     exit(0);
00785   }
00786 
00787   } // irepeat
00788   } // iev
00789   // if we've instantiated a plugin, delete it
00790   if (jet_def.strategy()==fj::plugin_strategy){
00791     delete jet_def.plugin();
00792   }
00793   // close any file that we've opened
00794   if (istr != &cin) delete istr;
00795   } // jet_defs
00796 
00797 }
00798 
00799 
00800 
00801 
00802 //------ HELPER ROUTINES -----------------------------------------------
00803 /// print a single jet
00804 void print_jet (const fj::PseudoJet & jet) {
00805   unsigned int n_constituents = jet.constituents().size();
00806   printf("%15.8f %15.8f %15.8f %8u\n",
00807          jet.rap(), jet.phi(), jet.perp(), n_constituents);
00808 }
00809 
00810 
00811 //----------------------------------------------------------------------
00812 void print_jets(const vector<fj::PseudoJet> & jets_in, bool show_constituents) {
00813   vector<fj::PseudoJet> jets;
00814   if (ee_print) {
00815     jets = sorted_by_E(jets_in);
00816     for (unsigned int j = 0; j < jets.size(); j++) {
00817       printf("%5u %15.8f %15.8f %15.8f %15.8f\n",
00818              j,jets[j].px(),jets[j].py(),jets[j].pz(),jets[j].E());
00819       if (show_constituents) {
00820         vector<fj::PseudoJet> const_jets = jets[j].constituents();
00821         for (unsigned int k = 0; k < const_jets.size(); k++) {
00822           printf("        jet%03u %15.8f %15.8f %15.8f %15.8f\n",j,const_jets[k].px(),
00823                  const_jets[k].py(),const_jets[k].pz(),const_jets[k].E());
00824         }
00825         cout << "\n\n";
00826     }
00827 
00828     }
00829   } else {
00830     jets = sorted_by_pt(jets_in);
00831     for (unsigned int j = 0; j < jets.size(); j++) {
00832       printf("%5u %15.8f %15.8f %15.8f",
00833              j,jets[j].rap(),jets[j].phi(),jets[j].perp());
00834       // also print out the scalar area and the perp component of the
00835       // 4-vector (just enough to check a reasonable 4-vector?)
00836       if (do_areas) printf(" %15.8f %15.8f", jets[j].area(),
00837                                              jets[j].area_4vector().perp());
00838       cout << "\n";
00839 
00840       if (show_constituents) {
00841         vector<fj::PseudoJet> const_jets = jets[j].constituents();
00842         for (unsigned int k = 0; k < const_jets.size(); k++) {
00843           printf("        jet%03u %15.8f %15.8f %15.8f %5d\n",j,const_jets[k].rap(),
00844                  const_jets[k].phi(),sqrt(const_jets[k].kt2()), const_jets[k].cluster_hist_index());
00845         }
00846         cout << "\n\n";
00847       }
00848     }
00849   }
00850 
00851   if (rootfile != "") {
00852     ofstream ostr(rootfile.c_str());
00853     ostr << "# " << cmdline_p->command_line() << endl;
00854     ostr << "# output for root" << endl;
00855     assert(jets.size() > 0);
00856     jets[0].validated_cs()->print_jets_for_root(jets,ostr);
00857   }
00858 
00859 }
00860 
00861 
00862 //----- SUBJETS --------------------------------------------------------
00863 /// a function that pretty prints a list of jets and the subjets for each
00864 /// one
00865 void print_jets_and_sub (const vector<fj::PseudoJet> & jets, double dcut) {
00866 
00867   // sort jets into increasing pt
00868   vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);  
00869 
00870   // label the columns
00871   printf("Printing jets and their subjets with subdcut = %10.5f\n",dcut);
00872   printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", 
00873          "phi", "pt", "n constituents");
00874 
00875   // have various kinds of subjet finding, to test consistency among them
00876   enum SubType {internal, newclust_dcut, newclust_R};
00877   SubType subtype = internal;
00878   //SubType subtype = newclust_dcut;
00879   //SubType subtype = newclust_R;
00880 
00881   // print out the details for each jet
00882   //for (unsigned int i = 0; i < sorted_jets.size(); i++) {
00883   for (vector<fj::PseudoJet>::const_iterator jet = sorted_jets.begin(); 
00884        jet != sorted_jets.end(); jet++) {
00885     const fj::JetDefinition & jet_def = jet->validated_cs()->jet_def();
00886 
00887     // if jet pt^2 < dcut with kt alg, then some methods of
00888     // getting subjets will return nothing -- so skip the jet
00889     if (jet_def.jet_algorithm() == fj::kt_algorithm 
00890         && jet->perp2() < dcut) continue;
00891 
00892 
00893     printf("%5u       ",(unsigned int) (jet - sorted_jets.begin()));
00894     print_jet(*jet);
00895     vector<fj::PseudoJet> subjets;
00896     fj::ClusterSequence * cspoint;
00897     if (subtype == internal) {
00898       cspoint = 0;
00899       subjets = jet->exclusive_subjets(dcut);
00900       double ddnp1 = jet->exclusive_subdmerge_max(subjets.size());
00901       double ddn   = jet->exclusive_subdmerge_max(subjets.size()-1);
00902       cout << "     for " << ddnp1 << " < d < " << ddn << " one has " << endl;
00903     } else if (subtype == newclust_dcut) {
00904       cspoint = new fj::ClusterSequence(jet->constituents(), jet_def);
00905       subjets = cspoint->exclusive_jets(dcut);
00906     } else if (subtype == newclust_R) {
00907       assert(jet_def.jet_algorithm() == fj::cambridge_algorithm);
00908       fj::JetDefinition subjd(jet_def.jet_algorithm(), jet_def.R()*sqrt(dcut));
00909       cspoint = new fj::ClusterSequence(jet->constituents(), subjd);
00910       subjets = cspoint->inclusive_jets();
00911     } else {
00912       cerr << "unrecognized subtype for subjet finding" << endl;
00913       exit(-1);
00914     }
00915 
00916     subjets = sorted_by_pt(subjets);
00917     for (unsigned int j = 0; j < subjets.size(); j++) {
00918       printf("    -sub-%02u ",j);
00919       print_jet(subjets[j]);
00920     }
00921 
00922     if (cspoint != 0) delete cspoint;
00923 
00924     //fj::ClusterSequence subseq(clust_seq->constituents(sorted_jets[i]),
00925     //                          fj::JetDefinition(fj::cambridge_algorithm, 0.4));
00926     //vector<fj::PseudoJet> subjets = sorted_by_pt(subseq.inclusive_jets());
00927     //for (unsigned int j = 0; j < subjets.size(); j++) {
00928     //  printf("    -sub-%02u ",j);
00929     //  print_jet(subseq, subjets[j]);
00930     //}
00931   }
00932 
00933 }
00934 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends