FastJet 3.0.2
ClusterSequence.hh
00001 //STARTHEADER
00002 // $Id: ClusterSequence.hh 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 
00030 #ifndef __FASTJET_CLUSTERSEQUENCE_HH__
00031 #define __FASTJET_CLUSTERSEQUENCE_HH__
00032 
00033 #include<vector>
00034 #include<map>
00035 #include "fastjet/internal/DynamicNearestNeighbours.hh"
00036 #include "fastjet/PseudoJet.hh"
00037 #include<memory>
00038 #include<cassert>
00039 #include<iostream>
00040 #include<string>
00041 #include<set>
00042 #include<cmath> // needed to get double std::abs(double)
00043 #include "fastjet/Error.hh"
00044 #include "fastjet/JetDefinition.hh"
00045 #include "fastjet/SharedPtr.hh"
00046 #include "fastjet/LimitedWarning.hh"
00047 #include "fastjet/FunctionOfPseudoJet.hh"
00048 #include "fastjet/ClusterSequenceStructure.hh"
00049 
00050 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00051 
00052 
00053 // forward declaration
00054 class ClusterSequenceStructure;
00055 
00056 /// @ingroup basic_classes
00057 /// \class ClusterSequence
00058 /// deals with clustering
00059 class ClusterSequence {
00060 
00061 
00062  public: 
00063 
00064   /// default constructor
00065   ClusterSequence () : _deletes_self_when_unused(false) {}
00066 
00067 //   /// create a clustersequence starting from the supplied set
00068 //   /// of pseudojets and clustering them with the long-invariant
00069 //   /// kt algorithm (E-scheme recombination) with the supplied
00070 //   /// value for R.
00071 //   ///
00072 //   /// If strategy=DumbN3 a very stupid N^3 algorithm is used for the
00073 //   /// clustering; otherwise strategy = NlnN* uses cylinders algorithms
00074 //   /// with some number of pi coverage. If writeout_combinations=true a
00075 //   /// summary of the recombination sequence is written out
00076 //   template<class L> ClusterSequence (const std::vector<L> & pseudojets, 
00077 //                 const double & R = 1.0,
00078 //                 const Strategy & strategy = Best,
00079 //                 const bool & writeout_combinations = false);
00080 
00081 
00082   /// create a clustersequence starting from the supplied set
00083   /// of pseudojets and clustering them with jet definition specified
00084   /// by jet_def (which also specifies the clustering strategy)
00085   template<class L> ClusterSequence (
00086                                   const std::vector<L> & pseudojets,
00087                                   const JetDefinition & jet_def,
00088                                   const bool & writeout_combinations = false);
00089   
00090   /// copy constructor for a ClusterSequence
00091   ClusterSequence (const ClusterSequence & cs) : _deletes_self_when_unused(false) {
00092     transfer_from_sequence(cs);
00093   }
00094 
00095   // virtual ClusterSequence destructor, in case any derived class
00096   // thinks of needing a destructor at some point
00097   virtual ~ClusterSequence (); //{}
00098 
00099   // NB: in the routines that follow, for extracting lists of jets, a
00100   //     list structure might be more efficient, if sometimes a little
00101   //     more awkward to use (at least for old fortran hands).
00102 
00103   /// return a vector of all jets (in the sense of the inclusive
00104   /// algorithm) with pt >= ptmin. Time taken should be of the order
00105   /// of the number of jets returned.
00106   std::vector<PseudoJet> inclusive_jets (const double & ptmin = 0.0) const;
00107 
00108   /// return the number of jets (in the sense of the exclusive
00109   /// algorithm) that would be obtained when running the algorithm
00110   /// with the given dcut.
00111   int n_exclusive_jets (const double & dcut) const;
00112 
00113   /// return a vector of all jets (in the sense of the exclusive
00114   /// algorithm) that would be obtained when running the algorithm
00115   /// with the given dcut.
00116   std::vector<PseudoJet> exclusive_jets (const double & dcut) const;
00117 
00118   /// return a vector of all jets when the event is clustered (in the
00119   /// exclusive sense) to exactly njets. 
00120   ///
00121   /// If there are fewer than njets particles in the ClusterSequence
00122   /// an error is thrown
00123   std::vector<PseudoJet> exclusive_jets (const int & njets) const;
00124 
00125   /// return a vector of all jets when the event is clustered (in the
00126   /// exclusive sense) to exactly njets. 
00127   ///
00128   /// If there are fewer than njets particles in the ClusterSequence
00129   /// the function just returns however many particles there were.
00130   std::vector<PseudoJet> exclusive_jets_up_to (const int & njets) const;
00131 
00132   /// return the dmin corresponding to the recombination that went
00133   /// from n+1 to n jets (sometimes known as d_{n n+1}). If the number
00134   /// of particles in the event is <= njets, the function returns 0.
00135   double exclusive_dmerge (const int & njets) const;
00136 
00137   /// return the maximum of the dmin encountered during all recombinations 
00138   /// up to the one that led to an n-jet final state; identical to
00139   /// exclusive_dmerge, except in cases where the dmin do not increase
00140   /// monotonically.
00141   double exclusive_dmerge_max (const int & njets) const;
00142 
00143   /// return the ymin corresponding to the recombination that went from
00144   /// n+1 to n jets (sometimes known as y_{n n+1}).
00145   double exclusive_ymerge (int njets) const {return exclusive_dmerge(njets) / Q2();}
00146 
00147   /// same as exclusive_dmerge_max, but normalised to squared total energy
00148   double exclusive_ymerge_max (int njets) const {return exclusive_dmerge_max(njets)/Q2();}
00149 
00150   /// the number of exclusive jets at the given ycut
00151   int n_exclusive_jets_ycut (double ycut) const {return n_exclusive_jets(ycut*Q2());}
00152 
00153   /// the exclusive jets obtained at the given ycut
00154   std::vector<PseudoJet> exclusive_jets_ycut (double ycut) const {
00155     int njets = n_exclusive_jets_ycut(ycut);
00156     return exclusive_jets(njets);
00157   }
00158 
00159 
00160   //int n_exclusive_jets (const PseudoJet & jet, const double & dcut) const;
00161 
00162   /// return a vector of all subjets of the current jet (in the sense
00163   /// of the exclusive algorithm) that would be obtained when running
00164   /// the algorithm with the given dcut. 
00165   ///
00166   /// Time taken is O(m ln m), where m is the number of subjets that
00167   /// are found. If m gets to be of order of the total number of
00168   /// constituents in the jet, this could be substantially slower than
00169   /// just getting that list of constituents.
00170   std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet, 
00171                                             const double & dcut) const;
00172 
00173   /// return the size of exclusive_subjets(...); still n ln n with same
00174   /// coefficient, but marginally more efficient than manually taking
00175   /// exclusive_subjets.size()
00176   int n_exclusive_subjets(const PseudoJet & jet, 
00177                           const double & dcut) const;
00178 
00179   /// return the list of subjets obtained by unclustering the supplied
00180   /// jet down to nsub subjets. Throws an error if there are fewer than
00181   /// nsub particles in the jet.
00182   ///
00183   /// This requires nsub ln nsub time
00184   std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet, 
00185                                             int nsub) const;
00186 
00187   /// return the list of subjets obtained by unclustering the supplied
00188   /// jet down to nsub subjets (or all constituents if there are fewer
00189   /// than nsub).
00190   ///
00191   /// This requires nsub ln nsub time
00192   std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet & jet, 
00193                                                   int nsub) const;
00194 
00195   /// return the dij that was present in the merging nsub+1 -> nsub 
00196   /// subjets inside this jet.
00197   ///
00198   /// Returns 0 if there were nsub or fewer constituents in the jet.
00199   double exclusive_subdmerge(const PseudoJet & jet, int nsub) const;
00200 
00201   /// return the maximum dij that occurred in the whole event at the
00202   /// stage that the nsub+1 -> nsub merge of subjets occurred inside 
00203   /// this jet.
00204   ///
00205   /// Returns 0 if there were nsub or fewer constituents in the jet.
00206   double exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const;
00207 
00208   //std::vector<PseudoJet> exclusive_jets (const PseudoJet & jet, 
00209   //                                       const int & njets) const;
00210   //double exclusive_dmerge (const PseudoJet & jet, const int & njets) const;
00211 
00212   /// returns the sum of all energies in the event (relevant mainly for e+e-)
00213   double Q() const {return _Qtot;}
00214   /// return Q()^2
00215   double Q2() const {return _Qtot*_Qtot;}
00216 
00217   /// returns true iff the object is included in the jet. 
00218   ///
00219   /// NB: this is only sensible if the object is already registered
00220   /// within the cluster sequence, so you cannot use it with an input
00221   /// particle to the CS (since the particle won't have the history
00222   /// index set properly).
00223   ///
00224   /// For nice clustering structures it should run in O(ln(N)) time
00225   /// but in worst cases (certain cone plugins) it can take O(n) time,
00226   /// where n is the number of particles in the jet.
00227   bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;
00228 
00229   /// if the jet has parents in the clustering, it returns true
00230   /// and sets parent1 and parent2 equal to them.
00231   ///
00232   /// if it has no parents it returns false and sets parent1 and
00233   /// parent2 to zero
00234   bool has_parents(const PseudoJet & jet, PseudoJet & parent1, 
00235                PseudoJet & parent2) const;
00236 
00237   /// if the jet has a child then return true and give the child jet
00238   /// otherwise return false and set the child to zero
00239   bool has_child(const PseudoJet & jet, PseudoJet & child) const;
00240 
00241   /// Version of has_child that sets a pointer to the child if the child
00242   /// exists;
00243   bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;
00244 
00245   /// if this jet has a child (and so a partner) return true
00246   /// and give the partner, otherwise return false and set the
00247   /// partner to zero
00248   bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;
00249 
00250   
00251   /// return a vector of the particles that make up jet
00252   std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
00253 
00254 
00255   /// output the supplied vector of jets in a format that can be read
00256   /// by an appropriate root script; the format is:
00257   /// jet-n jet-px jet-py jet-pz jet-E 
00258   ///   particle-n particle-rap particle-phi particle-pt
00259   ///   particle-n particle-rap particle-phi particle-pt
00260   ///   ...
00261   /// #END
00262   /// ... [i.e. above repeated]
00263   void print_jets_for_root(const std::vector<PseudoJet> & jets, 
00264                            std::ostream & ostr = std::cout) const;
00265 
00266   /// print jets for root to the file labelled filename, with an
00267   /// optional comment at the beginning
00268   void print_jets_for_root(const std::vector<PseudoJet> & jets, 
00269                            const std::string & filename,
00270                            const std::string & comment = "") const;
00271 
00272 // Not yet. Perhaps in a future release.
00273 //   /// print out all inclusive jets with pt > ptmin
00274 //   virtual void print_jets (const double & ptmin=0.0) const;
00275 
00276   /// add on to subjet_vector the constituents of jet (for internal use mainly)
00277   void add_constituents (const PseudoJet & jet, 
00278                          std::vector<PseudoJet> & subjet_vector) const;
00279 
00280   /// return the enum value of the strategy used to cluster the event
00281   inline Strategy strategy_used () const {return _strategy;}
00282 
00283   /// return the name of the strategy used to cluster the event
00284   std::string strategy_string () const {return strategy_string(_strategy);}
00285 
00286   /// return the name of the strategy associated with the enum strategy_in
00287   std::string strategy_string (Strategy strategy_in) const;
00288 
00289 
00290   /// return a reference to the jet definition
00291   const JetDefinition & jet_def() const {return _jet_def;}
00292 
00293   /// by calling this routine you tell the ClusterSequence to delete
00294   /// itself when all the Pseudojets associated with it have gone out
00295   /// of scope. 
00296   ///
00297   /// At the time you call this, there must be at least one jet or
00298   /// other object outside the CS that is associated with the CS
00299   /// (e.g. the result of inclusive_jets()).
00300   ///
00301   /// NB: after having made this call, the user is still allowed to
00302   /// delete the CS or let it go out of scope. Jets associated with it
00303   /// will then simply not be able to access their substructure after
00304   /// that point.
00305   void delete_self_when_unused();
00306 
00307   /// return true if the object has been told to delete itself
00308   /// when unused
00309   bool will_delete_self_when_unused() const {return _deletes_self_when_unused;}
00310 
00311   /// tell the ClusterSequence it's about to be self deleted (internal use only)
00312   void signal_imminent_self_deletion() const;
00313 
00314   /// returns the scale associated with a jet as required for this
00315   /// clustering algorithm (kt^2 for the kt-algorithm, 1 for the 
00316   /// Cambridge algorithm). [May become virtual at some point]
00317   double jet_scale_for_algorithm(const PseudoJet & jet) const;
00318 
00319   ///
00320 
00321   //----- next follow functions designed specifically for plugins, which
00322   //      may only be called when plugin_activated() returns true
00323 
00324   /// record the fact that there has been a recombination between
00325   /// jets()[jet_i] and jets()[jet_k], with the specified dij, and
00326   /// return the index (newjet_k) allocated to the new jet, whose
00327   /// momentum is assumed to be the 4-vector sum of that of jet_i and
00328   /// jet_j
00329   void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, 
00330                                       int & newjet_k) {
00331     assert(plugin_activated());
00332     _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
00333   }
00334 
00335   /// as for the simpler variant of plugin_record_ij_recombination,
00336   /// except that the new jet is attributed the momentum and
00337   /// user_index of newjet
00338   void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, 
00339                                       const PseudoJet & newjet, 
00340                                       int & newjet_k);
00341 
00342   /// record the fact that there has been a recombination between
00343   /// jets()[jet_i] and the beam, with the specified diB; when looking
00344   /// for inclusive jets, any iB recombination will returned to the user 
00345   /// as a jet.
00346   void plugin_record_iB_recombination(int jet_i, double diB) {
00347     assert(plugin_activated());
00348     _do_iB_recombination_step(jet_i, diB);
00349   }
00350 
00351   /// @ingroup extra_info
00352   /// \class Extras
00353   /// base class to store extra information that plugins may provide
00354   /// 
00355   /// a class intended to serve as a base in case a plugin needs to
00356   /// associate extra information with a ClusterSequence (see
00357   /// SISConePlugin.* for an example).
00358   class Extras {
00359   public:
00360     virtual ~Extras() {}
00361     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";}
00362   };
00363 
00364   /// the plugin can associate some extra information with the
00365   /// ClusterSequence object by calling this function
00366   inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) {
00367     //_extras = extras_in;
00368     _extras.reset(extras_in.release());
00369   }
00370 
00371   /// returns true when the plugin is allowed to run the show.
00372   inline bool plugin_activated() const {return _plugin_activated;}
00373 
00374   /// returns a pointer to the extras object (may be null)
00375   const Extras * extras() const {return _extras.get();}
00376 
00377   /// allows a plugin to run a templated clustering (nearest-neighbour heuristic)
00378   ///
00379   /// This has N^2 behaviour on "good" distance, but a worst case behaviour
00380   /// of N^3 (and many algs trigger the worst case behaviour)
00381   ///
00382   /// 
00383   /// For more details on how this works, see GenBriefJet below
00384   template<class GBJ> void plugin_simple_N2_cluster () {
00385     assert(plugin_activated());
00386     _simple_N2_cluster<GBJ>();
00387   }
00388 
00389 
00390 public:
00391 //DEP   /// set the default (static) jet finder across all current and future
00392 //DEP   /// ClusterSequence objects -- deprecated and obsolescent (i.e. may be
00393 //DEP   /// suppressed in a future release).
00394 //DEP   static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
00395 //DEP   /// same as above for backward compatibility
00396 //DEP   static void set_jet_finder (JetAlgorithm jet_algorithm)    {_default_jet_algorithm = jet_algorithm;}
00397 
00398 
00399   /// \ingroup extra_info
00400   /// \struct history_element
00401   /// a single element in the clustering history
00402   /// 
00403   /// (see vector _history below).
00404   struct history_element{
00405     int parent1; /// index in _history where first parent of this jet
00406                  /// was created (InexistentParent if this jet is an
00407                  /// original particle)
00408 
00409     int parent2; /// index in _history where second parent of this jet
00410                  /// was created (InexistentParent if this jet is an
00411                  /// original particle); BeamJet if this history entry
00412                  /// just labels the fact that the jet has recombined
00413                  /// with the beam)
00414 
00415     int child;   /// index in _history where the current jet is
00416                  /// recombined with another jet to form its child. It
00417                  /// is Invalid if this jet does not further
00418                  /// recombine.
00419 
00420     int jetp_index; /// index in the _jets vector where we will find the
00421                  /// PseudoJet object corresponding to this jet
00422                  /// (i.e. the jet created at this entry of the
00423                  /// history). NB: if this element of the history
00424                  /// corresponds to a beam recombination, then
00425                  /// jetp_index=Invalid.
00426 
00427     double dij;  /// the distance corresponding to the recombination
00428                  /// at this stage of the clustering.
00429 
00430     double max_dij_so_far; /// the largest recombination distance seen
00431                            /// so far in the clustering history.
00432   };
00433 
00434   enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
00435 
00436   /// allow the user to access the internally stored _jets() array,
00437   /// which contains both the initial particles and the various
00438   /// intermediate and final stages of recombination.
00439   ///
00440   /// The first n_particles() entries are the original particles,
00441   /// in the order in which they were supplied to the ClusterSequence
00442   /// constructor. It can be useful to access them for example when
00443   /// examining whether a given input object is part of a specific
00444   /// jet, via the objects_in_jet(...) member function (which only takes
00445   /// PseudoJets that are registered in the ClusterSequence).
00446   ///
00447   /// One of the other (internal uses) is related to the fact
00448   /// because we don't seem to be able to access protected elements of
00449   /// the class for an object that is not "this" (at least in case where
00450   /// "this" is of a slightly different kind from the object, both
00451   /// derived from ClusterSequence).
00452   const std::vector<PseudoJet> & jets()    const;
00453 
00454   /// allow the user to access the raw internal history.
00455   ///
00456   /// This is present (as for jets()) in part so that protected
00457   /// derived classes can access this information about other
00458   /// ClusterSequences.
00459   ///
00460   /// A user who wishes to follow the details of the ClusterSequence
00461   /// can also make use of this information (and should consult the
00462   /// history_element documentation for more information), but should
00463   /// be aware that these internal structures may evolve in future
00464   /// FastJet versions.
00465   const std::vector<history_element> & history() const;
00466 
00467   /// returns the number of particles that were provided to the
00468   /// clustering algorithm (helps the user find their way around the
00469   /// history and jets objects if they weren't paying attention
00470   /// beforehand).
00471   unsigned int n_particles() const;
00472 
00473   /// returns a vector of size n_particles() which indicates, for 
00474   /// each of the initial particles (in the order in which they were
00475   /// supplied), which of the supplied jets it belongs to; if it does
00476   /// not belong to any of the supplied jets, the index is set to -1;
00477   std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
00478 
00479   /// routine that returns an order in which to read the history
00480   /// such that clusterings that lead to identical jet compositions
00481   /// but different histories (because of degeneracies in the
00482   /// clustering order) will have matching constituents for each
00483   /// matching entry in the unique_history_order.
00484   ///
00485   /// The order has the property that an entry's parents will always
00486   /// appear prior to that entry itself. 
00487   ///
00488   /// Roughly speaking the order is such that we first provide all
00489   /// steps that lead to the final jet containing particle 1; then we
00490   /// have the steps that lead to reconstruction of the jet containing
00491   /// the next-lowest-numbered unclustered particle, etc...
00492   /// [see GPS CCN28-12 for more info -- of course a full explanation
00493   /// here would be better...]
00494   std::vector<int> unique_history_order() const;
00495 
00496   /// return the set of particles that have not been clustered. For 
00497   /// kt and cam/aachen algorithms this should always be null, but for
00498   /// cone type algorithms it can be non-null;
00499   std::vector<PseudoJet> unclustered_particles() const;
00500 
00501   /// Return the list of pseudojets in the ClusterSequence that do not
00502   /// have children (and are not among the inclusive jets). They may
00503   /// result from a clustering step or may be one of the pseudojets
00504   /// returned by unclustered_particles().
00505   std::vector<PseudoJet> childless_pseudojets() const;
00506 
00507   /// returns true if the object (jet or particle) is contained by (ie
00508   /// belongs to) this cluster sequence.
00509   ///
00510   /// Tests performed: if thejet's interface is this cluster sequence
00511   /// and its cluster history index is in a consistent range.
00512   bool contains(const PseudoJet & object) const;
00513 
00514   /// transfer the sequence contained in other_seq into our own;
00515   /// any plugin "extras" contained in the from_seq will be lost
00516   /// from there.
00517   ///
00518   /// It also sets the ClusterSequence pointers of the PseudoJets in
00519   /// the history to point to this ClusterSequence
00520   ///
00521   /// When specified, the second argument is an action that will be
00522   /// applied on every jets in the resulting ClusterSequence
00523   void transfer_from_sequence(const ClusterSequence & from_seq,
00524                               const FunctionOfPseudoJet<PseudoJet> * action_on_jets = 0);
00525 
00526   /// retrieve a shared pointer to the wrapper to this ClusterSequence
00527   ///
00528   /// this may turn useful if you want to track when this
00529   /// ClusterSequence goes out of scope
00530   const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const{
00531     return _structure_shared_ptr;
00532   }
00533 
00534   /// the structure type associated with a jet belonging to a ClusterSequence
00535   typedef ClusterSequenceStructure StructureType;
00536 
00537   /// This is the function that is automatically called during
00538   /// clustering to print the FastJet banner. Only the first call to
00539   /// this function will result in the printout of the banner. Users
00540   /// may wish to call this function themselves, during the
00541   /// initialization phase of their program, in order to ensure that
00542   /// the banner appears before other output. This call will not
00543   /// affect 3rd-party banners, e.g. those from plugins.
00544   static void print_banner();
00545 
00546   /// \cond internal_doc
00547   //  [this line must be left as is to hide the doxygen comment]
00548   /// A call to this function modifies the stream used to print
00549   /// banners (by default cout). If a null pointer is passed, banner
00550   /// printout is suppressed. This affects all banners, including
00551   /// those from plugins.
00552   ///
00553   /// Please note that if you distribute 3rd party code
00554   /// that links with FastJet, that 3rd party code must not
00555   /// use this call turn off the printing of FastJet banners
00556   /// by default. This requirement reflects the spirit of
00557   /// clause 2c of the GNU Public License (v2), under which
00558   /// FastJet and its plugins are distributed.
00559   static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
00560   //  [this line must be left as is to hide the doxygen comment]
00561   /// \endcond
00562 
00563   /// returns a pointer to the stream to be used to print banners
00564   /// (cout by default). This function is used by plugins to determine
00565   /// where to direct their banners. Plugins should properly handle
00566   /// the case where the pointer is null.
00567   static std::ostream * fastjet_banner_stream() {return _fastjet_banner_ostr;}
00568 
00569 private:
00570   /// \cond internal_doc
00571 
00572   /// contains the actual stream to use for banners 
00573   static std::ostream * _fastjet_banner_ostr;
00574 
00575   /// \endcond
00576 
00577 protected:
00578 //DEP  static JetAlgorithm _default_jet_algorithm;
00579   JetDefinition _jet_def;
00580 
00581   /// transfer the vector<L> of input jets into our own vector<PseudoJet>
00582   /// _jets (with some reserved space for future growth).
00583   template<class L> void _transfer_input_jets(
00584                                      const std::vector<L> & pseudojets);
00585 
00586   /// This is what is called to do all the initialisation and
00587   /// then run the clustering (may be called by various constructors).
00588   /// It assumes _jets contains the momenta to be clustered.
00589   void _initialise_and_run (const JetDefinition & jet_def,
00590                             const bool & writeout_combinations);
00591 
00592   //// this performs the initialisation, minus the option-decanting
00593   //// stage; for low multiplicity, initialising a few things in the
00594   //// constructor, calling the decant_options_partial() and then this
00595   //// is faster than going through _initialise_and_run.
00596   void _initialise_and_run_no_decant();
00597 
00598 //DEP   /// This is an alternative routine for initialising and running the
00599 //DEP   /// clustering, provided for legacy purposes. The jet finder is that
00600 //DEP   /// specified in the static member _default_jet_algorithm.
00601 //DEP   void _initialise_and_run (const double & R,
00602 //DEP                       const Strategy & strategy,
00603 //DEP                       const bool & writeout_combinations);
00604 
00605   /// fills in the various member variables with "decanted" options from
00606   /// the jet_definition and writeout_combinations variables
00607   void _decant_options(const JetDefinition & jet_def,
00608                        const bool & writeout_combinations);
00609 
00610   /// assuming that the jet definition, writeout_combinations and
00611   /// _structure_shared_ptr have been set (e.g. in an initialiser list
00612   /// in the constructor), it handles the remaining decanting of
00613   /// options.
00614   void _decant_options_partial();
00615 
00616   /// fill out the history (and jet cross refs) related to the initial
00617   /// set of jets (assumed already to have been "transferred"),
00618   /// without any clustering
00619   void _fill_initial_history();
00620 
00621   /// carry out the recombination between the jets numbered jet_i and
00622   /// jet_j, at distance scale dij; return the index newjet_k of the
00623   /// result of the recombination of i and j.
00624   void _do_ij_recombination_step(const int & jet_i, const int & jet_j, 
00625                                  const double & dij, int & newjet_k);
00626 
00627   /// carry out an recombination step in which _jets[jet_i] merges with
00628   /// the beam, 
00629   void _do_iB_recombination_step(const int & jet_i, const double & diB);
00630 
00631   /// every time a jet is added internally during clustering, this
00632   /// should be called to set the jet's structure shared ptr to point
00633   /// to the CS (and the count of internally associated objects is
00634   /// also updated). This should not be called outside construction of
00635   /// a CS object.
00636   void _set_structure_shared_ptr(PseudoJet & j);
00637 
00638   /// make sure that the CS's internal tally of the use count matches
00639   /// that of the _structure_shared_ptr
00640   void _update_structure_use_count();
00641   
00642 
00643   /// This contains the physical PseudoJets; for each PseudoJet one
00644   /// can find the corresponding position in the _history by looking
00645   /// at _jets[i].cluster_hist_index().
00646   std::vector<PseudoJet> _jets;
00647 
00648 
00649   /// this vector will contain the branching history; for each stage,
00650   /// _history[i].jetp_index indicates where to look in the _jets
00651   /// vector to get the physical PseudoJet.
00652   std::vector<history_element> _history;
00653 
00654   /// set subhist to be a set pointers to history entries corresponding to the
00655   /// subjets of this jet; one stops going working down through the
00656   /// subjets either when 
00657   ///   - there is no further to go
00658   ///   - one has found maxjet entries
00659   ///   - max_dij_so_far <= dcut
00660   /// By setting maxjet=0 one can use just dcut; by setting dcut<0
00661   /// one can use jet maxjet
00662   void get_subhist_set(std::set<const history_element*> & subhist,
00663                        const  PseudoJet & jet, double dcut, int maxjet) const;
00664 
00665   bool _writeout_combinations;
00666   int  _initial_n;
00667   double _Rparam, _R2, _invR2;
00668   double _Qtot;
00669   Strategy    _strategy;
00670   JetAlgorithm  _jet_algorithm;
00671 
00672   SharedPtr<PseudoJetStructureBase> _structure_shared_ptr; //< will actually be of type ClusterSequenceStructure
00673   int _structure_use_count_after_construction; //< info of use when CS handles its own memory
00674   /// if true then the CS will delete itself when the last external
00675   /// object referring to it disappears. It is mutable so as to ensure
00676   /// that signal_imminent_self_deletion() [const] can make relevant
00677   /// changes.
00678   mutable bool _deletes_self_when_unused;
00679 
00680  private:
00681 
00682   bool _plugin_activated;
00683   //std::auto_ptr<Extras> _extras; // things the plugin might want to add
00684   SharedPtr<Extras> _extras; // things the plugin might want to add
00685 
00686   void _really_dumb_cluster ();
00687   void _delaunay_cluster ();
00688   //void _simple_N2_cluster ();
00689   template<class BJ> void _simple_N2_cluster ();
00690   void _tiled_N2_cluster ();
00691   void _faster_tiled_N2_cluster ();
00692 
00693   //
00694   void _minheap_faster_tiled_N2_cluster();
00695 
00696   // things needed specifically for Cambridge with Chan's 2D closest
00697   // pairs method
00698   void _CP2DChan_cluster();
00699   void _CP2DChan_cluster_2pi2R ();
00700   void _CP2DChan_cluster_2piMultD ();
00701   void _CP2DChan_limited_cluster(double D);
00702   void _do_Cambridge_inclusive_jets();
00703 
00704   // NSqrtN method for C/A
00705   void _fast_NsqrtN_cluster();
00706 
00707   void _add_step_to_history(const int & step_number, const int & parent1, 
00708                                const int & parent2, const int & jetp_index,
00709                                const double & dij);
00710 
00711   /// internal routine associated with the construction of the unique
00712   /// history order (following children in the tree)
00713   void _extract_tree_children(int pos, std::valarray<bool> &, 
00714                 const std::valarray<int> &, std::vector<int> &) const;
00715 
00716   /// internal routine associated with the construction of the unique
00717   /// history order (following parents in the tree)
00718   void _extract_tree_parents (int pos, std::valarray<bool> &, 
00719                 const std::valarray<int> &,  std::vector<int> &) const;
00720 
00721 
00722   // these will be useful shorthands in the Voronoi-based code
00723   typedef std::pair<int,int> TwoVertices;
00724   typedef std::pair<double,TwoVertices> DijEntry;
00725   typedef std::multimap<double,TwoVertices> DistMap;
00726 
00727   /// currently used only in the Voronoi based code
00728   void _add_ktdistance_to_map(const int & ii, 
00729                               DistMap & DijMap,
00730                               const DynamicNearestNeighbours * DNN);
00731 
00732 
00733   /// will be set by default to be true for the first run
00734   static bool _first_time;
00735 
00736   /// record the number of warnings provided about the exclusive
00737   /// algorithm -- so that we don't print it out more than a few
00738   /// times.
00739   static int _n_exclusive_warnings;
00740 
00741   /// the limited warning member for notification of user that 
00742   /// their requested strategy has been overridden (usually because
00743   /// they have R>2pi and not all strategies work then)
00744   static LimitedWarning _changed_strategy_warning;
00745 
00746   //----------------------------------------------------------------------
00747   /// the fundamental structure which contains the minimal info about
00748   /// a jet, as needed for our plain N^2 algorithm -- the idea is to
00749   /// put all info that will be accessed N^2 times into an array of
00750   /// BriefJets...
00751   struct BriefJet {
00752     double     eta, phi, kt2, NN_dist;
00753     BriefJet * NN;
00754     int        _jets_index;
00755   };
00756 
00757 
00758   /// structure analogous to BriefJet, but with the extra information
00759   /// needed for dealing with tiles
00760   class TiledJet {
00761   public:
00762     double     eta, phi, kt2, NN_dist;
00763     TiledJet * NN, *previous, * next; 
00764     int        _jets_index, tile_index, diJ_posn;
00765     // routines that are useful in the minheap version of tiled
00766     // clustering ("misuse" the otherwise unused diJ_posn, so as
00767     // to indicate whether jets need to have their minheap entries
00768     // updated).
00769     inline void label_minheap_update_needed() {diJ_posn = 1;}
00770     inline void label_minheap_update_done()   {diJ_posn = 0;}
00771     inline bool minheap_update_needed() const {return diJ_posn==1;}
00772   };
00773 
00774   //-- some of the functions that follow are templates and will work
00775   //as well for briefjet and tiled jets
00776 
00777   /// set the kinematic and labelling info for jeta so that it corresponds
00778   /// to _jets[_jets_index]
00779   template <class J> void _bj_set_jetinfo( J * const jet, 
00780                                                  const int _jets_index) const;
00781 
00782   /// "remove" this jet, which implies updating links of neighbours and
00783   /// perhaps modifying the tile structure
00784   void _bj_remove_from_tiles( TiledJet * const jet) const;
00785 
00786   /// return the distance between two BriefJet objects
00787   template <class J> double _bj_dist(const J * const jeta, 
00788                         const J * const jetb) const;
00789 
00790   // return the diJ (multiplied by _R2) for this jet assuming its NN
00791   // info is correct
00792   template <class J> double _bj_diJ(const J * const jeta) const;
00793 
00794   /// for testing purposes only: if in the range head--tail-1 there is a
00795   /// a jet which corresponds to hist_index in the history, then
00796   /// return a pointer to that jet; otherwise return tail.
00797   template <class J> inline J * _bj_of_hindex(
00798                           const int hist_index, 
00799                           J * const head, J * const tail) 
00800     const {
00801     J * res;
00802     for(res = head; res<tail; res++) {
00803       if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
00804     }
00805     return res;
00806   }
00807 
00808 
00809   //-- remaining functions are different in various cases, so we
00810   //   will use templates but are not sure if they're useful...
00811 
00812   /// updates (only towards smaller distances) the NN for jeta without checking
00813   /// whether in the process jeta itself might be a new NN of one of
00814   /// the jets being scanned -- span the range head to tail-1 with
00815   /// assumption that jeta is not contained in that range
00816   template <class J> void _bj_set_NN_nocross(J * const jeta, 
00817             J * const head, const J * const tail) const;
00818 
00819   /// reset the NN for jeta and DO check whether in the process jeta
00820   /// itself might be a new NN of one of the jets being scanned --
00821   /// span the range head to tail-1 with assumption that jeta is not
00822   /// contained in that range
00823   template <class J> void _bj_set_NN_crosscheck(J * const jeta, 
00824             J * const head, const J * const tail) const;
00825   
00826 
00827 
00828   /// number of neighbours that a tile will have (rectangular geometry
00829   /// gives 9 neighbours).
00830   static const int n_tile_neighbours = 9;
00831   //----------------------------------------------------------------------
00832   /// The fundamental structures to be used for the tiled N^2 algorithm
00833   /// (see CCN27-44 for some discussion of pattern of tiling)
00834   struct Tile {
00835     /// pointers to neighbouring tiles, including self
00836     Tile *   begin_tiles[n_tile_neighbours]; 
00837     /// neighbouring tiles, excluding self
00838     Tile **  surrounding_tiles; 
00839     /// half of neighbouring tiles, no self
00840     Tile **  RH_tiles;  
00841     /// just beyond end of tiles
00842     Tile **  end_tiles; 
00843     /// start of list of BriefJets contained in this tile
00844     TiledJet * head;    
00845     /// sometimes useful to be able to tag a tile
00846     bool     tagged;    
00847   };
00848   std::vector<Tile> _tiles;
00849   double _tiles_eta_min, _tiles_eta_max;
00850   double _tile_size_eta, _tile_size_phi;
00851   int    _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
00852 
00853   // reasonably robust return of tile index given ieta and iphi, in particular
00854   // it works even if iphi is negative
00855   inline int _tile_index (int ieta, int iphi) const {
00856     // note that (-1)%n = -1 so that we have to add _n_tiles_phi
00857     // before performing modulo operation
00858     return (ieta-_tiles_ieta_min)*_n_tiles_phi
00859                   + (iphi+_n_tiles_phi) % _n_tiles_phi;
00860   }
00861 
00862   // routines for tiled case, including some overloads of the plain
00863   // BriefJet cases
00864   int  _tile_index(const double & eta, const double & phi) const;
00865   void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
00866   void  _bj_remove_from_tiles(TiledJet * const jet);
00867   void _initialise_tiles();
00868   void _print_tiles(TiledJet * briefjets ) const;
00869   void _add_neighbours_to_tile_union(const int tile_index, 
00870                  std::vector<int> & tile_union, int & n_near_tiles) const;
00871   void _add_untagged_neighbours_to_tile_union(const int tile_index, 
00872                  std::vector<int> & tile_union, int & n_near_tiles);
00873 
00874 
00875   //----------------------------------------------------------------------
00876   /// fundamental structure for e+e- clustering
00877   struct EEBriefJet {
00878     double NN_dist;  // obligatorily present
00879     double kt2;      // obligatorily present == E^2 in general
00880     EEBriefJet * NN; // must be present too
00881     int    _jets_index; // must also be present!
00882     //...........................................................
00883     double nx, ny, nz;  // our internal storage for fast distance calcs
00884   };
00885 
00886   /// to help instantiation (fj 2.4.0; did not quite work on gcc 33 and os x 10.3?)
00887   //void _dummy_N2_cluster_instantiation();
00888 
00889 
00890   /// to avoid issues with template instantiation (OS X 10.3, gcc 3.3)
00891   void _simple_N2_cluster_BriefJet();
00892   /// to avoid issues with template instantiation (OS X 10.3, gcc 3.3)
00893   void _simple_N2_cluster_EEBriefJet();
00894 };
00895 
00896 
00897 //**********************************************************************
00898 //**************    START   OF   INLINE   MATERIAL    ******************
00899 //**********************************************************************
00900 
00901 
00902 //----------------------------------------------------------------------
00903 // Transfer the initial jets into our internal structure
00904 template<class L> void ClusterSequence::_transfer_input_jets(
00905                                        const std::vector<L> & pseudojets) {
00906 
00907   // this will ensure that we can point to jets without difficulties
00908   // arising.
00909   _jets.reserve(pseudojets.size()*2);
00910 
00911   // insert initial jets this way so that any type L that can be
00912   // converted to a pseudojet will work fine (basically PseudoJet
00913   // and any type that has [] subscript access to the momentum
00914   // components, such as CLHEP HepLorentzVector).
00915   for (unsigned int i = 0; i < pseudojets.size(); i++) {
00916     _jets.push_back(pseudojets[i]);}
00917   
00918 }
00919 
00920 // //----------------------------------------------------------------------
00921 // // initialise from some generic type... Has to be made available
00922 // // here in order for it the template aspect of it to work...
00923 // template<class L> ClusterSequence::ClusterSequence (
00924 //                                const std::vector<L> & pseudojets,
00925 //                                const double & R,
00926 //                                const Strategy & strategy,
00927 //                                const bool & writeout_combinations) {
00928 // 
00929 //   // transfer the initial jets (type L) into our own array
00930 //   _transfer_input_jets(pseudojets);
00931 // 
00932 //   // run the clustering
00933 //   _initialise_and_run(R,strategy,writeout_combinations);
00934 // }
00935 
00936 
00937 //----------------------------------------------------------------------
00938 /// constructor of a jet-clustering sequence from a vector of
00939 /// four-momenta, with the jet definition specified by jet_def
00940 template<class L> ClusterSequence::ClusterSequence (
00941                                   const std::vector<L> & pseudojets,
00942                                   const JetDefinition & jet_def_in,
00943                                   const bool & writeout_combinations) :
00944   _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
00945   _structure_shared_ptr(new ClusterSequenceStructure(this))
00946 {
00947 
00948   // transfer the initial jets (type L) into our own array
00949   _transfer_input_jets(pseudojets);
00950 
00951   // transfer the remaining options
00952   _decant_options_partial();
00953 
00954   // run the clustering
00955   _initialise_and_run_no_decant();
00956 }
00957 
00958 
00959 inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
00960   return _jets;
00961 }
00962 
00963 inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
00964   return _history;
00965 }
00966 
00967 inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
00968 
00969 
00970 
00971 //----------------------------------------------------------------------
00972 template <class J> inline void ClusterSequence::_bj_set_jetinfo(
00973                             J * const jetA, const int _jets_index) const {
00974     jetA->eta  = _jets[_jets_index].rap();
00975     jetA->phi  = _jets[_jets_index].phi_02pi();
00976     jetA->kt2  = jet_scale_for_algorithm(_jets[_jets_index]);
00977     jetA->_jets_index = _jets_index;
00978     // initialise NN info as well
00979     jetA->NN_dist = _R2;
00980     jetA->NN      = NULL;
00981 }
00982 
00983 
00984 
00985 
00986 //----------------------------------------------------------------------
00987 template <class J> inline double ClusterSequence::_bj_dist(
00988                 const J * const jetA, const J * const jetB) const {
00989   double dphi = std::abs(jetA->phi - jetB->phi);
00990   double deta = (jetA->eta - jetB->eta);
00991   if (dphi > pi) {dphi = twopi - dphi;}
00992   return dphi*dphi + deta*deta;
00993 }
00994 
00995 //----------------------------------------------------------------------
00996 template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
00997   double kt2 = jet->kt2;
00998   if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
00999   return jet->NN_dist * kt2;
01000 }
01001 
01002 
01003 //----------------------------------------------------------------------
01004 // set the NN for jet without checking whether in the process you might
01005 // have discovered a new nearest neighbour for another jet
01006 template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
01007                  J * const jet, J * const head, const J * const tail) const {
01008   double NN_dist = _R2;
01009   J * NN  = NULL;
01010   if (head < jet) {
01011     for (J * jetB = head; jetB != jet; jetB++) {
01012       double dist = _bj_dist(jet,jetB);
01013       if (dist < NN_dist) {
01014         NN_dist = dist;
01015         NN = jetB;
01016       }
01017     }
01018   }
01019   if (tail > jet) {
01020     for (J * jetB = jet+1; jetB != tail; jetB++) {
01021       double dist = _bj_dist(jet,jetB);
01022       if (dist < NN_dist) {
01023         NN_dist = dist;
01024         NN = jetB;
01025       }
01026     }
01027   }
01028   jet->NN = NN;
01029   jet->NN_dist = NN_dist;
01030 }
01031 
01032 
01033 //----------------------------------------------------------------------
01034 template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet, 
01035                     J * const head, const J * const tail) const {
01036   double NN_dist = _R2;
01037   J * NN  = NULL;
01038   for (J * jetB = head; jetB != tail; jetB++) {
01039     double dist = _bj_dist(jet,jetB);
01040     if (dist < NN_dist) {
01041       NN_dist = dist;
01042       NN = jetB;
01043     }
01044     if (dist < jetB->NN_dist) {
01045       jetB->NN_dist = dist;
01046       jetB->NN = jet;
01047     }
01048   }
01049   jet->NN = NN;
01050   jet->NN_dist = NN_dist;
01051 }
01052 
01053 FASTJET_END_NAMESPACE
01054 
01055 #endif // __FASTJET_CLUSTERSEQUENCE_HH__
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends