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