FastJet 3.0.6
|
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