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