ClusterSequence.hh

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: ClusterSequence.hh 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 
00032 //----------------------------------------------------------------------
00033 // here's where we put the main page for fastjet (as explained in the
00034 // Doxygen faq)
00035 //......................................................................
00044 //----------------------------------------------------------------------
00045 
00046 #ifndef __FASTJET_CLUSTERSEQUENCE_HH__
00047 #define __FASTJET_CLUSTERSEQUENCE_HH__
00048 
00049 #include<vector>
00050 #include<map>
00051 #include "fastjet/internal/DynamicNearestNeighbours.hh"
00052 #include "fastjet/PseudoJet.hh"
00053 #include<memory>
00054 #include<cassert>
00055 #include<iostream>
00056 #include<string>
00057 #include<cmath> // needed to get double std::abs(double)
00058 #include "fastjet/Error.hh"
00059 #include "fastjet/JetDefinition.hh"
00060 
00061 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00062 
00063 
00065 class ClusterSequence {
00066 
00067 
00068  public: 
00069 
00071   ClusterSequence () {}
00072 
00082   template<class L> ClusterSequence (const std::vector<L> & pseudojets, 
00083                    const double & R = 1.0,
00084                    const Strategy & strategy = Best,
00085                    const bool & writeout_combinations = false);
00086 
00087 
00091   template<class L> ClusterSequence (
00092                                   const std::vector<L> & pseudojets,
00093                                   const JetDefinition & jet_def,
00094                                   const bool & writeout_combinations = false);
00095 
00096 
00097   // NB: in the routines that follow, for extracting lists of jets, a
00098   //     list structure might be more efficient, if sometimes a little
00099   //     more awkward to use (at least for old fortran hands).
00100 
00104   std::vector<PseudoJet> inclusive_jets (const double & ptmin = 0.0) const;
00105 
00109   int n_exclusive_jets (const double & dcut) const;
00110 
00114   std::vector<PseudoJet> exclusive_jets (const double & dcut) const;
00115 
00118   std::vector<PseudoJet> exclusive_jets (const int & njets) const;
00119 
00122   double exclusive_dmerge (const int & njets) const;
00123 
00128   double exclusive_dmerge_max (const int & njets) const;
00129 
00140   bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;
00141 
00147   bool has_parents(const PseudoJet & jet, PseudoJet & parent1, 
00148                PseudoJet & parent2) const;
00149 
00152   bool has_child(const PseudoJet & jet, PseudoJet & child) const;
00153 
00156   bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;
00157 
00161   bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;
00162 
00164   std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
00165 
00166 
00175   void print_jets_for_root(const std::vector<PseudoJet> & jets, 
00176                            std::ostream & ostr = std::cout) const;
00177 
00178 // Not yet. Perhaps in a future release.
00179 //   /// print out all inclusive jets with pt > ptmin
00180 //   virtual void print_jets (const double & ptmin=0.0) const;
00181 
00183   void add_constituents (const PseudoJet & jet, 
00184                          std::vector<PseudoJet> & subjet_vector) const;
00185 
00187   inline Strategy strategy_used () const {return _strategy;}
00188   std::string strategy_string () const;
00189 
00191   const JetDefinition & jet_def() const {return _jet_def;}
00192 
00196   double jet_scale_for_algorithm(const PseudoJet & jet) const;
00197 
00198   //----- next follow functions designed specifically for plugins, which
00199   //      may only be called when plugin_activated() returns true
00200 
00206   void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, 
00207                                       int & newjet_k) {
00208     assert(plugin_activated());
00209     _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
00210   }
00211 
00215   void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, 
00216                                       const PseudoJet & newjet, 
00217                                       int & newjet_k);
00218 
00223   void plugin_record_iB_recombination(int jet_i, double diB) {
00224     assert(plugin_activated());
00225     _do_iB_recombination_step(jet_i, diB);
00226   }
00227 
00231   class Extras {
00232   public:
00233     virtual ~Extras() {}
00234     virtual std::string description() const {return "This is a dummy extras class that contains no extra information! Derive from it if you want to use it to provide extra information from a plugin jet finder";}
00235   };
00236 
00239   inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) {
00240     _extras = extras_in;
00241   }
00242 
00244   inline bool plugin_activated() const {return _plugin_activated;}
00245 
00247   const Extras * extras() const {return _extras.get();}
00248 
00249 public:
00253   static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
00255   static void set_jet_finder (JetAlgorithm jet_algorithm)    {_default_jet_algorithm = jet_algorithm;}
00256 
00257 
00260   struct history_element{
00261     int parent1; 
00262 
00263 
00264 
00265     int parent2; 
00266 
00267 
00268 
00269 
00270 
00271     int child;   
00272 
00273 
00274 
00275 
00276     int jetp_index; 
00277 
00278 
00279 
00280 
00281 
00282 
00283     double dij;  
00284 
00285 
00286     double max_dij_so_far; 
00287 
00288   };
00289 
00290   enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
00291 
00297   const std::vector<PseudoJet> & jets()    const;
00298 
00301   const std::vector<history_element> & history() const;
00302 
00307   unsigned int n_particles() const;
00308 
00313   std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
00314 
00330   std::vector<int> unique_history_order() const;
00331 
00335   std::vector<PseudoJet> unclustered_particles() const;
00336 
00340   void transfer_from_sequence(ClusterSequence & from_seq);
00341 
00342 protected:
00343   static JetAlgorithm _default_jet_algorithm;
00344   JetDefinition _jet_def;
00345 
00348   bool _potentially_valid(const PseudoJet & jet) const {
00349     return jet.cluster_hist_index() >= 0 
00350       && jet.cluster_hist_index() < int(_history.size());
00351   }
00352 
00355   template<class L> void _transfer_input_jets(
00356                                      const std::vector<L> & pseudojets);
00357 
00361   void _initialise_and_run (const JetDefinition & jet_def,
00362                             const bool & writeout_combinations);
00363 
00367   void _initialise_and_run (const double & R,
00368                             const Strategy & strategy,
00369                             const bool & writeout_combinations);
00370 
00373   void _decant_options(const JetDefinition & jet_def,
00374                        const bool & writeout_combinations);
00375 
00379   void _fill_initial_history();
00380 
00384   void _do_ij_recombination_step(const int & jet_i, const int & jet_j, 
00385                                  const double & dij, int & newjet_k);
00386 
00389   void _do_iB_recombination_step(const int & jet_i, const double & diB);
00390 
00391 
00395   std::vector<PseudoJet> _jets;
00396 
00397 
00401   std::vector<history_element> _history;
00402 
00403   bool _writeout_combinations;
00404   int  _initial_n;
00405   double _Rparam, _R2, _invR2;
00406   Strategy    _strategy;
00407   JetAlgorithm  _jet_algorithm;
00408 
00409  private:
00410 
00411   bool _plugin_activated;
00412   std::auto_ptr<Extras> _extras; // things the plugin might want to add
00413 
00414   void _really_dumb_cluster ();
00415   void _delaunay_cluster ();
00416   void _simple_N2_cluster ();
00417   void _tiled_N2_cluster ();
00418   void _faster_tiled_N2_cluster ();
00419 
00420   //
00421   void _minheap_faster_tiled_N2_cluster();
00422 
00423   // things needed specifically for Cambridge with Chan's 2D closest
00424   // pairs method
00425   void _CP2DChan_cluster();
00426   void _CP2DChan_cluster_2pi2R ();
00427   void _CP2DChan_cluster_2piMultD ();
00428   void _CP2DChan_limited_cluster(double D);
00429   void _do_Cambridge_inclusive_jets();
00430 
00431   void _add_step_to_history(const int & step_number, const int & parent1, 
00432                                const int & parent2, const int & jetp_index,
00433                                const double & dij);
00434 
00437   void _extract_tree_children(int pos, std::valarray<bool> &, 
00438                 const std::valarray<int> &, std::vector<int> &) const;
00439 
00442   void _extract_tree_parents (int pos, std::valarray<bool> &, 
00443                 const std::valarray<int> &,  std::vector<int> &) const;
00444 
00445 
00446   // these will be useful shorthands in the Voronoi-based code
00447   typedef std::pair<int,int> TwoVertices;
00448   typedef std::pair<double,TwoVertices> DijEntry;
00449   typedef std::multimap<double,TwoVertices> DistMap;
00450 
00452   void _add_ktdistance_to_map(const int & ii, 
00453                               DistMap & DijMap,
00454                               const DynamicNearestNeighbours * DNN);
00455 
00457   void _print_banner();
00459   static bool _first_time;
00460 
00464   static int _n_exclusive_warnings;
00465 
00466   //----------------------------------------------------------------------
00471   struct BriefJet {
00472     double     eta, phi, kt2, NN_dist;
00473     BriefJet * NN;
00474     int        _jets_index;
00475   };
00478   class TiledJet {
00479   public:
00480     double     eta, phi, kt2, NN_dist;
00481     TiledJet * NN, *previous, * next; 
00482     int        _jets_index, tile_index, diJ_posn;
00483     // routines that are useful in the minheap version of tiled
00484     // clustering ("misuse" the otherwise unused diJ_posn, so as
00485     // to indicate whether jets need to have their minheap entries
00486     // updated).
00487     inline void label_minheap_update_needed() {diJ_posn = 1;}
00488     inline void label_minheap_update_done()   {diJ_posn = 0;}
00489     inline bool minheap_update_needed() const {return diJ_posn==1;}
00490   };
00491 
00492   //-- some of the functions that follow are templates and will work
00493   //as well for briefjet and tiled jets
00494 
00497   template <class J> void _bj_set_jetinfo( J * const jet, 
00498                                                  const int _jets_index) const;
00499 
00502   void _bj_remove_from_tiles( TiledJet * const jet) const;
00503 
00505   template <class J> double _bj_dist(const J * const jeta, 
00506                         const J * const jetb) const;
00507 
00508   // return the diJ (multiplied by _R2) for this jet assuming its NN
00509   // info is correct
00510   template <class J> double _bj_diJ(const J * const jeta) const;
00511 
00515   template <class J> inline J * _bj_of_hindex(
00516                           const int hist_index, 
00517                           J * const head, J * const tail) 
00518     const {
00519     J * res;
00520     for(res = head; res<tail; res++) {
00521       if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
00522     }
00523     return res;
00524   }
00525 
00526 
00527   //-- remaining functions are different in various cases, so we
00528   //   will use templates but are not sure if they're useful...
00529 
00534   template <class J> void _bj_set_NN_nocross(J * const jeta, 
00535             J * const head, const J * const tail) const;
00536 
00541   template <class J> void _bj_set_NN_crosscheck(J * const jeta, 
00542             J * const head, const J * const tail) const;
00543   
00544 
00545 
00548   static const int n_tile_neighbours = 9;
00549   //----------------------------------------------------------------------
00552   struct Tile {
00554     Tile *   begin_tiles[n_tile_neighbours]; 
00556     Tile **  surrounding_tiles; 
00558     Tile **  RH_tiles;  
00560     Tile **  end_tiles; 
00562     TiledJet * head;    
00564     bool     tagged;    
00565   };
00566   std::vector<Tile> _tiles;
00567   double _tiles_eta_min, _tiles_eta_max;
00568   double _tile_size_eta, _tile_size_phi;
00569   int    _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
00570 
00571   // reasonably robust return of tile index given ieta and iphi, in particular
00572   // it works even if iphi is negative
00573   inline int _tile_index (int ieta, int iphi) const {
00574     // note that (-1)%n = -1 so that we have to add _n_tiles_phi
00575     // before performing modulo operation
00576     return (ieta-_tiles_ieta_min)*_n_tiles_phi
00577                   + (iphi+_n_tiles_phi) % _n_tiles_phi;
00578   }
00579 
00580   // routines for tiled case, including some overloads of the plain
00581   // BriefJet cases
00582   int  _tile_index(const double & eta, const double & phi) const;
00583   void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
00584   void  _bj_remove_from_tiles(TiledJet * const jet);
00585   void _initialise_tiles();
00586   void _print_tiles(TiledJet * briefjets ) const;
00587   void _add_neighbours_to_tile_union(const int tile_index, 
00588                  std::vector<int> & tile_union, int & n_near_tiles) const;
00589   void _add_untagged_neighbours_to_tile_union(const int tile_index, 
00590                  std::vector<int> & tile_union, int & n_near_tiles);
00591 
00592 
00593 };
00594 
00595 
00596 
00597 //**********************************************************************
00598 //**************    START   OF   INLINE   MATERIAL    ******************
00599 //**********************************************************************
00600 
00601 
00602 //----------------------------------------------------------------------
00603 // Transfer the initial jets into our internal structure
00604 template<class L> void ClusterSequence::_transfer_input_jets(
00605                                        const std::vector<L> & pseudojets) {
00606 
00607   // this will ensure that we can point to jets without difficulties
00608   // arising.
00609   _jets.reserve(pseudojets.size()*2);
00610 
00611   // insert initial jets this way so that any type L that can be
00612   // converted to a pseudojet will work fine (basically PseudoJet
00613   // and any type that has [] subscript access to the momentum
00614   // components, such as CLHEP HepLorentzVector).
00615   for (unsigned int i = 0; i < pseudojets.size(); i++) {
00616     _jets.push_back(pseudojets[i]);}
00617   
00618 }
00619 
00620 //----------------------------------------------------------------------
00621 // initialise from some generic type... Has to be made available
00622 // here in order for it the template aspect of it to work...
00623 template<class L> ClusterSequence::ClusterSequence (
00624                                   const std::vector<L> & pseudojets,
00625                                   const double & R,
00626                                   const Strategy & strategy,
00627                                   const bool & writeout_combinations) {
00628 
00629   // transfer the initial jets (type L) into our own array
00630   _transfer_input_jets(pseudojets);
00631 
00632   // run the clustering
00633   _initialise_and_run(R,strategy,writeout_combinations);
00634 }
00635 
00636 
00637 //----------------------------------------------------------------------
00640 template<class L> ClusterSequence::ClusterSequence (
00641                                   const std::vector<L> & pseudojets,
00642                                   const JetDefinition & jet_def,
00643                                   const bool & writeout_combinations) {
00644 
00645   // transfer the initial jets (type L) into our own array
00646   _transfer_input_jets(pseudojets);
00647 
00648   // run the clustering
00649   _initialise_and_run(jet_def,writeout_combinations);
00650 }
00651 
00652 
00653 inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
00654   return _jets;
00655 }
00656 
00657 inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
00658   return _history;
00659 }
00660 
00661 inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
00662 
00663 
00664 
00665 //----------------------------------------------------------------------
00666 template <class J> inline void ClusterSequence::_bj_set_jetinfo(
00667                             J * const jetA, const int _jets_index) const {
00668     jetA->eta  = _jets[_jets_index].rap();
00669     jetA->phi  = _jets[_jets_index].phi_02pi();
00670     jetA->kt2  = jet_scale_for_algorithm(_jets[_jets_index]);
00671     jetA->_jets_index = _jets_index;
00672     // initialise NN info as well
00673     jetA->NN_dist = _R2;
00674     jetA->NN      = NULL;
00675 }
00676 
00677 
00678 
00679 
00680 //----------------------------------------------------------------------
00681 template <class J> inline double ClusterSequence::_bj_dist(
00682                 const J * const jetA, const J * const jetB) const {
00683   double dphi = std::abs(jetA->phi - jetB->phi);
00684   double deta = (jetA->eta - jetB->eta);
00685   if (dphi > pi) {dphi = twopi - dphi;}
00686   return dphi*dphi + deta*deta;
00687 }
00688 
00689 //----------------------------------------------------------------------
00690 template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
00691   double kt2 = jet->kt2;
00692   if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
00693   return jet->NN_dist * kt2;
00694 }
00695 
00696 
00697 //----------------------------------------------------------------------
00698 // set the NN for jet without checking whether in the process you might
00699 // have discovered a new nearest neighbour for another jet
00700 template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
00701                  J * const jet, J * const head, const J * const tail) const {
00702   double NN_dist = _R2;
00703   J * NN  = NULL;
00704   if (head < jet) {
00705     for (J * jetB = head; jetB != jet; jetB++) {
00706       double dist = _bj_dist(jet,jetB);
00707       if (dist < NN_dist) {
00708         NN_dist = dist;
00709         NN = jetB;
00710       }
00711     }
00712   }
00713   if (tail > jet) {
00714     for (J * jetB = jet+1; jetB != tail; jetB++) {
00715       double dist = _bj_dist(jet,jetB);
00716       if (dist < NN_dist) {
00717         NN_dist = dist;
00718         NN = jetB;
00719       }
00720     }
00721   }
00722   jet->NN = NN;
00723   jet->NN_dist = NN_dist;
00724 }
00725 
00726 
00727 //----------------------------------------------------------------------
00728 template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet, 
00729                     J * const head, const J * const tail) const {
00730   double NN_dist = _R2;
00731   J * NN  = NULL;
00732   for (J * jetB = head; jetB != tail; jetB++) {
00733     double dist = _bj_dist(jet,jetB);
00734     if (dist < NN_dist) {
00735       NN_dist = dist;
00736       NN = jetB;
00737     }
00738     if (dist < jetB->NN_dist) {
00739       jetB->NN_dist = dist;
00740       jetB->NN = jet;
00741     }
00742   }
00743   jet->NN = NN;
00744   jet->NN_dist = NN_dist;
00745 }
00746 
00747 
00748 
00749 
00750 FASTJET_END_NAMESPACE
00751 
00752 #endif // __FASTJET_CLUSTERSEQUENCE_HH__

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