FastJet 3.0.4
fastjet_timing_plugins.cc
00001 //STARTHEADER
00002 // $Id: fastjet_timing_plugins.cc 3132 2013-06-05 09:27:27Z 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 etamax = cmdline.double_val("-etamax",1.0e305);
00307   bool   show_constituents = cmdline.present("-const");
00308   bool   massless = cmdline.present("-massless");
00309   int    nev     = cmdline.int_val("-nev",1);
00310   bool   add_dense_coverage = cmdline.present("-dense");
00311   double ghost_maxrap = cmdline.value("-ghost-maxrap",5.0);
00312   bool   all_algs = cmdline.present("-all-algs");
00313 
00314   fj::Selector particles_sel = (cmdline.present("-nhardest"))
00315     ? fj::SelectorNHardest(cmdline.value<unsigned int>("-nhardest"))
00316     : fj::SelectorIdentity();
00317 
00318   do_areas = cmdline.present("-area");
00319   fj::AreaDefinition area_def;
00320   if (do_areas) {
00321     assert(!write); // it's incompatible
00322     fj::GhostedAreaSpec ghost_spec(ghost_maxrap, 
00323                                    cmdline.value("-area:repeat", 1),
00324                                    cmdline.value("-ghost-area", 0.01));
00325     if (cmdline.present("-area:fj2")) ghost_spec.set_fj2_placement(true);
00326     if (cmdline.present("-area:explicit")) {
00327       area_def = fj::AreaDefinition(fj::active_area_explicit_ghosts, ghost_spec);
00328     } else if (cmdline.present("-area:passive")) {
00329       area_def = fj::AreaDefinition(fj::passive_area, ghost_spec);
00330     } else if (cmdline.present("-area:voronoi")) {
00331       double Rfact = cmdline.value<double>("-area:voronoi");
00332       area_def = fj::AreaDefinition(fj::voronoi_area, 
00333                                     fj::VoronoiAreaSpec(Rfact));
00334     } else {
00335       cmdline.present("-area:active"); // allow, but do not require, arg
00336       area_def = fj::AreaDefinition(fj::active_area, ghost_spec);
00337     }
00338   }
00339   bool do_bkgd = cmdline.present("-bkgd"); // background estimation
00340   bool do_bkgd_csab = false, do_bkgd_jetmedian = false, do_bkgd_fj2 = false;
00341   bool do_bkgd_gridmedian = false;
00342   fj::Selector bkgd_range;
00343   if (do_bkgd) {
00344     bkgd_range = fj::SelectorAbsRapMax(ghost_maxrap - ktR); 
00345     if      (cmdline.present("-bkgd:csab"))      {do_bkgd_csab = true;}
00346     else if (cmdline.present("-bkgd:jetmedian")) {do_bkgd_jetmedian = true;
00347       do_bkgd_fj2 = cmdline.present("-bkgd:fj2");
00348     } else if (cmdline.present("-bkgd:gridmedian")) {do_bkgd_gridmedian = true;
00349     } else {
00350       throw fj::Error("with the -bkgd option, some particular background must be specified (csab or jetmedian)");
00351     }
00352     assert(do_areas || do_bkgd_gridmedian);
00353   }
00354 
00355   bool show_cones = cmdline.present("-cones"); // only works for siscone
00356 
00357   // for cone algorithms
00358   // allow -f and -overlap
00359   double overlap_threshold = cmdline.double_val("-overlap",0.5);
00360   overlap_threshold = cmdline.double_val("-f",overlap_threshold); 
00361   double seed_threshold = cmdline.double_val("-seed",1.0);
00362 
00363   // for ee algorithms, allow to specify ycut
00364   double ycut = cmdline.double_val("-ycut",0.08);
00365 
00366   // for printing jets to a file for reading by root
00367   rootfile = cmdline.value<string>("-root","");
00368 
00369   // out default scheme is the E_scheme
00370   fj::RecombinationScheme scheme = fj::E_scheme;
00371 
00372   // The following option causes the Cambridge algo to be used.
00373   // Note that currently the only output that works sensibly here is
00374   // "-incl 0"
00375   vector<fj::JetDefinition> jet_defs;
00376   if (all_algs || cmdline.present("-cam") || cmdline.present("-CA")) {
00377     jet_defs.push_back( fj::JetDefinition(fj::cambridge_algorithm, ktR, scheme, strategy));
00378   } 
00379   if (all_algs || cmdline.present("-antikt")) {
00380     jet_defs.push_back( fj::JetDefinition(fj::antikt_algorithm, ktR, scheme, strategy));
00381   } 
00382   if (all_algs || cmdline.present("-genkt")) {
00383     double p;
00384     if (cmdline.present("-genkt")) p = cmdline.value<double>("-genkt");
00385     else                           p = -0.5;
00386     jet_defs.push_back( fj::JetDefinition(fj::genkt_algorithm, ktR, p, scheme, strategy));
00387   } 
00388   if (all_algs || cmdline.present("-eekt")) {
00389     jet_defs.push_back( fj::JetDefinition(fj::ee_kt_algorithm));
00390   } 
00391   if (all_algs || cmdline.present("-eegenkt")) {
00392     double p;
00393     if (cmdline.present("-eegenkt")) p = cmdline.value<double>("-eegenkt");
00394     else                             p = -0.5;
00395     jet_defs.push_back( fj::JetDefinition(fj::ee_genkt_algorithm, ktR, p, scheme, strategy));
00396 
00397 // checking if one asks to run a plugin (don't delete this line)
00398   } 
00399   if (all_algs || cmdline.present("-midpoint")) {
00400 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
00401     typedef fj::CDFMidPointPlugin MPPlug; // for brevity
00402     double cone_area_fraction = 1.0;
00403     int    max_pair_size = 2;
00404     int    max_iterations = 100;
00405     MPPlug::SplitMergeScale sm_scale = MPPlug::SM_pt;
00406     if (cmdline.present("-sm-pttilde")) sm_scale = MPPlug::SM_pttilde;
00407     if (cmdline.present("-sm-pt")) sm_scale = MPPlug::SM_pt; // default
00408     if (cmdline.present("-sm-mt")) sm_scale = MPPlug::SM_mt;
00409     if (cmdline.present("-sm-Et")) sm_scale = MPPlug::SM_Et;
00410     jet_defs.push_back( fj::JetDefinition( new fj::CDFMidPointPlugin (
00411                                       seed_threshold, ktR, 
00412                                       cone_area_fraction, max_pair_size,
00413                                       max_iterations, overlap_threshold,
00414                                       sm_scale)));
00415 #else  // FASTJET_ENABLE_PLUGIN_CDFCONES
00416     is_unavailable("MidPoint");
00417 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
00418   } 
00419   if (all_algs || cmdline.present("-pxcone")) {
00420 #ifdef FASTJET_ENABLE_PLUGIN_PXCONE
00421     double min_jet_energy = 5.0;
00422     jet_defs.push_back( fj::JetDefinition( new fj::PxConePlugin (
00423                                       ktR, min_jet_energy,
00424                                       overlap_threshold)));
00425 #else  // FASTJET_ENABLE_PLUGIN_PXCONE
00426     is_unavailable("PxCone");
00427 #endif // FASTJET_ENABLE_PLUGIN_PXCONE
00428   } 
00429   if (all_algs || cmdline.present("-jetclu")) {
00430 #ifdef FASTJET_ENABLE_PLUGIN_CDFCONES
00431     jet_defs.push_back( fj::JetDefinition( new fj::CDFJetCluPlugin (
00432                                                                     ktR, overlap_threshold, seed_threshold)));
00433 #else  // FASTJET_ENABLE_PLUGIN_CDFCONES
00434     is_unavailable("JetClu");
00435 #endif // FASTJET_ENABLE_PLUGIN_CDFCONES
00436   } 
00437   if (all_algs || cmdline.present("-siscone") || cmdline.present("-sisconespheri")) {
00438 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
00439     typedef fj::SISConePlugin SISPlug; // for brevity
00440     int npass = cmdline.value("-npass",0);
00441     if (all_algs || cmdline.present("-siscone")) {
00442       double sisptmin = cmdline.value("-sisptmin",0.0);
00443       SISPlug * plugin = new SISPlug (ktR, overlap_threshold,npass,sisptmin);
00444       if (cmdline.present("-sm-pt")) plugin->set_split_merge_scale(SISPlug::SM_pt);
00445       if (cmdline.present("-sm-mt")) plugin->set_split_merge_scale(SISPlug::SM_mt);
00446       if (cmdline.present("-sm-Et")) plugin->set_split_merge_scale(SISPlug::SM_Et);
00447       if (cmdline.present("-sm-pttilde")) plugin->set_split_merge_scale(SISPlug::SM_pttilde);
00448       // cause it to use the jet-definition's own recombiner
00449       plugin->set_use_jet_def_recombiner(true);
00450       jet_defs.push_back( fj::JetDefinition(plugin));
00451     } 
00452     if (all_algs || cmdline.present("-sisconespheri")) {
00453       double sisEmin = cmdline.value("-sisEmin",0.0);
00454       fj::SISConeSphericalPlugin * plugin = 
00455         new fj::SISConeSphericalPlugin(ktR, overlap_threshold,npass,sisEmin);
00456       if (cmdline.present("-ghost-sep")) {
00457         plugin->set_ghost_separation_scale(cmdline.value<double>("-ghost-sep"));
00458       }
00459       jet_defs.push_back( fj::JetDefinition(plugin));
00460     }
00461 #else  // FASTJET_ENABLE_PLUGIN_SISCONE
00462     is_unavailable("SISCone");
00463 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
00464   } 
00465   if (all_algs || cmdline.present("-d0runiicone")) {
00466 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNIICONE
00467     double min_jet_Et = 6.0; // was 8 GeV in earlier work
00468     jet_defs.push_back( fj::JetDefinition(new fj::D0RunIIConePlugin(ktR,min_jet_Et)));
00469 #else  // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
00470     is_unavailable("D0RunIICone");
00471 #endif // FASTJET_ENABLE_PLUGIN_D0RUNIICONE
00472   } 
00473   if (all_algs || cmdline.present("-trackjet")) {
00474 #ifdef FASTJET_ENABLE_PLUGIN_TRACKJET
00475     jet_defs.push_back( fj::JetDefinition(new fj::TrackJetPlugin(ktR)));
00476 #else  // FASTJET_ENABLE_PLUGIN_TRACKJET
00477     is_unavailable("TrackJet");
00478 #endif // FASTJET_ENABLE_PLUGIN_TRACKJET
00479   } 
00480   if (all_algs || cmdline.present("-atlascone")) {
00481 #ifdef FASTJET_ENABLE_PLUGIN_ATLASCONE
00482     jet_defs.push_back( fj::JetDefinition(new fj::ATLASConePlugin(ktR)));
00483 #else  // FASTJET_ENABLE_PLUGIN_ATLASCONE
00484     is_unavailable("ATLASCone");
00485 #endif // FASTJET_ENABLE_PLUGIN_ATLASCONE
00486   } 
00487   if (all_algs || cmdline.present("-eecambridge")) {
00488 #ifdef FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
00489     jet_defs.push_back( fj::JetDefinition(new fj::EECambridgePlugin(ycut)));
00490 #else  // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
00491     is_unavailable("EECambridge");
00492 #endif // FASTJET_ENABLE_PLUGIN_EECAMBRIDGE
00493   } 
00494   if (all_algs || cmdline.present("-jade")) {
00495 #ifdef FASTJET_ENABLE_PLUGIN_JADE
00496     jet_defs.push_back( fj::JetDefinition(new fj::JadePlugin()));
00497 #else  // FASTJET_ENABLE_PLUGIN_JADE
00498     is_unavailable("Jade");
00499 #endif // FASTJET_ENABLE_PLUGIN_JADE
00500   } 
00501   if (all_algs || cmdline.present("-cmsiterativecone")) {
00502 #ifdef FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
00503     jet_defs.push_back( fj::JetDefinition(new fj::CMSIterativeConePlugin(ktR,seed_threshold)));
00504 #else  // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
00505     is_unavailable("CMSIterativeCone");
00506 #endif // FASTJET_ENABLE_PLUGIN_CMSITERATIVECONE
00507   } 
00508   if (all_algs || cmdline.present("-d0runipre96cone")) {
00509 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
00510     jet_defs.push_back( fj::JetDefinition(new fj::D0RunIpre96ConePlugin(ktR, seed_threshold, overlap_threshold)));
00511 #else  // FASTJET_ENABLE_PLUGIN_D0RUNICONE
00512     is_unavailable("D0RunIpre96Cone");
00513 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
00514   } 
00515   if (all_algs || cmdline.present("-d0runicone")) {
00516 #ifdef FASTJET_ENABLE_PLUGIN_D0RUNICONE
00517     jet_defs.push_back( fj::JetDefinition(new fj::D0RunIConePlugin(ktR, seed_threshold, overlap_threshold)));
00518 #else  // FASTJET_ENABLE_PLUGIN_D0RUNICONE
00519     is_unavailable("D0RunICone");
00520 #endif // FASTJET_ENABLE_PLUGIN_D0RUNICONE
00521   } 
00522   if (all_algs || cmdline.present("-gridjet")) {
00523 #ifdef FASTJET_ENABLE_PLUGIN_GRIDJET
00524     // we want a grid_ymax of 5.0, but when using R=0.4 (i.e. grid
00525     // spacing of 0.8), this leads to 12.5 grid cells; depending on
00526     // whether this is 12.499999999999 or 12.5000000....1 this gets
00527     // converted either to 12 or 13, making the results sensitive to
00528     // rounding errors.
00529     //
00530     // Instead we therefore take 4.9999999999, which avoids this problem.
00531     double grid_ymax = 4.9999999999;
00532     jet_defs.push_back( fj::JetDefinition(new fj::GridJetPlugin(grid_ymax, ktR*2.0)));
00533 #else  // FASTJET_ENABLE_PLUGIN_GRIDJET
00534     is_unavailable("GridJet");
00535 #endif // FASTJET_ENABLE_PLUGIN_GRIDJET
00536 // end of checking if one asks to run a plugin (don't delete this line)
00537   } 
00538   if (all_algs || 
00539       cmdline.present("-kt") || 
00540       (jet_defs.size() == 0 && !found_unavailable))  {
00541     jet_defs.push_back( fj::JetDefinition(fj::kt_algorithm, ktR, strategy));
00542   }
00543 
00544   string filename = cmdline.value<string>("-file", "");
00545 
00546 
00547   if (!cmdline.all_options_used()) {cerr << 
00548       "Error: some options were not recognized"<<endl; 
00549     exit(-1);}
00550 
00551   for (unsigned idef = 0; idef < jet_defs.size(); idef++) {
00552   fj::JetDefinition & jet_def = jet_defs[idef];
00553   istream * istr;
00554   if (filename == "") istr = &cin;
00555   else                istr = new ifstream(filename.c_str());
00556 
00557   for (int iev = 0; iev < nev; iev++) {
00558   vector<fj::PseudoJet> jets;
00559   vector<fj::PseudoJet> particles;
00560   string line;
00561   int  ndone = 0;
00562   while (getline(*istr, line)) {
00563       //cout << line<<endl;
00564     istringstream linestream(line);
00565     if (line == "#END") {
00566       ndone += 1;
00567       if (ndone == combine) {break;}
00568     }
00569     if (line.substr(0,1) == "#") {continue;}
00570     valarray<double> fourvec(4);
00571     if (hydjet) {
00572       // special reading from hydjet.txt event record (though actually
00573       // this is supposed to be a standard pythia event record, so
00574       // being able to read from it is perhaps not so bad an idea...)
00575       int ii, istat,id,m1,m2,d1,d2;
00576       double mass;
00577       linestream >> ii>> istat >> id >> m1 >> m2 >> d1 >> d2
00578                  >> fourvec[0] >> fourvec[1] >> fourvec[2] >> mass;
00579       // current file contains mass of particle as 4th entry
00580       if (istat == 1) {
00581         fourvec[3] = sqrt(+pow2(fourvec[0])+pow2(fourvec[1])
00582                           +pow2(fourvec[2])+pow2(mass));
00583       }
00584     } else {
00585       if (massless) {
00586         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2];
00587         fourvec[3] = sqrt(pow2(fourvec[0])+pow2(fourvec[1])+pow2(fourvec[2]));}
00588       else {
00589         linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00590       }
00591     }
00592     fj::PseudoJet psjet(fourvec);
00593     if (abs(psjet.rap() < etamax)) {particles.push_back(psjet);}
00594   }
00595 
00596   // add a fake underlying event which is very soft, uniformly distributed
00597   // in eta,phi so as to allow one to reconstruct the area that is associated
00598   // with each jet.
00599   if (add_dense_coverage) {
00600     fj::GhostedAreaSpec ghosted_area_spec(ghost_maxrap);
00601     //fj::GhostedAreaSpec ghosted_area_spec(-2.0,4.0); // asymmetric range
00602     // for plots, reduce the scatter default of 1, to avoid "holes"
00603     // in the subsequent calorimeter view
00604     ghosted_area_spec.set_grid_scatter(0.5); 
00605     ghosted_area_spec.add_ghosts(particles);
00606     //----- old code ------------------
00607     // srand(2);
00608     // int nphi = 60;
00609     // int neta = 100;
00610     // double kt = 1e-1;
00611     // for (int iphi = 0; iphi<nphi; iphi++) {
00612     //   for (int ieta = -neta; ieta<neta+1; ieta++) {
00613     //  double phi = (iphi+0.5) * (fj::twopi/nphi) + rand()*0.001/RAND_MAX;
00614     //  double eta = ieta * (10.0/neta)  + rand()*0.001/RAND_MAX;
00615     //  kt = 1e-20*(1+rand()*0.1/RAND_MAX);
00616     //  double pminus = kt*exp(-eta);
00617     //  double pplus  = kt*exp(+eta);
00618     //  double px = kt*sin(phi);
00619     //  double py = kt*cos(phi);
00620     //  //cout << kt<<" "<<eta<<" "<<phi<<"\n";
00621     //  fj::PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00622     //  particles.push_back(mom);
00623     //   }
00624     // }
00625   }
00626 
00627   // select the particles that pass the selection cut
00628   particles = particles_sel(particles);
00629   
00630   for (int irepeat = 0; irepeat < repeat ; irepeat++) {
00631     int nparticles = particles.size();
00632     try {
00633     auto_ptr<fj::ClusterSequence> clust_seq;
00634     if (do_areas) {
00635       clust_seq.reset(new fj::ClusterSequenceArea(particles,jet_def,area_def));
00636     } else {
00637       clust_seq.reset(new fj::ClusterSequence(particles,jet_def,write));
00638     }
00639 
00640     // repetitive output
00641     if (repeatinclkt >= 0.0) {
00642       vector<fj::PseudoJet> jets_local = sorted_by_pt(clust_seq->inclusive_jets(repeatinclkt));
00643     }
00644 
00645     if (irepeat != 0) {continue;}
00646     cout << "iev "<<iev<< ": number of particles = "<< nparticles << endl;
00647     cout << "strategy used =  "<< clust_seq->strategy_string()<< endl;
00648     if (iev == 0) cout << "Jet Definition: " << jet_def.description() << " (" << fj::fastjet_version_string() << ")" << endl;
00649     if (do_areas && iev == 0) cout << "Area definition: " << area_def.description() << endl;
00650 
00651     // now provide some nice output...
00652     if (inclkt >= 0.0) {
00653       vector<fj::PseudoJet> jets_local = sorted_by_pt(clust_seq->inclusive_jets(inclkt));
00654       print_jets(jets_local, show_constituents);
00655 
00656     }
00657 
00658     if (excln > 0) {
00659       cout << "Printing "<<excln<<" exclusive jets\n";
00660       print_jets(clust_seq->exclusive_jets(excln), show_constituents);
00661     }
00662 
00663     if (excld > 0.0) {
00664       cout << "Printing exclusive jets for d = "<<excld<<"\n";
00665       print_jets(clust_seq->exclusive_jets(excld), show_constituents);
00666     }
00667 
00668     if (excly > 0.0) {
00669       cout << "Printing exclusive jets for ycut = "<<excly<<"\n";
00670       print_jets(clust_seq->exclusive_jets_ycut(excly), show_constituents);
00671     }
00672 
00673     if (get_all_dij) {
00674       for (int i = nparticles-1; i >= 0; i--) {
00675         printf("d for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_dmerge(i));
00676       }
00677     }
00678     if (get_all_yij) {
00679       for (int i = nparticles-1; i >= 0; i--) {
00680         printf("y for n = %4d -> %4d is %14.5e\n", i+1, i, clust_seq->exclusive_ymerge(i));
00681       }
00682     }
00683 
00684     // have the option of printing out the subjets (at scale dcut) of
00685     // each inclusive jet
00686     if (subdcut >= 0.0) {
00687       print_jets_and_sub(clust_seq->inclusive_jets(), subdcut);
00688     }
00689     
00690     // useful for testing that recombination sequences are unique
00691     if (unique_write) {
00692       vector<int> unique_history = clust_seq->unique_history_order();
00693       // construct the inverse of the above mapping
00694       vector<int> inv_unique_history(clust_seq->history().size());
00695       for (unsigned int i = 0; i < unique_history.size(); i++) {
00696         inv_unique_history[unique_history[i]] = i;}
00697 
00698       for (unsigned int i = 0; i < unique_history.size(); i++) {
00699         fj::ClusterSequence::history_element el = 
00700           clust_seq->history()[unique_history[i]];
00701         int uhp1 = el.parent1>=0 ? inv_unique_history[el.parent1] : el.parent1;
00702         int uhp2 = el.parent2>=0 ? inv_unique_history[el.parent2] : el.parent2;
00703         printf("%7d u %15.8e %7d u %7d u\n",i,el.dij,uhp1, uhp2);
00704       }
00705     }
00706 
00707 
00708 #ifdef FASTJET_ENABLE_PLUGIN_SISCONE
00709     // provide some complementary information for SISCone 
00710     if (show_cones) {
00711       const fj::SISConeExtras * extras = 
00712         dynamic_cast<const fj::SISConeExtras *>(clust_seq->extras());
00713       if (extras == 0) 
00714         throw fastjet::Error("extras object for SISCone was null (this can happen with certain area types)");
00715       cout << "most ambiguous split (difference in squared dist) = "
00716            << extras->most_ambiguous_split() << endl;
00717       vector<fastjet::PseudoJet> stable_cones(extras->stable_cones()); 
00718       stable_cones = sorted_by_rapidity(stable_cones);
00719       for (unsigned int i = 0; i < stable_cones.size(); i++) {
00720       //if (stable_cones[i].phi() < 5.0 && stable_cones[i].phi() > 4.0) {
00721         printf("%5u %15.8f %15.8f %15.8f\n",
00722                i,stable_cones[i].rap(),stable_cones[i].phi(),
00723                stable_cones[i].perp() );
00724       //}
00725       }
00726       
00727       // also show passes for jets
00728       vector<fj::PseudoJet> sisjets = clust_seq->inclusive_jets();
00729       printf("\n%15s %15s %15s %12s %8s %8s\n","rap","phi","pt","user-index","pass","nconst");
00730       for (unsigned i = 0; i < sisjets.size(); i++) {
00731         printf("%15.8f %15.8f %15.8f %12d %8d %8u\n",
00732                sisjets[i].rap(), sisjets[i].phi(), sisjets[i].perp(), 
00733                sisjets[i].user_index(), extras->pass(sisjets[i]),
00734                (unsigned int) clust_seq->constituents(sisjets[i]).size()
00735                );
00736         
00737       }
00738     }
00739 #endif // FASTJET_ENABLE_PLUGIN_SISCONE
00740 
00741     if (do_bkgd) {
00742       double rho, sigma, mean_area, empty_area, n_empty_jets;
00743       fj::ClusterSequenceAreaBase * csab = 
00744         dynamic_cast<fj::ClusterSequenceAreaBase *>(clust_seq.get());
00745       if (do_bkgd_csab) {
00746         csab->get_median_rho_and_sigma(bkgd_range, true, rho, sigma, mean_area);
00747         empty_area = csab->empty_area(bkgd_range);
00748         n_empty_jets = csab->n_empty_jets(bkgd_range);
00749       } else if (do_bkgd_jetmedian) {
00750         fj::JetMedianBackgroundEstimator bge(bkgd_range);
00751         bge.set_provide_fj2_sigma(do_bkgd_fj2);
00752         bge.set_cluster_sequence(*csab);
00753         rho = bge.rho();
00754         sigma = bge.sigma();
00755         mean_area = bge.mean_area();
00756         empty_area = bge.empty_area();
00757         n_empty_jets = bge.n_empty_jets();
00758       } else {
00759         assert(do_bkgd_gridmedian);
00760         double rapmin, rapmax;
00761         bkgd_range.get_rapidity_extent(rapmin, rapmax);
00762         fj::GridMedianBackgroundEstimator bge(rapmax, 2*ktR);
00763         bge.set_particles(particles);
00764         rho = bge.rho();
00765         sigma = bge.sigma();
00766         mean_area = bge.mean_area();
00767         empty_area = 0;
00768         n_empty_jets = 0;
00769       }
00770       cout << "  rho = " << rho 
00771            << ", sigma = " << sigma 
00772            << ", mean_area = " << mean_area
00773            << ", empty_area = " << empty_area
00774            << ", n_empty_jets = " << n_empty_jets
00775            << endl;
00776     }
00777   } // try
00778   catch (fastjet::Error fjerr) {
00779     cout << "Caught fastjet error, exiting gracefully" << endl;
00780     exit(0);
00781   }
00782 
00783   } // irepeat
00784   } // iev
00785   // if we've instantiated a plugin, delete it
00786   if (jet_def.strategy()==fj::plugin_strategy){
00787     delete jet_def.plugin();
00788   }
00789   // close any file that we've opened
00790   if (istr != &cin) delete istr;
00791   } // jet_defs
00792 
00793 }
00794 
00795 
00796 
00797 
00798 //------ HELPER ROUTINES -----------------------------------------------
00799 /// print a single jet
00800 void print_jet (const fj::PseudoJet & jet) {
00801   unsigned int n_constituents = jet.constituents().size();
00802   printf("%15.8f %15.8f %15.8f %8u\n",
00803          jet.rap(), jet.phi(), jet.perp(), n_constituents);
00804 }
00805 
00806 
00807 //----------------------------------------------------------------------
00808 void print_jets(const vector<fj::PseudoJet> & jets_in, bool show_constituents) {
00809   vector<fj::PseudoJet> jets;
00810   if (ee_print) {
00811     jets = sorted_by_E(jets_in);
00812     for (unsigned int j = 0; j < jets.size(); j++) {
00813       printf("%5u %15.8f %15.8f %15.8f %15.8f\n",
00814              j,jets[j].px(),jets[j].py(),jets[j].pz(),jets[j].E());
00815       if (show_constituents) {
00816         vector<fj::PseudoJet> const_jets = jets[j].constituents();
00817         for (unsigned int k = 0; k < const_jets.size(); k++) {
00818           printf("        jet%03u %15.8f %15.8f %15.8f %15.8f\n",j,const_jets[k].px(),
00819                  const_jets[k].py(),const_jets[k].pz(),const_jets[k].E());
00820         }
00821         cout << "\n\n";
00822     }
00823 
00824     }
00825   } else {
00826     jets = sorted_by_pt(jets_in);
00827     for (unsigned int j = 0; j < jets.size(); j++) {
00828       printf("%5u %15.8f %15.8f %15.8f",
00829              j,jets[j].rap(),jets[j].phi(),jets[j].perp());
00830       // also print out the scalar area and the perp component of the
00831       // 4-vector (just enough to check a reasonable 4-vector?)
00832       if (do_areas) printf(" %15.8f %15.8f", jets[j].area(),
00833                                              jets[j].area_4vector().perp());
00834       cout << "\n";
00835 
00836       if (show_constituents) {
00837         vector<fj::PseudoJet> const_jets = jets[j].constituents();
00838         for (unsigned int k = 0; k < const_jets.size(); k++) {
00839           printf("        jet%03u %15.8f %15.8f %15.8f %5d\n",j,const_jets[k].rap(),
00840                  const_jets[k].phi(),sqrt(const_jets[k].kt2()), const_jets[k].cluster_hist_index());
00841         }
00842         cout << "\n\n";
00843       }
00844     }
00845   }
00846 
00847   if (rootfile != "") {
00848     ofstream ostr(rootfile.c_str());
00849     ostr << "# " << cmdline_p->command_line() << endl;
00850     ostr << "# output for root" << endl;
00851     assert(jets.size() > 0);
00852     jets[0].validated_cs()->print_jets_for_root(jets,ostr);
00853   }
00854 
00855 }
00856 
00857 
00858 //----- SUBJETS --------------------------------------------------------
00859 /// a function that pretty prints a list of jets and the subjets for each
00860 /// one
00861 void print_jets_and_sub (const vector<fj::PseudoJet> & jets, double dcut) {
00862 
00863   // sort jets into increasing pt
00864   vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);  
00865 
00866   // label the columns
00867   printf("Printing jets and their subjets with subdcut = %10.5f\n",dcut);
00868   printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", 
00869          "phi", "pt", "n constituents");
00870 
00871   // have various kinds of subjet finding, to test consistency among them
00872   enum SubType {internal, newclust_dcut, newclust_R};
00873   SubType subtype = internal;
00874   //SubType subtype = newclust_dcut;
00875   //SubType subtype = newclust_R;
00876 
00877   // print out the details for each jet
00878   //for (unsigned int i = 0; i < sorted_jets.size(); i++) {
00879   for (vector<fj::PseudoJet>::const_iterator jet = sorted_jets.begin(); 
00880        jet != sorted_jets.end(); jet++) {
00881     const fj::JetDefinition & jet_def = jet->validated_cs()->jet_def();
00882 
00883     // if jet pt^2 < dcut with kt alg, then some methods of
00884     // getting subjets will return nothing -- so skip the jet
00885     if (jet_def.jet_algorithm() == fj::kt_algorithm 
00886         && jet->perp2() < dcut) continue;
00887 
00888 
00889     printf("%5u       ",(unsigned int) (jet - sorted_jets.begin()));
00890     print_jet(*jet);
00891     vector<fj::PseudoJet> subjets;
00892     fj::ClusterSequence * cspoint;
00893     if (subtype == internal) {
00894       cspoint = 0;
00895       subjets = jet->exclusive_subjets(dcut);
00896       double ddnp1 = jet->exclusive_subdmerge_max(subjets.size());
00897       double ddn   = jet->exclusive_subdmerge_max(subjets.size()-1);
00898       cout << "     for " << ddnp1 << " < d < " << ddn << " one has " << endl;
00899     } else if (subtype == newclust_dcut) {
00900       cspoint = new fj::ClusterSequence(jet->constituents(), jet_def);
00901       subjets = cspoint->exclusive_jets(dcut);
00902     } else if (subtype == newclust_R) {
00903       assert(jet_def.jet_algorithm() == fj::cambridge_algorithm);
00904       fj::JetDefinition subjd(jet_def.jet_algorithm(), jet_def.R()*sqrt(dcut));
00905       cspoint = new fj::ClusterSequence(jet->constituents(), subjd);
00906       subjets = cspoint->inclusive_jets();
00907     } else {
00908       cerr << "unrecognized subtype for subjet finding" << endl;
00909       exit(-1);
00910     }
00911 
00912     subjets = sorted_by_pt(subjets);
00913     for (unsigned int j = 0; j < subjets.size(); j++) {
00914       printf("    -sub-%02u ",j);
00915       print_jet(subjets[j]);
00916     }
00917 
00918     if (cspoint != 0) delete cspoint;
00919 
00920     //fj::ClusterSequence subseq(clust_seq->constituents(sorted_jets[i]),
00921     //                          fj::JetDefinition(fj::cambridge_algorithm, 0.4));
00922     //vector<fj::PseudoJet> subjets = sorted_by_pt(subseq.inclusive_jets());
00923     //for (unsigned int j = 0; j < subjets.size(); j++) {
00924     //  printf("    -sub-%02u ",j);
00925     //  print_jet(subseq, subjets[j]);
00926     //}
00927   }
00928 
00929 }
00930 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends