FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: ClusterSequence.cc 2806 2011-12-01 17:21:00Z 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 #include "fastjet/Error.hh" 00030 #include "fastjet/PseudoJet.hh" 00031 #include "fastjet/ClusterSequence.hh" 00032 #include "fastjet/ClusterSequenceStructure.hh" 00033 #include "fastjet/version.hh" // stores the current version number 00034 #include<iostream> 00035 #include<sstream> 00036 #include<fstream> 00037 #include<cmath> 00038 #include<cstdlib> 00039 #include<cassert> 00040 #include<string> 00041 #include<set> 00042 00043 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00044 00045 //---------------------------------------------------------------------- 00046 // here's where we put the main page for fastjet (as explained in the 00047 // Doxygen FAQ) 00048 // We put it inside the fastjet namespace to have the links without 00049 // having to specify (fastjet::) 00050 //...................................................................... 00051 /** \mainpage FastJet code documentation 00052 * 00053 * These pages provide automatically generated documentation for the 00054 * FastJet package. 00055 * 00056 * \section useful_classes The most useful classes 00057 * 00058 * Many of the facilities of FastJet can be accessed through the three 00059 * following classes: 00060 * 00061 * - PseudoJet: the basic class for holding the 4-momentum of a 00062 * particle or a jet. 00063 * 00064 * - JetDefinition: the combination of a #JetAlgorithm and its 00065 * associated parameters. Can also be initialised with a \ref plugins "plugin". 00066 * 00067 * - ClusterSequence: constructed with a vector of input (PseudoJet) 00068 * particles and a JetDefinition, it computes and stores the 00069 * information on how the input particles are clustered into jets. 00070 * 00071 * \section advanced_classes Selected more advanced classes 00072 * 00073 * - ClusterSequenceArea: with the help of an AreaDefinition, provides 00074 * jets that also contain information about their area. 00075 * 00076 * \section Tools Selected additional tools 00077 * 00078 * - JetMedianBackgroundEstimator: with the help of a Selector, a JetDefinition and 00079 * an AreaDefinition, allows one to estimate the background noise density in an event; for a simpler, quicker, effective alternative, use GridMedianBackgroundEstimator 00080 * 00081 * - Transformer: class from which are derived various tools for 00082 * manipulating jets and accessing their substructure. Examples are 00083 * Subtractor, Filter, Pruner and various taggers (e.g. JHTopTagger 00084 * and MassDropTagger). 00085 * 00086 * \section further_info Further information 00087 * 00088 * - Selected classes ordered by topics can be found under the <a 00089 * href="modules.html">modules</a> tab. 00090 * 00091 * - The complete list of classes is available under the <a 00092 * href="annotated.html">classes</a> tab. 00093 * 00094 * - For non-class material (<a href="namespacefastjet.html#enum-members">enums</a>, 00095 * <a href="namespacefastjet.html#typedef-members">typedefs</a>, 00096 * <a href="namespacefastjet.html#func-members">functions</a>), see the 00097 * #fastjet documentation 00098 * 00099 * - For further information and normal documentation, see the main <a 00100 * href="http://fastjet.fr/">FastJet</a> page. 00101 * 00102 * \section examples Examples 00103 * See our \subpage Examples page 00104 */ 00105 00106 // define the doxygen groups 00107 /// \defgroup basic_classes Fundamental FastJet classes 00108 /// \defgroup area_classes Area-related classes 00109 /// \defgroup sec_area_classes Secondary area-related classes 00110 /// \defgroup plugins Plugins for non-native jet definitions 00111 /// \defgroup selectors Selectors 00112 /// \defgroup tools FastJet tools 00113 /// \{ \defgroup tools_generic Generic tools 00114 /// \defgroup tools_background Background subtraction 00115 /// \defgroup tools_taggers Taggers 00116 /// \} 00117 /// \defgroup extra_info Access to extra information 00118 /// \defgroup error_handling Error handling 00119 /// \defgroup advanced_usage Advanced usage 00120 /// \if internal_doc 00121 /// \defgroup internal 00122 /// \endif 00123 00124 //---------------------------------------------------------------------- 00125 00126 00127 using namespace std; 00128 00129 00130 // The following variable can be modified from within user code 00131 // so as to redirect banners to an ostream other than cout. 00132 // 00133 // Please note that if you distribute 3rd party code 00134 // that links with FastJet, that 3rd party code is NOT 00135 // allowed to turn off the printing of FastJet banners 00136 // by default. This requirement reflects the spirit of 00137 // clause 2c of the GNU Public License (v2), under which 00138 // FastJet and its plugins are distributed. 00139 std::ostream * ClusterSequence::_fastjet_banner_ostr = &cout; 00140 00141 00142 // destructor that guarantees proper bookkeeping for the CS Structure 00143 ClusterSequence::~ClusterSequence () { 00144 // set the pointer in the wrapper to this object to NULL to say that 00145 // we're going out of scope 00146 if (_structure_shared_ptr()){ 00147 ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr()); 00148 // normally the csi is purely internal so it really should not be 00149 // NULL i.e assert should be OK 00150 // (we assert rather than throw an error, since failure here is a 00151 // sign of major internal problems) 00152 assert(csi != NULL); 00153 csi->set_associated_cs(NULL); 00154 00155 // if the user had given the CS responsibility to delete itself, 00156 // but then deletes the CS themselves, the following lines of 00157 // code will ensure that the structure_shared_ptr will have 00158 // a proper object count (so that jets associated with the CS will 00159 // throw the correct error if the user tries to access their 00160 // constituents). 00161 if (_deletes_self_when_unused) { 00162 _structure_shared_ptr.set_count(_structure_shared_ptr.use_count() 00163 + _structure_use_count_after_construction); 00164 } 00165 } 00166 } 00167 00168 //----------- 00169 void ClusterSequence::signal_imminent_self_deletion() const { 00170 // normally if the destructor is called when 00171 // _deletes_self_when_unused is true, it assumes that it's been 00172 // called by the user (and it therefore resets the shared pointer 00173 // count to the true count). 00174 // 00175 // for self deletion (called from the destructor of the CSstructure, 00176 // the shared_ptr to which has just had its pointer -> 0) you do 00177 // _not_ want to reset the pointer count (otherwise you will end up 00178 // with a double delete on the shared pointer once you start 00179 // deleting the internal structure of the CS). 00180 // 00181 // the following modification ensures that the count reset will not 00182 // take place in the destructor 00183 assert(_deletes_self_when_unused); 00184 _deletes_self_when_unused = false; 00185 } 00186 00187 //DEP //---------------------------------------------------------------------- 00188 //DEP void ClusterSequence::_initialise_and_run ( 00189 //DEP const double & R, 00190 //DEP const Strategy & strategy, 00191 //DEP const bool & writeout_combinations) { 00192 //DEP 00193 //DEP JetDefinition jet_def(_default_jet_algorithm, R, strategy); 00194 //DEP _initialise_and_run(jet_def, writeout_combinations); 00195 //DEP } 00196 00197 00198 //---------------------------------------------------------------------- 00199 void ClusterSequence::_initialise_and_run ( 00200 const JetDefinition & jet_def_in, 00201 const bool & writeout_combinations) { 00202 00203 // transfer all relevant info into internal variables 00204 _decant_options(jet_def_in, writeout_combinations); 00205 00206 // now run 00207 _initialise_and_run_no_decant(); 00208 } 00209 00210 //---------------------------------------------------------------------- 00211 void ClusterSequence::_initialise_and_run_no_decant () { 00212 00213 // set up the history entries for the initial particles (those 00214 // currently in _jets) 00215 _fill_initial_history(); 00216 00217 // don't run anything if the event is empty 00218 if (n_particles() == 0) return; 00219 00220 // ----- deal with special cases: plugins & e+e- ------ 00221 if (_jet_algorithm == plugin_algorithm) { 00222 // allows plugin_xyz() functions to modify cluster sequence 00223 _plugin_activated = true; 00224 // let the plugin do its work here 00225 _jet_def.plugin()->run_clustering( (*this) ); 00226 _plugin_activated = false; 00227 _update_structure_use_count(); 00228 return; 00229 } else if (_jet_algorithm == ee_kt_algorithm || 00230 _jet_algorithm == ee_genkt_algorithm) { 00231 // ignore requested strategy 00232 _strategy = N2Plain; 00233 if (_jet_algorithm == ee_kt_algorithm) { 00234 // make sure that R is large enough so that "beam" recomb only 00235 // occurs when a single particle is left 00236 // Normally, this should be automatically set to 4 from JetDefinition 00237 assert(_Rparam > 2.0); 00238 // this is used to renormalise the dij to get a "standard" form 00239 // and our convention in e+e- will be different from that 00240 // in long.inv case; NB: _invR2 name should be changed -> _renorm_dij? 00241 _invR2 = 1.0; 00242 } else { 00243 // as of 2009-01-09, choose R to be an angular distance, in 00244 // radians. Since the algorithm uses 2(1-cos(theta)) as its 00245 // squared angular measure, make sure that the _R2 is defined 00246 // in a similar way. 00247 if (_Rparam > pi) { 00248 // choose a value that ensures that back-to-back particles will 00249 // always recombine 00250 //_R2 = 4.0000000000001; 00251 _R2 = 2 * ( 3.0 + cos(_Rparam) ); 00252 } else { 00253 _R2 = 2 * ( 1.0 - cos(_Rparam) ); 00254 } 00255 _invR2 = 1.0/_R2; 00256 } 00257 _simple_N2_cluster_EEBriefJet(); 00258 return; 00259 } else if (_jet_algorithm == undefined_jet_algorithm) { 00260 throw Error("A ClusterSequence cannot be created with an uninitialised JetDefinition"); 00261 } 00262 00263 00264 // automatically redefine the strategy according to N if that is 00265 // what the user requested -- transition points (and especially 00266 // their R-dependence) are based on empirical observations for a 00267 // R=0.4, 0.7 and 1.0, running on toth (3.4GHz, Pentium IV D [dual 00268 // core] with 2MB of cache). 00269 //------------- 00270 // 2011-11-15: lowered N2Plain -> N2Tiled switchover based on some 00271 // new tests on an Intel Core 2 Duo T9400 @ 2.53 GHz 00272 // with 6MB cache; tests performed with lines such as 00273 // ./fastjet_timing_plugins -kt -nhardest 30 -repeat 50000 -strategy -3 -R 0.5 -nev 1 < ../../data/Pythia-PtMin1000-LHC-1000ev.dat 00274 if (_strategy == Best) { 00275 int N = _jets.size(); 00276 //if (N <= 55*max(0.5,min(1.0,_Rparam))) {// old empirical scaling with R 00277 //---------------------- 00278 // 2011-11-15: new empirical scaling with R; NB: low-R N2Tiled 00279 // could be significantly improved at low N by limiting the 00280 // minimum size of tiles when R is small 00281 if (min(1.0,max(0.1,_Rparam)*3.3)*N <= 30) { 00282 _strategy = N2Plain; 00283 } else if (N > 6200/pow(_Rparam,2.0) && _jet_def.jet_algorithm() == cambridge_algorithm) { 00284 _strategy = NlnNCam; 00285 #ifndef DROP_CGAL 00286 } else if ((N > 16000/pow(_Rparam,1.15) && _jet_def.jet_algorithm() != antikt_algorithm) 00287 || N > 35000/pow(_Rparam,1.15)) { 00288 _strategy = NlnN; 00289 #endif // DROP_CGAL 00290 } else if (N <= 450) { 00291 _strategy = N2Tiled; 00292 } else { 00293 _strategy = N2MinHeapTiled; 00294 } 00295 } 00296 00297 // R >= 2pi is not supported by all clustering strategies owing to 00298 // periodicity issues (a particle might cluster with itself). When 00299 // R>=2pi, we therefore automatically switch to a strategy that is 00300 // known to work. 00301 if (_Rparam >= twopi) { 00302 if ( _strategy == NlnN 00303 || _strategy == NlnN3pi 00304 || _strategy == NlnNCam 00305 || _strategy == NlnNCam2pi2R 00306 || _strategy == NlnNCam4pi) { 00307 #ifdef DROP_CGAL 00308 _strategy = N2MinHeapTiled; 00309 #else 00310 _strategy = NlnN4pi; 00311 #endif 00312 } 00313 if (_jet_def.strategy() != Best && _strategy != _jet_def.strategy()) { 00314 ostringstream oss; 00315 oss << "Cluster strategy " << strategy_string(_jet_def.strategy()) 00316 << " automatically changed to " << strategy_string() 00317 << " because the former is not supported for R = " << _Rparam 00318 << " >= 2pi"; 00319 _changed_strategy_warning.warn(oss.str()); 00320 } 00321 } 00322 00323 00324 // run the code containing the selected strategy 00325 // 00326 // We order the strategies stqrting from the ones used by the Best 00327 // strategy in the order of increasing N, then the remaining ones 00328 // again in the order of increasing N. 00329 if (_strategy == N2Plain) { 00330 // BriefJet provides standard long.invariant kt alg. 00331 this->_simple_N2_cluster_BriefJet(); 00332 } else if (_strategy == N2Tiled) { 00333 this->_faster_tiled_N2_cluster(); 00334 } else if (_strategy == N2MinHeapTiled) { 00335 this->_minheap_faster_tiled_N2_cluster(); 00336 } else if (_strategy == NlnN) { 00337 this->_delaunay_cluster(); 00338 } else if (_strategy == NlnNCam) { 00339 this->_CP2DChan_cluster_2piMultD(); 00340 } else if (_strategy == NlnN3pi || _strategy == NlnN4pi ) { 00341 this->_delaunay_cluster(); 00342 } else if (_strategy == N3Dumb ) { 00343 this->_really_dumb_cluster(); 00344 } else if (_strategy == N2PoorTiled) { 00345 this->_tiled_N2_cluster(); 00346 } else if (_strategy == NlnNCam4pi) { 00347 this->_CP2DChan_cluster(); 00348 } else if (_strategy == NlnNCam2pi2R) { 00349 this->_CP2DChan_cluster_2pi2R(); 00350 } else { 00351 ostringstream err; 00352 err << "Unrecognised value for strategy: "<<_strategy; 00353 throw Error(err.str()); 00354 } 00355 00356 } 00357 00358 00359 // these needs to be defined outside the class definition. 00360 bool ClusterSequence::_first_time = true; 00361 int ClusterSequence::_n_exclusive_warnings = 0; 00362 00363 00364 //---------------------------------------------------------------------- 00365 // the version string 00366 string fastjet_version_string() { 00367 return "FastJet version "+string(fastjet_version); 00368 } 00369 00370 00371 //---------------------------------------------------------------------- 00372 // prints a banner on the first call 00373 void ClusterSequence::print_banner() { 00374 00375 if (!_first_time) {return;} 00376 _first_time = false; 00377 00378 // make sure the user has not set the banner stream to NULL 00379 ostream * ostr = _fastjet_banner_ostr; 00380 if (!ostr) return; 00381 00382 (*ostr) << "#--------------------------------------------------------------------------\n"; 00383 (*ostr) << "# FastJet release " << fastjet_version << endl; 00384 (*ostr) << "# M. Cacciari, G.P. Salam and G. Soyez \n"; 00385 (*ostr) << "# A software package for jet finding and analysis at colliders \n"; 00386 (*ostr) << "# http://fastjet.fr \n"; 00387 (*ostr) << "# \n"; 00388 (*ostr) << "# Please cite arXiv:1111.6097 if you use this package for scientific \n"; 00389 (*ostr) << "# work and optionally also Phys. Lett. B641 (2006) [hep-ph/0512210]. \n"; 00390 (*ostr) << "# \n"; 00391 (*ostr) << "# FastJet is provided without warranty under the terms of the GNU GPLv2.\n"; 00392 (*ostr) << "# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code"; 00393 #ifndef DROP_CGAL 00394 (*ostr) << ",\n# CGAL "; 00395 #else 00396 (*ostr) << "\n# "; 00397 #endif // DROP_CGAL 00398 (*ostr) << "and 3rd party plugin jet algorithms. See COPYING file for details.\n"; 00399 (*ostr) << "#--------------------------------------------------------------------------\n"; 00400 // make sure we really have the output done. 00401 ostr->flush(); 00402 } 00403 00404 //---------------------------------------------------------------------- 00405 // transfer all relevant info into internal variables 00406 void ClusterSequence::_decant_options(const JetDefinition & jet_def_in, 00407 const bool & writeout_combinations) { 00408 // make a local copy of the jet definition (for future use) 00409 _jet_def = jet_def_in; 00410 _writeout_combinations = writeout_combinations; 00411 // initialised the wrapper to the current CS 00412 _structure_shared_ptr.reset(new ClusterSequenceStructure(this)); 00413 00414 _decant_options_partial(); 00415 } 00416 00417 //---------------------------------------------------------------------- 00418 // transfer all relevant info into internal variables 00419 void ClusterSequence::_decant_options_partial() { 00420 // let the user know what's going on 00421 print_banner(); 00422 00423 _jet_algorithm = _jet_def.jet_algorithm(); 00424 _Rparam = _jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2; 00425 _strategy = _jet_def.strategy(); 00426 00427 // disallow interference from the plugin 00428 _plugin_activated = false; 00429 00430 // initialised the wrapper to the current CS 00431 //_structure_shared_ptr.reset(new ClusterSequenceStructure(this)); 00432 _update_structure_use_count(); // make sure it's correct already here 00433 } 00434 00435 00436 //---------------------------------------------------------------------- 00437 // initialise the history in a standard way 00438 void ClusterSequence::_fill_initial_history () { 00439 00440 //if (_jets.size() == 0) {throw Error("Cannot run jet-finder on empty event");} 00441 00442 // reserve sufficient space for everything 00443 _jets.reserve(_jets.size()*2); 00444 _history.reserve(_jets.size()*2); 00445 00446 _Qtot = 0; 00447 00448 for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) { 00449 history_element element; 00450 element.parent1 = InexistentParent; 00451 element.parent2 = InexistentParent; 00452 element.child = Invalid; 00453 element.jetp_index = i; 00454 element.dij = 0.0; 00455 element.max_dij_so_far = 0.0; 00456 00457 _history.push_back(element); 00458 00459 // do any momentum preprocessing needed by the recombination scheme 00460 _jet_def.recombiner()->preprocess(_jets[i]); 00461 00462 // get cross-referencing right from PseudoJets 00463 _jets[i].set_cluster_hist_index(i); 00464 _set_structure_shared_ptr(_jets[i]); 00465 00466 // determine the total energy in the event 00467 _Qtot += _jets[i].E(); 00468 } 00469 _initial_n = _jets.size(); 00470 _deletes_self_when_unused = false; 00471 } 00472 00473 00474 //---------------------------------------------------------------------- 00475 // Return the component corresponding to the specified index. 00476 // taken from CLHEP 00477 string ClusterSequence::strategy_string (Strategy strategy_in) const { 00478 string strategy; 00479 switch(strategy_in) { 00480 case NlnN: 00481 strategy = "NlnN"; break; 00482 case NlnN3pi: 00483 strategy = "NlnN3pi"; break; 00484 case NlnN4pi: 00485 strategy = "NlnN4pi"; break; 00486 case N2Plain: 00487 strategy = "N2Plain"; break; 00488 case N2Tiled: 00489 strategy = "N2Tiled"; break; 00490 case N2MinHeapTiled: 00491 strategy = "N2MinHeapTiled"; break; 00492 case N2PoorTiled: 00493 strategy = "N2PoorTiled"; break; 00494 case N3Dumb: 00495 strategy = "N3Dumb"; break; 00496 case NlnNCam4pi: 00497 strategy = "NlnNCam4pi"; break; 00498 case NlnNCam2pi2R: 00499 strategy = "NlnNCam2pi2R"; break; 00500 case NlnNCam: 00501 strategy = "NlnNCam"; break; // 2piMultD 00502 case plugin_strategy: 00503 strategy = "plugin strategy"; break; 00504 default: 00505 strategy = "Unrecognized"; 00506 } 00507 return strategy; 00508 } 00509 00510 00511 double ClusterSequence::jet_scale_for_algorithm( 00512 const PseudoJet & jet) const { 00513 if (_jet_algorithm == kt_algorithm) {return jet.kt2();} 00514 else if (_jet_algorithm == cambridge_algorithm) {return 1.0;} 00515 else if (_jet_algorithm == antikt_algorithm) { 00516 double kt2=jet.kt2(); 00517 return kt2 > 1e-300 ? 1.0/kt2 : 1e300; 00518 } else if (_jet_algorithm == genkt_algorithm) { 00519 double kt2 = jet.kt2(); 00520 double p = jet_def().extra_param(); 00521 if (p <= 0 && kt2 < 1e-300) kt2 = 1e-300; // dodgy safety check 00522 return pow(kt2, p); 00523 } else if (_jet_algorithm == cambridge_for_passive_algorithm) { 00524 double kt2 = jet.kt2(); 00525 double lim = _jet_def.extra_param(); 00526 if (kt2 < lim*lim && kt2 != 0.0) { 00527 return 1.0/kt2; 00528 } else {return 1.0;} 00529 } else {throw Error("Unrecognised jet algorithm");} 00530 } 00531 00532 00533 // //---------------------------------------------------------------------- 00534 // /// transfer the sequence contained in other_seq into our own; 00535 // /// any plugin "extras" contained in the from_seq will be lost 00536 // /// from there. 00537 // void ClusterSequence::transfer_from_sequence(ClusterSequence & from_seq) { 00538 // 00539 // if (will_delete_self_when_unused()) 00540 // throw(Error("cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()")); 00541 // 00542 // // the metadata 00543 // _jet_def = from_seq._jet_def ; 00544 // _writeout_combinations = from_seq._writeout_combinations ; 00545 // _initial_n = from_seq._initial_n ; 00546 // _Rparam = from_seq._Rparam ; 00547 // _R2 = from_seq._R2 ; 00548 // _invR2 = from_seq._invR2 ; 00549 // _strategy = from_seq._strategy ; 00550 // _jet_algorithm = from_seq._jet_algorithm ; 00551 // _plugin_activated = from_seq._plugin_activated ; 00552 // 00553 // // the data 00554 // _jets = from_seq._jets; 00555 // _history = from_seq._history; 00556 // // the following transfers ownership of the extras from the from_seq 00557 // _extras = from_seq._extras; 00558 // 00559 // // transfer of ownership 00560 // if (_structure_shared_ptr()) { 00561 // // anything that is currently associated with the cluster sequence 00562 // // should be told that its cluster sequence no longer exists 00563 // ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr()); 00564 // assert(csi != NULL); 00565 // csi->set_associated_cs(NULL); 00566 // } 00567 // // create a new _structure_shared_ptr to reflect the fact that 00568 // // this CS is essentially a new one 00569 // _structure_shared_ptr.reset(new ClusterSequenceStructure(this)); 00570 // _update_structure_use_count(); 00571 // 00572 // for (vector<PseudoJet>::iterator jit = _jets.begin(); jit != _jets.end(); jit++) 00573 // _set_structure_shared_ptr(*jit); 00574 // } 00575 00576 00577 //---------------------------------------------------------------------- 00578 // transfer the sequence contained in other_seq into our own; 00579 // any plugin "extras" contained in the from_seq will be lost 00580 // from there. 00581 // 00582 // It also sets the ClusterSequence pointers of the PseudoJets in 00583 // the history to point to this ClusterSequence 00584 // 00585 // The second argument is an action that will be applied on every 00586 // jets in the resulting ClusterSequence 00587 void ClusterSequence::transfer_from_sequence(const ClusterSequence & from_seq, 00588 const FunctionOfPseudoJet<PseudoJet> * action_on_jets){ 00589 00590 if (will_delete_self_when_unused()) 00591 throw(Error("cannot use CS::transfer_from_sequence after a call to delete_self_when_unused()")); 00592 00593 // the metadata 00594 _jet_def = from_seq._jet_def ; 00595 _writeout_combinations = from_seq._writeout_combinations ; 00596 _initial_n = from_seq._initial_n ; 00597 _Rparam = from_seq._Rparam ; 00598 _R2 = from_seq._R2 ; 00599 _invR2 = from_seq._invR2 ; 00600 _strategy = from_seq._strategy ; 00601 _jet_algorithm = from_seq._jet_algorithm ; 00602 _plugin_activated = from_seq._plugin_activated ; 00603 00604 // the data 00605 00606 // apply the transformation on the jets if needed 00607 if (action_on_jets) 00608 _jets = (*action_on_jets)(from_seq._jets); 00609 else 00610 _jets = from_seq._jets; 00611 _history = from_seq._history; 00612 // the following shares ownership of the extras with the from_seq; 00613 // no transformations will be applied to the extras 00614 _extras = from_seq._extras; 00615 00616 // clean up existing structure 00617 if (_structure_shared_ptr()) { 00618 // If there are jets associated with an old version of the CS and 00619 // a new one, keeping track of when to delete the CS becomes more 00620 // complex; so we don't allow this situation to occur. 00621 if (_deletes_self_when_unused) throw Error("transfer_from_sequence cannot be used for a cluster sequence that deletes self when unused"); 00622 00623 // anything that is currently associated with the cluster sequence 00624 // should be told that its cluster sequence no longer exists 00625 ClusterSequenceStructure* csi = dynamic_cast<ClusterSequenceStructure*>(_structure_shared_ptr()); 00626 assert(csi != NULL); 00627 csi->set_associated_cs(NULL); 00628 } 00629 // create a new _structure_shared_ptr to reflect the fact that 00630 // this CS is essentially a new one 00631 _structure_shared_ptr.reset(new ClusterSequenceStructure(this)); 00632 _update_structure_use_count(); 00633 00634 for (unsigned int i=0; i<_jets.size(); i++){ 00635 // we reset the cluster history index in case action_on_jets 00636 // messed up with it 00637 _jets[i].set_cluster_hist_index(from_seq._jets[i].cluster_hist_index()); 00638 00639 // reset the structure pointer 00640 _set_structure_shared_ptr(_jets[i]); 00641 } 00642 } 00643 00644 00645 //---------------------------------------------------------------------- 00646 // record an ij recombination and reset the _jets[newjet_k] momentum and 00647 // user index to be those of newjet 00648 void ClusterSequence::plugin_record_ij_recombination( 00649 int jet_i, int jet_j, double dij, 00650 const PseudoJet & newjet, int & newjet_k) { 00651 00652 plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k); 00653 00654 // now transfer newjet into place 00655 int tmp_index = _jets[newjet_k].cluster_hist_index(); 00656 _jets[newjet_k] = newjet; 00657 _jets[newjet_k].set_cluster_hist_index(tmp_index); 00658 _set_structure_shared_ptr(_jets[newjet_k]); 00659 } 00660 00661 00662 //---------------------------------------------------------------------- 00663 // return all inclusive jets with pt > ptmin 00664 vector<PseudoJet> ClusterSequence::inclusive_jets (const double & ptmin) const{ 00665 double dcut = ptmin*ptmin; 00666 int i = _history.size() - 1; // last jet 00667 vector<PseudoJet> jets_local; 00668 if (_jet_algorithm == kt_algorithm) { 00669 while (i >= 0) { 00670 // with our specific definition of dij and diB (i.e. R appears only in 00671 // dij), then dij==diB is the same as the jet.perp2() and we can exploit 00672 // this in selecting the jets... 00673 if (_history[i].max_dij_so_far < dcut) {break;} 00674 if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) { 00675 // for beam jets 00676 int parent1 = _history[i].parent1; 00677 jets_local.push_back(_jets[_history[parent1].jetp_index]);} 00678 i--; 00679 } 00680 } else if (_jet_algorithm == cambridge_algorithm) { 00681 while (i >= 0) { 00682 // inclusive jets are all at end of clustering sequence in the 00683 // Cambridge algorithm -- so if we find a non-exclusive jet, then 00684 // we can exit 00685 if (_history[i].parent2 != BeamJet) {break;} 00686 int parent1 = _history[i].parent1; 00687 const PseudoJet & jet = _jets[_history[parent1].jetp_index]; 00688 if (jet.perp2() >= dcut) {jets_local.push_back(jet);} 00689 i--; 00690 } 00691 } else if (_jet_algorithm == plugin_algorithm 00692 || _jet_algorithm == ee_kt_algorithm 00693 || _jet_algorithm == antikt_algorithm 00694 || _jet_algorithm == genkt_algorithm 00695 || _jet_algorithm == ee_genkt_algorithm 00696 || _jet_algorithm == cambridge_for_passive_algorithm) { 00697 // for inclusive jets with a plugin algorithm, we make no 00698 // assumptions about anything (relation of dij to momenta, 00699 // ordering of the dij, etc.) 00700 while (i >= 0) { 00701 if (_history[i].parent2 == BeamJet) { 00702 int parent1 = _history[i].parent1; 00703 const PseudoJet & jet = _jets[_history[parent1].jetp_index]; 00704 if (jet.perp2() >= dcut) {jets_local.push_back(jet);} 00705 } 00706 i--; 00707 } 00708 } else {throw Error("cs::inclusive_jets(...): Unrecognized jet algorithm");} 00709 return jets_local; 00710 } 00711 00712 00713 //---------------------------------------------------------------------- 00714 // return the number of exclusive jets that would have been obtained 00715 // running the algorithm in exclusive mode with the given dcut 00716 int ClusterSequence::n_exclusive_jets (const double & dcut) const { 00717 00718 // first locate the point where clustering would have stopped (i.e. the 00719 // first time max_dij_so_far > dcut) 00720 int i = _history.size() - 1; // last jet 00721 while (i >= 0) { 00722 if (_history[i].max_dij_so_far <= dcut) {break;} 00723 i--; 00724 } 00725 int stop_point = i + 1; 00726 // relation between stop_point, njets assumes one extra jet disappears 00727 // at each clustering. 00728 int njets = 2*_initial_n - stop_point; 00729 return njets; 00730 } 00731 00732 //---------------------------------------------------------------------- 00733 // return all exclusive jets that would have been obtained running 00734 // the algorithm in exclusive mode with the given dcut 00735 vector<PseudoJet> ClusterSequence::exclusive_jets (const double & dcut) const { 00736 int njets = n_exclusive_jets(dcut); 00737 return exclusive_jets(njets); 00738 } 00739 00740 00741 //---------------------------------------------------------------------- 00742 // return the jets obtained by clustering the event to n jets. 00743 // Throw an error if there are fewer than n particles. 00744 vector<PseudoJet> ClusterSequence::exclusive_jets (const int & njets) const { 00745 00746 // make sure the user does not ask for more than jets than there 00747 // were particles in the first place. 00748 if (njets > _initial_n) { 00749 ostringstream err; 00750 err << "Requested " << njets << " exclusive jets, but there were only " 00751 << _initial_n << " particles in the event"; 00752 throw Error(err.str()); 00753 } 00754 00755 return exclusive_jets_up_to(njets); 00756 } 00757 00758 //---------------------------------------------------------------------- 00759 // return the jets obtained by clustering the event to n jets. 00760 // If there are fewer than n particles, simply return all particles 00761 vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int & njets) const { 00762 00763 // provide a warning when extracting exclusive jets for algorithms 00764 // that does not support it explicitly. 00765 // Native algorithm that support it are: kt, ee_kt, cambridge, 00766 // genkt and ee_genkt (both with p>=0) 00767 // For plugins, we check Plugin::exclusive_sequence_meaningful() 00768 if (( _jet_def.jet_algorithm() != kt_algorithm) && 00769 ( _jet_def.jet_algorithm() != cambridge_algorithm) && 00770 ( _jet_def.jet_algorithm() != ee_kt_algorithm) && 00771 (((_jet_def.jet_algorithm() != genkt_algorithm) && 00772 (_jet_def.jet_algorithm() != ee_genkt_algorithm)) || 00773 (_jet_def.extra_param() <0)) && 00774 ((_jet_def.jet_algorithm() != plugin_algorithm) || 00775 (!_jet_def.plugin()->exclusive_sequence_meaningful())) && 00776 (_n_exclusive_warnings < 5)) { 00777 _n_exclusive_warnings++; 00778 cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl; 00779 } 00780 00781 00782 // calculate the point where we have to stop the clustering. 00783 // relation between stop_point, njets assumes one extra jet disappears 00784 // at each clustering. 00785 int stop_point = 2*_initial_n - njets; 00786 // make sure it's safe when more jets are requested than there are particles 00787 if (stop_point < _initial_n) stop_point = _initial_n; 00788 00789 // some sanity checking to make sure that e+e- does not give us 00790 // surprises (should we ever implement e+e-)... 00791 if (2*_initial_n != static_cast<int>(_history.size())) { 00792 ostringstream err; 00793 err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n"; 00794 throw Error(err.str()); 00795 //assert(false); 00796 } 00797 00798 // now go forwards and reconstitute the jets that we have -- 00799 // basically for any history element, see if the parent jets to 00800 // which it refers were created before the stopping point -- if they 00801 // were then add them to the list, otherwise they are subsequent 00802 // recombinations of the jets that we are looking for. 00803 vector<PseudoJet> jets_local; 00804 for (unsigned int i = stop_point; i < _history.size(); i++) { 00805 int parent1 = _history[i].parent1; 00806 if (parent1 < stop_point) { 00807 jets_local.push_back(_jets[_history[parent1].jetp_index]); 00808 } 00809 int parent2 = _history[i].parent2; 00810 if (parent2 < stop_point && parent2 > 0) { 00811 jets_local.push_back(_jets[_history[parent2].jetp_index]); 00812 } 00813 00814 } 00815 00816 // sanity check... 00817 if (int(jets_local.size()) != min(_initial_n, njets)) { 00818 ostringstream err; 00819 err << "ClusterSequence::exclusive_jets: size of returned vector (" 00820 <<jets_local.size()<<") does not coincide with requested number of jets (" 00821 <<njets<<")"; 00822 throw Error(err.str()); 00823 } 00824 00825 return jets_local; 00826 } 00827 00828 //---------------------------------------------------------------------- 00829 /// return the dmin corresponding to the recombination that went from 00830 /// n+1 to n jets 00831 double ClusterSequence::exclusive_dmerge (const int & njets) const { 00832 assert(njets >= 0); 00833 if (njets >= _initial_n) {return 0.0;} 00834 return _history[2*_initial_n-njets-1].dij; 00835 } 00836 00837 00838 //---------------------------------------------------------------------- 00839 /// return the maximum of the dmin encountered during all recombinations 00840 /// up to the one that led to an n-jet final state; identical to 00841 /// exclusive_dmerge, except in cases where the dmin do not increase 00842 /// monotonically. 00843 double ClusterSequence::exclusive_dmerge_max (const int & njets) const { 00844 assert(njets >= 0); 00845 if (njets >= _initial_n) {return 0.0;} 00846 return _history[2*_initial_n-njets-1].max_dij_so_far; 00847 } 00848 00849 00850 //---------------------------------------------------------------------- 00851 /// return a vector of all subjets of the current jet (in the sense 00852 /// of the exclusive algorithm) that would be obtained when running 00853 /// the algorithm with the given dcut. 00854 std::vector<PseudoJet> ClusterSequence::exclusive_subjets 00855 (const PseudoJet & jet, const double & dcut) const { 00856 00857 set<const history_element*> subhist; 00858 00859 // get the set of history elements that correspond to subjets at 00860 // scale dcut 00861 get_subhist_set(subhist, jet, dcut, 0); 00862 00863 // now transfer this into a sequence of jets 00864 vector<PseudoJet> subjets; 00865 subjets.reserve(subhist.size()); 00866 for (set<const history_element*>::iterator elem = subhist.begin(); 00867 elem != subhist.end(); elem++) { 00868 subjets.push_back(_jets[(*elem)->jetp_index]); 00869 } 00870 return subjets; 00871 } 00872 00873 //---------------------------------------------------------------------- 00874 /// return the size of exclusive_subjets(...); still n ln n with same 00875 /// coefficient, but marginally more efficient than manually taking 00876 /// exclusive_subjets.size() 00877 int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet, 00878 const double & dcut) const { 00879 set<const history_element*> subhist; 00880 // get the set of history elements that correspond to subjets at 00881 // scale dcut 00882 get_subhist_set(subhist, jet, dcut, 0); 00883 return subhist.size(); 00884 } 00885 00886 //---------------------------------------------------------------------- 00887 /// return the list of subjets obtained by unclustering the supplied 00888 /// jet down to nsub subjets. Throws an error if there are fewer than 00889 /// nsub particles in the jet. 00890 std::vector<PseudoJet> ClusterSequence::exclusive_subjets 00891 (const PseudoJet & jet, int nsub) const { 00892 vector<PseudoJet> subjets = exclusive_subjets_up_to(jet, nsub); 00893 if (int(subjets.size()) < nsub) { 00894 ostringstream err; 00895 err << "Requested " << nsub << " exclusive subjets, but there were only " 00896 << subjets.size() << " particles in the jet"; 00897 throw Error(err.str()); 00898 } 00899 return subjets; 00900 00901 } 00902 00903 //---------------------------------------------------------------------- 00904 /// return the list of subjets obtained by unclustering the supplied 00905 /// jet down to nsub subjets (or all constituents if there are fewer 00906 /// than nsub). 00907 std::vector<PseudoJet> ClusterSequence::exclusive_subjets_up_to 00908 (const PseudoJet & jet, int nsub) const { 00909 00910 set<const history_element*> subhist; 00911 00912 // prepare the vector into which we'll put the result 00913 vector<PseudoJet> subjets; 00914 if (nsub < 0) throw Error("Requested a negative number of subjets. This is nonsensical."); 00915 if (nsub == 0) return subjets; 00916 00917 // get the set of history elements that correspond to subjets at 00918 // scale dcut 00919 get_subhist_set(subhist, jet, -1.0, nsub); 00920 00921 // now transfer this into a sequence of jets 00922 subjets.reserve(subhist.size()); 00923 for (set<const history_element*>::iterator elem = subhist.begin(); 00924 elem != subhist.end(); elem++) { 00925 subjets.push_back(_jets[(*elem)->jetp_index]); 00926 } 00927 return subjets; 00928 } 00929 00930 00931 //---------------------------------------------------------------------- 00932 /// return the dij that was present in the merging nsub+1 -> nsub 00933 /// subjets inside this jet. 00934 /// 00935 /// If the jet has nsub or fewer constituents, it will return 0. 00936 double ClusterSequence::exclusive_subdmerge(const PseudoJet & jet, int nsub) const { 00937 set<const history_element*> subhist; 00938 00939 // get the set of history elements that correspond to subjets at 00940 // scale dcut 00941 get_subhist_set(subhist, jet, -1.0, nsub); 00942 00943 set<const history_element*>::iterator highest = subhist.end(); 00944 highest--; 00945 /// will be zero if nconst <= nsub, since highest will be an original 00946 /// particle have zero dij 00947 return (*highest)->dij; 00948 } 00949 00950 00951 //---------------------------------------------------------------------- 00952 /// return the maximum dij that occurred in the whole event at the 00953 /// stage that the nsub+1 -> nsub merge of subjets occurred inside 00954 /// this jet. 00955 /// 00956 /// If the jet has nsub or fewer constituents, it will return 0. 00957 double ClusterSequence::exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const { 00958 00959 set<const history_element*> subhist; 00960 00961 // get the set of history elements that correspond to subjets at 00962 // scale dcut 00963 get_subhist_set(subhist, jet, -1.0, nsub); 00964 00965 set<const history_element*>::iterator highest = subhist.end(); 00966 highest--; 00967 /// will be zero if nconst <= nsub, since highest will be an original 00968 /// particle have zero dij 00969 return (*highest)->max_dij_so_far; 00970 } 00971 00972 00973 00974 //---------------------------------------------------------------------- 00975 /// return a set of pointers to history entries corresponding to the 00976 /// subjets of this jet; one stops going working down through the 00977 /// subjets either when 00978 /// - there is no further to go 00979 /// - one has found maxjet entries 00980 /// - max_dij_so_far <= dcut 00981 void ClusterSequence::get_subhist_set(set<const history_element*> & subhist, 00982 const PseudoJet & jet, 00983 double dcut, int maxjet) const { 00984 assert(contains(jet)); 00985 00986 subhist.clear(); 00987 subhist.insert(&(_history[jet.cluster_hist_index()])); 00988 00989 // establish the set of jets that are relevant 00990 int njet = 1; 00991 while (true) { 00992 // first find out if we need to probe deeper into jet. 00993 // Get history element closest to end of sequence 00994 set<const history_element*>::iterator highest = subhist.end(); 00995 assert (highest != subhist.begin()); 00996 highest--; 00997 const history_element* elem = *highest; 00998 // make sure we haven't got too many jets 00999 if (njet == maxjet) break; 01000 // make sure it has parents 01001 if (elem->parent1 < 0) break; 01002 // make sure that we still resolve it at scale dcut 01003 if (elem->max_dij_so_far <= dcut) break; 01004 01005 // then do so: replace "highest" with its two parents 01006 subhist.erase(highest); 01007 subhist.insert(&(_history[elem->parent1])); 01008 subhist.insert(&(_history[elem->parent2])); 01009 njet++; 01010 } 01011 } 01012 01013 //---------------------------------------------------------------------- 01014 // work through the object's history until 01015 bool ClusterSequence::object_in_jet(const PseudoJet & object, 01016 const PseudoJet & jet) const { 01017 01018 // make sure the object conceivably belongs to this clustering 01019 // sequence 01020 assert(contains(object) && contains(jet)); 01021 01022 const PseudoJet * this_object = &object; 01023 const PseudoJet * childp; 01024 while(true) { 01025 if (this_object->cluster_hist_index() == jet.cluster_hist_index()) { 01026 return true; 01027 } else if (has_child(*this_object, childp)) { 01028 this_object = childp; 01029 } else { 01030 return false; 01031 } 01032 } 01033 } 01034 01035 //---------------------------------------------------------------------- 01036 /// if the jet has parents in the clustering, it returns true 01037 /// and sets parent1 and parent2 equal to them. 01038 /// 01039 /// if it has no parents it returns false and sets parent1 and 01040 /// parent2 to zero 01041 bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1, 01042 PseudoJet & parent2) const { 01043 01044 const history_element & hist = _history[jet.cluster_hist_index()]; 01045 01046 // make sure we do not run into any unexpected situations -- 01047 // i.e. both parents valid, or neither 01048 assert ((hist.parent1 >= 0 && hist.parent2 >= 0) || 01049 (hist.parent1 < 0 && hist.parent2 < 0)); 01050 01051 if (hist.parent1 < 0) { 01052 parent1 = PseudoJet(0.0,0.0,0.0,0.0); 01053 parent2 = parent1; 01054 return false; 01055 } else { 01056 parent1 = _jets[_history[hist.parent1].jetp_index]; 01057 parent2 = _jets[_history[hist.parent2].jetp_index]; 01058 // order the parents in decreasing pt 01059 if (parent1.perp2() < parent2.perp2()) std::swap(parent1,parent2); 01060 return true; 01061 } 01062 } 01063 01064 //---------------------------------------------------------------------- 01065 /// if the jet has a child then return true and give the child jet 01066 /// otherwise return false and set the child to zero 01067 bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const { 01068 01069 //const history_element & hist = _history[jet.cluster_hist_index()]; 01070 // 01071 //if (hist.child >= 0) { 01072 // child = _jets[_history[hist.child].jetp_index]; 01073 // return true; 01074 //} else { 01075 // child = PseudoJet(0.0,0.0,0.0,0.0); 01076 // return false; 01077 //} 01078 const PseudoJet * childp; 01079 bool res = has_child(jet, childp); 01080 if (res) { 01081 child = *childp; 01082 return true; 01083 } else { 01084 child = PseudoJet(0.0,0.0,0.0,0.0); 01085 return false; 01086 } 01087 } 01088 01089 bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const { 01090 01091 const history_element & hist = _history[jet.cluster_hist_index()]; 01092 01093 // check that this jet has a child and that the child corresponds to 01094 // a true jet [RETHINK-IF-CHANGE-NUMBERING: what is the right 01095 // behaviour if the child is the same jet but made inclusive...?] 01096 if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) { 01097 childp = &(_jets[_history[hist.child].jetp_index]); 01098 return true; 01099 } else { 01100 childp = NULL; 01101 return false; 01102 } 01103 } 01104 01105 01106 //---------------------------------------------------------------------- 01107 /// if this jet has a child (and so a partner) return true 01108 /// and give the partner, otherwise return false and set the 01109 /// partner to zero 01110 bool ClusterSequence::has_partner(const PseudoJet & jet, 01111 PseudoJet & partner) const { 01112 01113 const history_element & hist = _history[jet.cluster_hist_index()]; 01114 01115 // make sure we have a child and that the child does not correspond 01116 // to a clustering with the beam (or some other invalid quantity) 01117 if (hist.child >= 0 && _history[hist.child].parent2 >= 0) { 01118 const history_element & child_hist = _history[hist.child]; 01119 if (child_hist.parent1 == jet.cluster_hist_index()) { 01120 // partner will be child's parent2 -- for iB clustering 01121 // parent2 will not be valid 01122 partner = _jets[_history[child_hist.parent2].jetp_index]; 01123 } else { 01124 // partner will be child's parent1 01125 partner = _jets[_history[child_hist.parent1].jetp_index]; 01126 } 01127 return true; 01128 } else { 01129 partner = PseudoJet(0.0,0.0,0.0,0.0); 01130 return false; 01131 } 01132 } 01133 01134 01135 //---------------------------------------------------------------------- 01136 // return a vector of the particles that make up a jet 01137 vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const { 01138 vector<PseudoJet> subjets; 01139 add_constituents(jet, subjets); 01140 return subjets; 01141 } 01142 01143 //---------------------------------------------------------------------- 01144 /// output the supplied vector of jets in a format that can be read 01145 /// by an appropriate root script; the format is: 01146 /// jet-n jet-px jet-py jet-pz jet-E 01147 /// particle-n particle-rap particle-phi particle-pt 01148 /// particle-n particle-rap particle-phi particle-pt 01149 /// ... 01150 /// #END 01151 /// ... [i.e. above repeated] 01152 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in, 01153 ostream & ostr) const { 01154 for (unsigned i = 0; i < jets_in.size(); i++) { 01155 ostr << i << " " 01156 << jets_in[i].px() << " " 01157 << jets_in[i].py() << " " 01158 << jets_in[i].pz() << " " 01159 << jets_in[i].E() << endl; 01160 vector<PseudoJet> cst = constituents(jets_in[i]); 01161 for (unsigned j = 0; j < cst.size() ; j++) { 01162 ostr << " " << j << " " 01163 << cst[j].rap() << " " 01164 << cst[j].phi() << " " 01165 << cst[j].perp() << endl; 01166 } 01167 ostr << "#END" << endl; 01168 } 01169 } 01170 01171 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets_in, 01172 const std::string & filename, 01173 const std::string & comment ) const { 01174 std::ofstream ostr(filename.c_str()); 01175 if (comment != "") ostr << "# " << comment << endl; 01176 print_jets_for_root(jets_in, ostr); 01177 } 01178 01179 01180 // Not yet. Perhaps in a future release 01181 // //---------------------------------------------------------------------- 01182 // // print out all inclusive jets with pt > ptmin 01183 // void ClusterSequence::print_jets (const double & ptmin) const{ 01184 // vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin)); 01185 // 01186 // for (size_t j = 0; j < jets.size(); j++) { 01187 // printf("%5u %7.3f %7.3f %9.3f\n", 01188 // j,jets[j].rap(),jets[j].phi(),jets[j].perp()); 01189 // } 01190 // } 01191 01192 //---------------------------------------------------------------------- 01193 /// returns a vector of size n_particles() which indicates, for 01194 /// each of the initial particles (in the order in which they were 01195 /// supplied), which of the supplied jets it belongs to; if it does 01196 /// not belong to any of the supplied jets, the index is set to -1; 01197 vector<int> ClusterSequence::particle_jet_indices( 01198 const vector<PseudoJet> & jets_in) const { 01199 01200 vector<int> indices(n_particles()); 01201 01202 // first label all particles as not belonging to any jets 01203 for (unsigned ipart = 0; ipart < n_particles(); ipart++) 01204 indices[ipart] = -1; 01205 01206 // then for each of the jets relabel its consituents as belonging to 01207 // that jet 01208 for (unsigned ijet = 0; ijet < jets_in.size(); ijet++) { 01209 01210 vector<PseudoJet> jet_constituents(constituents(jets_in[ijet])); 01211 01212 for (unsigned ip = 0; ip < jet_constituents.size(); ip++) { 01213 // a safe (if slightly redundant) way of getting the particle 01214 // index (for initial particles it is actually safe to assume 01215 // ipart=iclust). 01216 unsigned iclust = jet_constituents[ip].cluster_hist_index(); 01217 unsigned ipart = history()[iclust].jetp_index; 01218 indices[ipart] = ijet; 01219 } 01220 } 01221 01222 return indices; 01223 } 01224 01225 01226 //---------------------------------------------------------------------- 01227 // recursive routine that adds on constituents of jet to the subjet_vector 01228 void ClusterSequence::add_constituents ( 01229 const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const { 01230 // find out position in cluster history 01231 int i = jet.cluster_hist_index(); 01232 int parent1 = _history[i].parent1; 01233 int parent2 = _history[i].parent2; 01234 01235 if (parent1 == InexistentParent) { 01236 // It is an original particle (labelled by its parent having value 01237 // InexistentParent), therefore add it on to the subjet vector 01238 // Note: we add the initial particle and not simply 'jet' so that 01239 // calling add_constituents with a subtracted jet containing 01240 // only one particle will work. 01241 subjet_vector.push_back(_jets[i]); 01242 return; 01243 } 01244 01245 // add parent 1 01246 add_constituents(_jets[_history[parent1].jetp_index], subjet_vector); 01247 01248 // see if parent2 is a real jet; if it is then add its constituents 01249 if (parent2 != BeamJet) { 01250 add_constituents(_jets[_history[parent2].jetp_index], subjet_vector); 01251 } 01252 } 01253 01254 01255 01256 //---------------------------------------------------------------------- 01257 // initialise the history in a standard way 01258 void ClusterSequence::_add_step_to_history ( 01259 const int & step_number, const int & parent1, 01260 const int & parent2, const int & jetp_index, 01261 const double & dij) { 01262 01263 history_element element; 01264 element.parent1 = parent1; 01265 element.parent2 = parent2; 01266 element.jetp_index = jetp_index; 01267 element.child = Invalid; 01268 element.dij = dij; 01269 element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far); 01270 _history.push_back(element); 01271 01272 int local_step = _history.size()-1; 01273 assert(local_step == step_number); 01274 01275 assert(parent1 >= 0); 01276 _history[parent1].child = local_step; 01277 if (parent2 >= 0) {_history[parent2].child = local_step;} 01278 01279 // get cross-referencing right from PseudoJets 01280 if (jetp_index != Invalid) { 01281 assert(jetp_index >= 0); 01282 //cout << _jets.size() <<" "<<jetp_index<<"\n"; 01283 _jets[jetp_index].set_cluster_hist_index(local_step); 01284 _set_structure_shared_ptr(_jets[jetp_index]); 01285 } 01286 01287 if (_writeout_combinations) { 01288 cout << local_step << ": " 01289 << parent1 << " with " << parent2 01290 << "; y = "<< dij<<endl; 01291 } 01292 01293 } 01294 01295 01296 01297 01298 //====================================================================== 01299 // Return an order in which to read the history such that _history[order[i]] 01300 // will always correspond to the same set of consituent particles if 01301 // two branching histories are equivalent in terms of the particles 01302 // contained in any given pseudojet. 01303 vector<int> ClusterSequence::unique_history_order() const { 01304 01305 // first construct an array that will tell us the lowest constituent 01306 // of a given jet -- this will always be one of the original 01307 // particles, whose order is well defined and so will help us to 01308 // follow the tree in a unique manner. 01309 valarray<int> lowest_constituent(_history.size()); 01310 int hist_n = _history.size(); 01311 lowest_constituent = hist_n; // give it a large number 01312 for (int i = 0; i < hist_n; i++) { 01313 // sets things up for the initial partons 01314 lowest_constituent[i] = min(lowest_constituent[i],i); 01315 // propagates them through to the children of this parton 01316 if (_history[i].child > 0) lowest_constituent[_history[i].child] 01317 = min(lowest_constituent[_history[i].child],lowest_constituent[i]); 01318 } 01319 01320 // establish an array for what we have and have not extracted so far 01321 valarray<bool> extracted(_history.size()); extracted = false; 01322 vector<int> unique_tree; 01323 unique_tree.reserve(_history.size()); 01324 01325 // now work our way through the tree 01326 for (unsigned i = 0; i < n_particles(); i++) { 01327 if (!extracted[i]) { 01328 unique_tree.push_back(i); 01329 extracted[i] = true; 01330 _extract_tree_children(i, extracted, lowest_constituent, unique_tree); 01331 } 01332 } 01333 01334 return unique_tree; 01335 } 01336 01337 //====================================================================== 01338 // helper for unique_history_order 01339 void ClusterSequence::_extract_tree_children( 01340 int position, 01341 valarray<bool> & extracted, 01342 const valarray<int> & lowest_constituent, 01343 vector<int> & unique_tree) const { 01344 if (!extracted[position]) { 01345 // that means we may have unidentified parents around, so go and 01346 // collect them (extracted[position]) will then be made true) 01347 _extract_tree_parents(position,extracted,lowest_constituent,unique_tree); 01348 } 01349 01350 // now look after the children... 01351 int child = _history[position].child; 01352 if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree); 01353 } 01354 01355 01356 //====================================================================== 01357 // return the list of unclustered particles 01358 vector<PseudoJet> ClusterSequence::unclustered_particles() const { 01359 vector<PseudoJet> unclustered; 01360 for (unsigned i = 0; i < n_particles() ; i++) { 01361 if (_history[i].child == Invalid) 01362 unclustered.push_back(_jets[_history[i].jetp_index]); 01363 } 01364 return unclustered; 01365 } 01366 01367 //====================================================================== 01368 /// Return the list of pseudojets in the ClusterSequence that do not 01369 /// have children (and are not among the inclusive jets). They may 01370 /// result from a clustering step or may be one of the pseudojets 01371 /// returned by unclustered_particles(). 01372 vector<PseudoJet> ClusterSequence::childless_pseudojets() const { 01373 vector<PseudoJet> unclustered; 01374 for (unsigned i = 0; i < _history.size() ; i++) { 01375 if ((_history[i].child == Invalid) && (_history[i].parent2 != BeamJet)) 01376 unclustered.push_back(_jets[_history[i].jetp_index]); 01377 } 01378 return unclustered; 01379 } 01380 01381 01382 01383 //---------------------------------------------------------------------- 01384 // returns true if the cluster sequence contains this jet (i.e. jet's 01385 // structure is this cluster sequence's and the cluster history index 01386 // is in a consistent range) 01387 bool ClusterSequence::contains(const PseudoJet & jet) const { 01388 return jet.cluster_hist_index() >= 0 01389 && jet.cluster_hist_index() < int(_history.size()) 01390 && jet.has_valid_cluster_sequence() 01391 && jet.associated_cluster_sequence() == this; 01392 } 01393 01394 01395 01396 //====================================================================== 01397 // helper for unique_history_order 01398 void ClusterSequence::_extract_tree_parents( 01399 int position, 01400 valarray<bool> & extracted, 01401 const valarray<int> & lowest_constituent, 01402 vector<int> & unique_tree) const { 01403 01404 if (!extracted[position]) { 01405 int parent1 = _history[position].parent1; 01406 int parent2 = _history[position].parent2; 01407 // where relevant order parents so that we will first treat the 01408 // one containing the smaller "lowest_constituent" 01409 if (parent1 >= 0 && parent2 >= 0) { 01410 if (lowest_constituent[parent1] > lowest_constituent[parent2]) 01411 std::swap(parent1, parent2); 01412 } 01413 // then actually run through the parents to extract the constituents... 01414 if (parent1 >= 0 && !extracted[parent1]) 01415 _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree); 01416 if (parent2 >= 0 && !extracted[parent2]) 01417 _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree); 01418 // finally declare this position to be accounted for and push it 01419 // onto our list. 01420 unique_tree.push_back(position); 01421 extracted[position] = true; 01422 } 01423 } 01424 01425 01426 //====================================================================== 01427 /// carries out the bookkeeping associated with the step of recombining 01428 /// jet_i and jet_j (assuming a distance dij) and returns the index 01429 /// of the recombined jet, newjet_k. 01430 void ClusterSequence::_do_ij_recombination_step( 01431 const int & jet_i, const int & jet_j, 01432 const double & dij, 01433 int & newjet_k) { 01434 01435 // Create the new jet by recombining the first two. 01436 // 01437 // For efficiency reasons, use a ctr that initialises only the 01438 // shared pointers, since the rest of the info will anyway be dealt 01439 // with by the recombiner. 01440 PseudoJet newjet(false); 01441 _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet); 01442 _jets.push_back(newjet); 01443 // original version... 01444 //_jets.push_back(_jets[jet_i] + _jets[jet_j]); 01445 01446 // get its index 01447 newjet_k = _jets.size()-1; 01448 01449 // get history index 01450 int newstep_k = _history.size(); 01451 // and provide jet with the info 01452 _jets[newjet_k].set_cluster_hist_index(newstep_k); 01453 01454 // finally sort out the history 01455 int hist_i = _jets[jet_i].cluster_hist_index(); 01456 int hist_j = _jets[jet_j].cluster_hist_index(); 01457 01458 _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j), 01459 newjet_k, dij); 01460 01461 } 01462 01463 01464 //====================================================================== 01465 /// carries out the bookkeeping associated with the step of recombining 01466 /// jet_i with the beam 01467 void ClusterSequence::_do_iB_recombination_step( 01468 const int & jet_i, const double & diB) { 01469 // get history index 01470 int newstep_k = _history.size(); 01471 01472 // recombine the jet with the beam 01473 _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet, 01474 Invalid, diB); 01475 01476 } 01477 01478 01479 // make sure the static member _changed_strategy_warning is defined. 01480 LimitedWarning ClusterSequence::_changed_strategy_warning; 01481 01482 01483 //---------------------------------------------------------------------- 01484 void ClusterSequence::_set_structure_shared_ptr(PseudoJet & j) { 01485 j.set_structure_shared_ptr(_structure_shared_ptr); 01486 // record the use count of the structure shared point to help 01487 // in case we want to ask the CS to handle its own memory 01488 _update_structure_use_count(); 01489 } 01490 01491 01492 //---------------------------------------------------------------------- 01493 void ClusterSequence::_update_structure_use_count() { 01494 // record the use count of the structure shared point to help 01495 // in case we want to ask the CS to handle its own memory 01496 _structure_use_count_after_construction = _structure_shared_ptr.use_count(); 01497 } 01498 01499 //---------------------------------------------------------------------- 01500 /// by calling this routine you tell the ClusterSequence to delete 01501 /// itself when all the Pseudojets associated with it have gone out 01502 /// of scope. 01503 void ClusterSequence::delete_self_when_unused() { 01504 // the trick we use to handle this is to modify the use count; 01505 // that way the structure will be deleted when there are no external 01506 // objects left associated the CS and the structure's destructor will then 01507 // look after deleting the cluster sequence 01508 01509 // first make sure that there is at least one other object 01510 // associated with the CS 01511 int new_count = _structure_shared_ptr.use_count() - _structure_use_count_after_construction; 01512 if (new_count <= 0) { 01513 throw Error("delete_self_when_unused may only be called if at least one object outside the CS (e.g. a jet) is already associated with the CS"); 01514 } 01515 01516 _structure_shared_ptr.set_count(new_count); 01517 _deletes_self_when_unused = true; 01518 } 01519 01520 01521 FASTJET_END_NAMESPACE 01522