ClusterSequence.cc

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: ClusterSequence.cc 960 2007-11-29 17:03:47Z cacciari $
00003 //
00004 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
00026 //  Foundation, Inc.:
00027 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00028 //----------------------------------------------------------------------
00029 //ENDHEADER
00030 
00031 #include "fastjet/Error.hh"
00032 #include "fastjet/PseudoJet.hh"
00033 #include "fastjet/ClusterSequence.hh"
00034 #include "fastjet/version.hh" // stores the current version number
00035 #include<iostream>
00036 #include<sstream>
00037 #include<cmath>
00038 #include<cstdlib>
00039 #include<cassert>
00040 #include<string>
00041 
00042 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00043 
00044 using namespace std;
00045 
00047 JetAlgorithm ClusterSequence::_default_jet_algorithm = kt_algorithm;
00048 //
00049 
00050 
00051 //----------------------------------------------------------------------
00052 void ClusterSequence::_initialise_and_run (
00053                                   const double & R,
00054                                   const Strategy & strategy,
00055                                   const bool & writeout_combinations) {
00056 
00057   JetDefinition jet_def(_default_jet_algorithm, R, strategy);
00058   _initialise_and_run(jet_def, writeout_combinations);
00059 }
00060 
00061 
00062 //----------------------------------------------------------------------
00063 void ClusterSequence::_initialise_and_run (
00064                                   const JetDefinition & jet_def,
00065                                   const bool & writeout_combinations) {
00066 
00067   // transfer all relevant info into internal variables
00068   _decant_options(jet_def, writeout_combinations);
00069 
00070   // set up the history entries for the initial particles (those
00071   // currently in _jets)
00072   _fill_initial_history();
00073 
00074   // run the plugin if that's what's decreed
00075   if (_jet_algorithm == plugin_algorithm) {
00076     _plugin_activated = true;
00077     _jet_def.plugin()->run_clustering( (*this) );
00078     _plugin_activated = false;
00079     return;
00080   }
00081 
00082 
00083   // automatically redefine the strategy according to N if that is
00084   // what the user requested -- transition points (and especially
00085   // their R-dependence) are based on empirical observations for a
00086   // R=0.4, 0.7 and 1.0, running on toth (3.4GHz, Pentium IV D [dual
00087   // core] with 2MB of cache).
00088   if (_strategy == Best) {
00089     int N = _jets.size();
00090     if (N > 6200/pow(_Rparam,2.0) 
00091         && jet_def.jet_algorithm() == cambridge_algorithm) {
00092       _strategy = NlnNCam;}
00093     else
00094 #ifndef DROP_CGAL
00095     if (N > 16000/pow(_Rparam,1.15)) {
00096       _strategy = NlnN; }   
00097     else                    
00098 #endif  // DROP_CGAL
00099       if (N > 450) {
00100       _strategy = N2MinHeapTiled;
00101     }
00102     else if (N > 55*max(0.5,min(1.0,_Rparam))) {// empirical scaling with R
00103       _strategy = N2Tiled;
00104     } else {
00105       _strategy = N2Plain;
00106     }
00107   }
00108 
00109 
00110   // run the code containing the selected strategy
00111   if (_strategy == NlnN || _strategy == NlnN3pi 
00112       || _strategy == NlnN4pi ) {
00113     this->_delaunay_cluster();
00114   } else if (_strategy ==  N3Dumb ) {
00115     this->_really_dumb_cluster();
00116   } else if (_strategy == N2Tiled) {
00117     this->_faster_tiled_N2_cluster();
00118   } else if (_strategy == N2PoorTiled) {
00119     this->_tiled_N2_cluster();
00120   } else if (_strategy == N2Plain) {
00121     this->_simple_N2_cluster();
00122   } else if (_strategy == N2MinHeapTiled) {
00123     this->_minheap_faster_tiled_N2_cluster();
00124   } else if (_strategy == NlnNCam4pi) {
00125     this->_CP2DChan_cluster();
00126   } else if (_strategy == NlnNCam2pi2R) {
00127     this->_CP2DChan_cluster_2pi2R();
00128   } else if (_strategy == NlnNCam) {
00129     this->_CP2DChan_cluster_2piMultD();
00130   } else {
00131     ostringstream err;
00132     err << "Unrecognised value for strategy: "<<_strategy;
00133     throw Error(err.str());
00134     //assert(false);
00135   }
00136 }
00137 
00138 
00139 // these needs to be defined outside the class definition.
00140 bool ClusterSequence::_first_time = true;
00141 int ClusterSequence::_n_exclusive_warnings = 0;
00142 
00143 
00144 //----------------------------------------------------------------------
00145 // the version string
00146 string fastjet_version_string() {
00147   return "FastJet version "+string(fastjet_version);
00148 }
00149 
00150 
00151 //----------------------------------------------------------------------
00152 // prints a banner on the first call
00153 void ClusterSequence::_print_banner() {
00154 
00155   if (!_first_time) {return;}
00156   _first_time = false;
00157   
00158   
00159   //Symp. Discr. Alg, p.472 (2002) and  CGAL (http://www.cgal.org);
00160 
00161   cout << "#--------------------------------------------------------------------------\n";
00162   cout << "#                      FastJet release " << fastjet_version << endl;
00163   cout << "#            Written by M. Cacciari, G.P. Salam and G. Soyez            \n"; 
00164   cout << "#              http://www.lpthe.jussieu.fr/~salam/fastjet               \n"; 
00165   cout << "#                                                                       \n";
00166   cout << "# Longitudinally invariant Kt, anti-Kt, and inclusive Cambridge/Aachen  \n";
00167   cout << "# clustering using fast geometric algorithms, with area measures and optional\n";
00168   cout << "# external jet-finder plugins.                                          \n";
00169   cout << "# Please cite Phys. Lett. B641 (2006) [hep-ph/0512210] if you use this code.\n";
00170   cout << "#                                                                       \n";
00171   cout << "# This package uses T.Chan's closest pair algorithm, Proc.13th ACM-SIAM \n";
00172   cout << "# Symp. Discr. Alg, p.472 (2002), S.Fortune's Voronoi algorithm and code " ;
00173 #ifndef DROP_CGAL
00174   cout << endl << "# and CGAL: http://www.cgal.org/";
00175 #endif  // DROP_CGAL
00176   cout << ".\n";
00177   cout << "#-------------------------------------------------------------------------\n";
00178 }
00179 
00180 //----------------------------------------------------------------------
00181 // transfer all relevant info into internal variables
00182 void ClusterSequence::_decant_options(const JetDefinition & jet_def,
00183                                       const bool & writeout_combinations) {
00184 
00185   // let the user know what's going on
00186   _print_banner();
00187 
00188   // make a local copy of the jet definition (for future use?)
00189   _jet_def = jet_def;
00190   
00191   _writeout_combinations = writeout_combinations;
00192   _jet_algorithm = jet_def.jet_algorithm();
00193   _Rparam = jet_def.R();  _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
00194   _strategy = jet_def.strategy();
00195 
00196   // disallow interference from the plugin
00197   _plugin_activated = false;
00198   
00199 }
00200 
00201 
00202 //----------------------------------------------------------------------
00203 // initialise the history in a standard way
00204 void ClusterSequence::_fill_initial_history () {
00205 
00206   if (_jets.size() == 0) {throw Error("Cannot run jet-finder on empty event");}
00207 
00208   // reserve sufficient space for everything
00209   _jets.reserve(_jets.size()*2);
00210   _history.reserve(_jets.size()*2);
00211 
00212   for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
00213     history_element element;
00214     element.parent1 = InexistentParent;
00215     element.parent2 = InexistentParent;
00216     element.child   = Invalid;
00217     element.jetp_index = i;
00218     element.dij     = 0.0;
00219     element.max_dij_so_far = 0.0;
00220 
00221     _history.push_back(element);
00222     
00223     // do any momentum preprocessing needed by the recombination scheme
00224     _jet_def.recombiner()->preprocess(_jets[i]);
00225 
00226     // get cross-referencing right from PseudoJets
00227     _jets[i].set_cluster_hist_index(i);
00228   }
00229   _initial_n = _jets.size();
00230 }
00231 
00232 
00233 //----------------------------------------------------------------------
00234 // Return the component corresponding to the specified index.
00235 // taken from CLHEP
00236 string ClusterSequence::strategy_string ()  const {
00237   string strategy;
00238   switch(_strategy) {
00239   case NlnN:
00240     strategy = "NlnN"; break;
00241   case NlnN3pi:
00242     strategy = "NlnN3pi"; break;
00243   case NlnN4pi:
00244     strategy = "NlnN4pi"; break;
00245   case N2Plain:
00246     strategy = "N2Plain"; break;
00247   case N2Tiled:
00248     strategy = "N2Tiled"; break;
00249   case N2MinHeapTiled:
00250     strategy = "N2MinHeapTiled"; break;
00251   case N2PoorTiled:
00252     strategy = "N2PoorTiled"; break;
00253   case N3Dumb:
00254     strategy = "N3Dumb"; break;
00255   case NlnNCam4pi:
00256     strategy = "NlnNCam4pi"; break;
00257   case NlnNCam2pi2R:
00258     strategy = "NlnNCam2pi2R"; break;
00259   case NlnNCam:
00260     strategy = "NlnNCam"; break; // 2piMultD
00261   case plugin_strategy:
00262     strategy = "plugin strategy"; break;
00263   default:
00264     strategy = "Unrecognized";
00265   }
00266   return strategy;
00267 }  
00268 
00269 
00270 double ClusterSequence::jet_scale_for_algorithm(
00271                                   const PseudoJet & jet) const {
00272   if (_jet_algorithm == kt_algorithm)             {return jet.kt2();}
00273   else if (_jet_algorithm == cambridge_algorithm) {return 1.0;}
00274   else if (_jet_algorithm == antikt_algorithm) {
00275     double kt2=jet.kt2();
00276     return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
00277   } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
00278     double kt2 = jet.kt2();
00279     double lim = _jet_def.extra_param();
00280     if (kt2 < lim*lim && kt2 != 0.0) {
00281       return 1.0/kt2;
00282     } else {return 1.0;}
00283   } else {throw Error("Unrecognised jet finder");}
00284 }
00285 
00286 
00287 //----------------------------------------------------------------------
00291 void ClusterSequence::transfer_from_sequence(ClusterSequence & from_seq) {
00292 
00293   // the metadata
00294   _jet_def                 = from_seq._jet_def                ;
00295   _writeout_combinations   = from_seq._writeout_combinations  ;
00296   _initial_n               = from_seq._initial_n              ;
00297   _Rparam                  = from_seq._Rparam                 ;
00298   _R2                      = from_seq._R2                     ;
00299   _invR2                   = from_seq._invR2                  ;
00300   _strategy                = from_seq._strategy               ;
00301   _jet_algorithm           = from_seq._jet_algorithm          ;
00302   _plugin_activated        = from_seq._plugin_activated       ;
00303 
00304   // the data
00305   _jets     = from_seq._jets;
00306   _history  = from_seq._history;
00307   // the following transferse ownership of the extras from the from_seq
00308   _extras   = from_seq._extras;
00309 
00310 }
00311 
00312 //----------------------------------------------------------------------
00313 // record an ij recombination and reset the _jets[newjet_k] momentum and
00314 // user index to be those of newjet
00315 void ClusterSequence::plugin_record_ij_recombination(
00316            int jet_i, int jet_j, double dij, 
00317            const PseudoJet & newjet, int & newjet_k) {
00318 
00319   plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
00320 
00321   // now transfer newjet into place
00322   int tmp_index = _jets[newjet_k].cluster_hist_index();
00323   _jets[newjet_k] = newjet;
00324   _jets[newjet_k].set_cluster_hist_index(tmp_index);
00325 }
00326 
00327 
00328 //----------------------------------------------------------------------
00329 // return all inclusive jets with pt > ptmin
00330 vector<PseudoJet> ClusterSequence::inclusive_jets (const double & ptmin) const{
00331   double dcut = ptmin*ptmin;
00332   int i = _history.size() - 1; // last jet
00333   vector<PseudoJet> jets;
00334   if (_jet_algorithm == kt_algorithm) {
00335     while (i >= 0) {
00336       // with our specific definition of dij and diB (i.e. R appears only in 
00337       // dij), then dij==diB is the same as the jet.perp2() and we can exploit
00338       // this in selecting the jets...
00339       if (_history[i].max_dij_so_far < dcut) {break;}
00340       if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
00341         // for beam jets
00342         int parent1 = _history[i].parent1;
00343         jets.push_back(_jets[_history[parent1].jetp_index]);}
00344       i--;
00345     }
00346   } else if (_jet_algorithm == cambridge_algorithm) {
00347     while (i >= 0) {
00348       // inclusive jets are all at end of clustering sequence in the
00349       // Cambridge algorithm -- so if we find a non-exclusive jet, then
00350       // we can exit
00351       if (_history[i].parent2 != BeamJet) {break;}
00352       int parent1 = _history[i].parent1;
00353       const PseudoJet & jet = _jets[_history[parent1].jetp_index];
00354       if (jet.perp2() >= dcut) {jets.push_back(jet);}
00355       i--;
00356     }
00357   } else if (_jet_algorithm == plugin_algorithm 
00358              || _jet_algorithm == antikt_algorithm
00359              || _jet_algorithm == cambridge_for_passive_algorithm) {
00360     // for inclusive jets with a plugin algorithm, we make no
00361     // assumptions about anything (relation of dij to momenta,
00362     // ordering of the dij, etc.)
00363     while (i >= 0) {
00364       if (_history[i].parent2 == BeamJet) {
00365         int parent1 = _history[i].parent1;
00366         const PseudoJet & jet = _jets[_history[parent1].jetp_index];
00367         if (jet.perp2() >= dcut) {jets.push_back(jet);}
00368       }
00369       i--;
00370     }
00371   } else {throw Error("Unrecognized jet algorithm");}
00372   return jets;
00373 }
00374 
00375 
00376 //----------------------------------------------------------------------
00377 // return the number of exclusive jets that would have been obtained
00378 // running the algorithm in exclusive mode with the given dcut
00379 int ClusterSequence::n_exclusive_jets (const double & dcut) const {
00380 
00381   // first locate the point where clustering would have stopped (i.e. the
00382   // first time max_dij_so_far > dcut)
00383   int i = _history.size() - 1; // last jet
00384   while (i >= 0) {
00385     if (_history[i].max_dij_so_far <= dcut) {break;}
00386     i--;
00387   }
00388   int stop_point = i + 1;
00389   // relation between stop_point, njets assumes one extra jet disappears
00390   // at each clustering.
00391   int njets = 2*_initial_n - stop_point;
00392   return njets;
00393 }
00394 
00395 //----------------------------------------------------------------------
00396 // return all exclusive jets that would have been obtained running
00397 // the algorithm in exclusive mode with the given dcut
00398 vector<PseudoJet> ClusterSequence::exclusive_jets (const double & dcut) const {
00399   int njets = n_exclusive_jets(dcut);
00400   return exclusive_jets(njets);
00401 }
00402 
00403 
00404 //----------------------------------------------------------------------
00405 // return the jets obtained by clustering the event to n jets.
00406 vector<PseudoJet> ClusterSequence::exclusive_jets (const int & njets) const {
00407 
00408   // make sure the user does not ask for more than jets than there
00409   // were particles in the first place.
00410   assert (njets <= _initial_n);
00411 
00412   // provide a warning when extracting exclusive jets 
00413   if (_jet_def.jet_algorithm() != kt_algorithm && _n_exclusive_warnings < 5) {
00414     _n_exclusive_warnings++;
00415     cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
00416   }
00417 
00418 
00419   // calculate the point where we have to stop the clustering.
00420   // relation between stop_point, njets assumes one extra jet disappears
00421   // at each clustering.
00422   int stop_point = 2*_initial_n - njets;
00423 
00424   // some sanity checking to make sure that e+e- does not give us
00425   // surprises (should we ever implement e+e-)...
00426   if (2*_initial_n != static_cast<int>(_history.size())) {
00427     ostringstream err;
00428     err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
00429     throw Error(err.str());
00430     //assert(false);
00431   }
00432 
00433   // now go forwards and reconstitute the jets that we have --
00434   // basically for any history element, see if the parent jets to
00435   // which it refers were created before the stopping point -- if they
00436   // were then add them to the list, otherwise they are subsequent
00437   // recombinations of the jets that we are looking for.
00438   vector<PseudoJet> jets;
00439   for (unsigned int i = stop_point; i < _history.size(); i++) {
00440     int parent1 = _history[i].parent1;
00441     if (parent1 < stop_point) {
00442       jets.push_back(_jets[_history[parent1].jetp_index]);
00443     }
00444     int parent2 = _history[i].parent2;
00445     if (parent2 < stop_point && parent2 > 0) {
00446       jets.push_back(_jets[_history[parent2].jetp_index]);
00447     }
00448     
00449   }
00450 
00451   // sanity check...
00452   if (static_cast<int>(jets.size()) != njets) {
00453     ostringstream err;
00454     err << "ClusterSequence::exclusive_jets: size of returned vector ("
00455          <<jets.size()<<") does not coincide with requested number of jets ("
00456          <<njets<<")";
00457     throw Error(err.str());
00458   }
00459 
00460   return jets;
00461 }
00462 
00463 //----------------------------------------------------------------------
00466 double ClusterSequence::exclusive_dmerge (const int & njets) const {
00467   assert(njets > 0);
00468   if (njets >= _initial_n) {return 0.0;}
00469   return _history[2*_initial_n-njets-1].dij;
00470 }
00471 
00472 
00473 //----------------------------------------------------------------------
00478 double ClusterSequence::exclusive_dmerge_max (const int & njets) const {
00479   assert(njets > 0);
00480   if (njets >= _initial_n) {return 0.0;}
00481   return _history[2*_initial_n-njets-1].max_dij_so_far;
00482 }
00483 
00484 
00485 //----------------------------------------------------------------------
00486 // work through the object's history until
00487 bool ClusterSequence::object_in_jet(const PseudoJet & object, 
00488                                     const PseudoJet & jet) const {
00489 
00490   // make sure the object conceivably belongs to this clustering
00491   // sequence
00492   assert(_potentially_valid(object) && _potentially_valid(jet));
00493 
00494   const PseudoJet * this_object = &object;
00495   const PseudoJet * childp;
00496   while(true) {
00497     if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
00498       return true;
00499     } else if (has_child(*this_object, childp)) {this_object = childp;}
00500     else {return false;}
00501   }
00502 }
00503 
00504 //----------------------------------------------------------------------
00510 bool ClusterSequence::has_parents(const PseudoJet & jet, PseudoJet & parent1, 
00511                               PseudoJet & parent2) const {
00512 
00513   const history_element & hist = _history[jet.cluster_hist_index()];
00514 
00515   // make sure we do not run into any unexpected situations --
00516   // i.e. both parents valid, or neither
00517   assert ((hist.parent1 >= 0 && hist.parent2 >= 0) || 
00518           (hist.parent1 < 0 && hist.parent2 < 0));
00519 
00520   if (hist.parent1 < 0) {
00521     parent1 = PseudoJet(0.0,0.0,0.0,0.0);
00522     parent2 = parent1;
00523     return false;
00524   } else {
00525     parent1 = _jets[_history[hist.parent1].jetp_index];
00526     parent2 = _jets[_history[hist.parent2].jetp_index];
00527     // order the parents in decreasing pt
00528     if (parent1.perp2() < parent2.perp2()) swap(parent1,parent2);
00529     return true;
00530   }
00531 }
00532 
00533 //----------------------------------------------------------------------
00536 bool ClusterSequence::has_child(const PseudoJet & jet, PseudoJet & child) const {
00537 
00538   //const history_element & hist = _history[jet.cluster_hist_index()];
00539   //
00540   //if (hist.child >= 0) {
00541   //  child = _jets[_history[hist.child].jetp_index];
00542   //  return true;
00543   //} else {
00544   //  child = PseudoJet(0.0,0.0,0.0,0.0);
00545   //  return false;
00546   //}
00547   const PseudoJet * childp;
00548   bool res = has_child(jet, childp);
00549   if (res) {
00550     child = *childp;
00551     return true;
00552   } else {
00553     child = PseudoJet(0.0,0.0,0.0,0.0);
00554     return false;
00555   }
00556 }
00557 
00558 bool ClusterSequence::has_child(const PseudoJet & jet, const PseudoJet * & childp) const {
00559 
00560   const history_element & hist = _history[jet.cluster_hist_index()];
00561 
00562   // check that this jet has a child and that the child corresponds to
00563   // a true jet [RETHINK-IF-CHANGE-NUMBERING: what is the right
00564   // behaviour if the child is the same jet but made inclusive...?]
00565   if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
00566     childp = &(_jets[_history[hist.child].jetp_index]);
00567     return true;
00568   } else {
00569     childp = NULL;
00570     return false;
00571   }
00572 }
00573 
00574 
00575 //----------------------------------------------------------------------
00579 bool ClusterSequence::has_partner(const PseudoJet & jet, 
00580                               PseudoJet & partner) const {
00581 
00582   const history_element & hist = _history[jet.cluster_hist_index()];
00583 
00584   // make sure we have a child and that the child does not correspond
00585   // to a clustering with the beam (or some other invalid quantity)
00586   if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
00587     const history_element & child_hist = _history[hist.child];
00588     if (child_hist.parent1 == jet.cluster_hist_index()) {
00589       // partner will be child's parent2 -- for iB clustering
00590       // parent2 will not be valid
00591       partner = _jets[_history[child_hist.parent2].jetp_index];
00592     } else {
00593       // partner will be child's parent1
00594       partner = _jets[_history[child_hist.parent1].jetp_index];
00595     }
00596     return true;
00597   } else {
00598     partner = PseudoJet(0.0,0.0,0.0,0.0);
00599     return false;
00600   }
00601 }
00602 
00603 
00604 //----------------------------------------------------------------------
00605 // return a vector of the particles that make up a jet
00606 vector<PseudoJet> ClusterSequence::constituents (const PseudoJet & jet) const {
00607   vector<PseudoJet> subjets;
00608   add_constituents(jet, subjets);
00609   return subjets;
00610 }
00611 
00612 //----------------------------------------------------------------------
00621 void ClusterSequence::print_jets_for_root(const std::vector<PseudoJet> & jets, 
00622                                           ostream & ostr) const {
00623   for (unsigned i = 0; i < jets.size(); i++) {
00624     ostr << i  << " "
00625          << jets[i].px() << " "
00626          << jets[i].py() << " "
00627          << jets[i].pz() << " "
00628          << jets[i].E() << endl;
00629     vector<PseudoJet> cst = constituents(jets[i]);
00630     for (unsigned j = 0; j < cst.size() ; j++) {
00631       ostr << " " << j << " "
00632            << cst[j].rap() << " "
00633            << cst[j].phi() << " "
00634            << cst[j].perp() << endl;
00635     }
00636     ostr << "#END" << endl;
00637   }
00638 }
00639 
00640 // Not yet. Perhaps in a future release
00641 // //----------------------------------------------------------------------
00642 // // print out all inclusive jets with pt > ptmin
00643 // void ClusterSequence::print_jets (const double & ptmin) const{
00644 //     vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
00645 // 
00646 //     for (size_t j = 0; j < jets.size(); j++) {
00647 //        printf("%5u %7.3f %7.3f %9.3f\n",
00648 //        j,jets[j].rap(),jets[j].phi(),jets[j].perp());
00649 //     }
00650 // }
00651 
00652 //----------------------------------------------------------------------
00657 vector<int> ClusterSequence::particle_jet_indices(
00658                         const vector<PseudoJet> & jets) const {
00659 
00660   vector<int> indices(n_particles());
00661 
00662   // first label all particles as not belonging to any jets
00663   for (unsigned ipart = 0; ipart < n_particles(); ipart++) 
00664     indices[ipart] = -1;
00665 
00666   // then for each of the jets relabel its consituents as belonging to
00667   // that jet
00668   for (unsigned ijet = 0; ijet < jets.size(); ijet++) {
00669 
00670     vector<PseudoJet> jet_constituents(constituents(jets[ijet]));
00671 
00672     for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
00673       // a safe (if slightly redundant) way of getting the particle
00674       // index (for initial particles it is actually safe to assume
00675       // ipart=iclust).
00676       unsigned iclust = jet_constituents[ip].cluster_hist_index();
00677       unsigned ipart = history()[iclust].jetp_index;
00678       indices[ipart] = ijet;
00679     }
00680   }
00681 
00682   return indices;
00683 }
00684 
00685 
00686 //----------------------------------------------------------------------
00687 // recursive routine that adds on constituents of jet to the subjet_vector
00688 void ClusterSequence::add_constituents (
00689            const PseudoJet & jet, vector<PseudoJet> & subjet_vector) const {
00690   // find out position in cluster history
00691   int i = jet.cluster_hist_index();
00692   int parent1 = _history[i].parent1;
00693   int parent2 = _history[i].parent2;
00694 
00695   if (parent1 == InexistentParent) {
00696     // It is an original particle (labelled by its parent having value
00697     // InexistentParent), therefore add it on to the subjet vector
00698     subjet_vector.push_back(jet);
00699     return;
00700   } 
00701 
00702   // add parent 1
00703   add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
00704 
00705   // see if parent2 is a real jet; if it is then add its constituents
00706   if (parent2 != BeamJet) {
00707     add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
00708   }
00709 }
00710 
00711 
00712 
00713 //----------------------------------------------------------------------
00714 // initialise the history in a standard way
00715 void ClusterSequence::_add_step_to_history (
00716                const int & step_number, const int & parent1, 
00717                const int & parent2, const int & jetp_index,
00718                const double & dij) {
00719 
00720   history_element element;
00721   element.parent1 = parent1;
00722   element.parent2 = parent2;
00723   element.jetp_index = jetp_index;
00724   element.child = Invalid;
00725   element.dij   = dij;
00726   element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
00727   _history.push_back(element);
00728 
00729   int local_step = _history.size()-1;
00730   assert(local_step == step_number);
00731 
00732   assert(parent1 >= 0);
00733   _history[parent1].child = local_step;
00734   if (parent2 >= 0) {_history[parent2].child = local_step;}
00735 
00736   // get cross-referencing right from PseudoJets
00737   if (jetp_index != Invalid) {
00738     assert(jetp_index >= 0);
00739     //cout << _jets.size() <<" "<<jetp_index<<"\n";
00740     _jets[jetp_index].set_cluster_hist_index(local_step);
00741   }
00742 
00743   if (_writeout_combinations) {
00744     cout << local_step << ": " 
00745          << parent1 << " with " << parent2
00746          << "; y = "<< dij<<endl;
00747   }
00748 
00749 }
00750 
00751 
00752 
00753 
00754 //======================================================================
00755 // Return an order in which to read the history such that _history[order[i]] 
00756 // will always correspond to the same set of consituent particles if 
00757 // two branching histories are equivalent in terms of the particles
00758 // contained in any given pseudojet.
00759 vector<int> ClusterSequence::unique_history_order() const {
00760 
00761   // first construct an array that will tell us the lowest constituent
00762   // of a given jet -- this will always be one of the original
00763   // particles, whose order is well defined and so will help us to
00764   // follow the tree in a unique manner.
00765   valarray<int> lowest_constituent(_history.size());
00766   int hist_n = _history.size();
00767   lowest_constituent = hist_n; // give it a large number
00768   for (int i = 0; i < hist_n; i++) {
00769     // sets things up for the initial partons
00770     lowest_constituent[i] = min(lowest_constituent[i],i); 
00771     // propagates them through to the children of this parton
00772     if (_history[i].child > 0) lowest_constituent[_history[i].child] 
00773       = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
00774   }
00775 
00776   // establish an array for what we have and have not extracted so far
00777   valarray<bool> extracted(_history.size()); extracted = false;
00778   vector<int> unique_tree;
00779   unique_tree.reserve(_history.size());
00780 
00781   // now work our way through the tree
00782   for (unsigned i = 0; i < n_particles(); i++) {
00783     if (!extracted[i]) {
00784       unique_tree.push_back(i);
00785       extracted[i] = true;
00786       _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
00787     }
00788   }
00789 
00790   return unique_tree;
00791 }
00792 
00793 //======================================================================
00794 // helper for unique_history_order
00795 void ClusterSequence::_extract_tree_children(
00796        int position, 
00797        valarray<bool> & extracted, 
00798        const valarray<int> & lowest_constituent,
00799        vector<int> & unique_tree) const {
00800   if (!extracted[position]) {
00801     // that means we may have unidentified parents around, so go and
00802     // collect them (extracted[position]) will then be made true)
00803     _extract_tree_parents(position,extracted,lowest_constituent,unique_tree);
00804   } 
00805   
00806   // now look after the children...
00807   int child = _history[position].child;
00808   if (child  >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree);
00809 }
00810 
00811 
00812 //======================================================================
00813 // return the list of unclustered particles
00814 vector<PseudoJet> ClusterSequence::unclustered_particles() const {
00815   vector<PseudoJet> unclustered;
00816   for (unsigned i = 0; i < n_particles() ; i++) {
00817     if (_history[i].child == Invalid) 
00818       unclustered.push_back(_jets[_history[i].jetp_index]);
00819   }
00820   return unclustered;
00821 }
00822 
00823 
00824 
00825 //======================================================================
00826 // helper for unique_history_order
00827 void ClusterSequence::_extract_tree_parents(
00828        int position, 
00829        valarray<bool> & extracted, 
00830        const valarray<int> & lowest_constituent,
00831        vector<int> & unique_tree) const {
00832 
00833   if (!extracted[position]) {
00834     int parent1 = _history[position].parent1;
00835     int parent2 = _history[position].parent2;
00836     // where relevant order parents so that we will first treat the
00837     // one containing the smaller "lowest_constituent"
00838     if (parent1 >= 0 && parent2 >= 0) {
00839       if (lowest_constituent[parent1] > lowest_constituent[parent2]) 
00840         swap(parent1, parent2);
00841     }
00842     // then actually run through the parents to extract the constituents...
00843     if (parent1 >= 0 && !extracted[parent1]) 
00844       _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree);
00845     if (parent2 >= 0 && !extracted[parent2]) 
00846       _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree);
00847     // finally declare this position to be accounted for and push it
00848     // onto our list.
00849     unique_tree.push_back(position);
00850     extracted[position] = true;
00851   }
00852 }
00853 
00854 
00855 //======================================================================
00859 void ClusterSequence::_do_ij_recombination_step(
00860                                const int & jet_i, const int & jet_j, 
00861                                const double & dij, 
00862                                int & newjet_k) {
00863 
00864   // create the new jet by recombining the first two
00865   PseudoJet newjet;
00866   _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
00867   _jets.push_back(newjet);
00868   // original version...
00869   //_jets.push_back(_jets[jet_i] + _jets[jet_j]);
00870 
00871   // get its index
00872   newjet_k = _jets.size()-1;
00873 
00874   // get history index
00875   int newstep_k = _history.size();
00876   // and provide jet with the info
00877   _jets[newjet_k].set_cluster_hist_index(newstep_k);
00878 
00879   // finally sort out the history 
00880   int hist_i = _jets[jet_i].cluster_hist_index();
00881   int hist_j = _jets[jet_j].cluster_hist_index();
00882 
00883   _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
00884                        newjet_k, dij);
00885 
00886 }
00887 
00888 
00889 //======================================================================
00892 void ClusterSequence::_do_iB_recombination_step(
00893                                   const int & jet_i, const double & diB) {
00894   // get history index
00895   int newstep_k = _history.size();
00896 
00897   // recombine the jet with the beam
00898   _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
00899                        Invalid, diB);
00900 
00901 }
00902 
00903 FASTJET_END_NAMESPACE
00904 

Generated on Tue Dec 18 17:05:03 2007 for fastjet by  doxygen 1.5.2