FastJet 3.0.2
ClusterSequence.cc
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends