added pdet-ppart over ppart histogram for detector response

[u/mrichter/AliRoot.git] / JETAN / fastjet / fastjet / ClusterSequence.hh

//STARTHEADER
// $Id: ClusterSequence.hh 1435 2009-02-12 21:11:04Z salam $
//
// Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
//  FastJet is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  The algorithms that underlie FastJet have required considerable
//  development and are described in hep-ph/0512210. If you use
//  FastJet as part of work towards a scientific publication, please
//  include a citation to the FastJet paper.
//
//  FastJet is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with FastJet; if not, write to the Free Software
//  Foundation, Inc.:
//      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//----------------------------------------------------------------------
//ENDHEADER


//----------------------------------------------------------------------
// here's where we put the main page for fastjet (as explained in the
// Doxygen faq)
//......................................................................
/*! \mainpage FastJet code documentation
 *
 * These pages provide automatically generated documentation for the 
 * FastJet package.
 * 
 * For further information and normal documentation, see the main <a
 * href="http://www.lpthe.jussieu.fr/~salam/fastjet">FastJet</a> page.
 */
//----------------------------------------------------------------------

#ifndef __FASTJET_CLUSTERSEQUENCE_HH__
#define __FASTJET_CLUSTERSEQUENCE_HH__

#include<vector>
#include<map>
#include "fastjet/internal/DynamicNearestNeighbours.hh"
#include "fastjet/PseudoJet.hh"
#include<memory>
#include<cassert>
#include<iostream>
#include<string>
#include<set>
#include<cmath> // needed to get double std::abs(double)
#include "fastjet/Error.hh"
#include "fastjet/JetDefinition.hh"

namespace fastjet {      // defined in fastjet/internal/base.hh


/// deals with clustering
class ClusterSequence {


 public: 

  /// default constructor
  ClusterSequence () {}

  /// create a clustersequence starting from the supplied set
  /// of pseudojets and clustering them with the long-invariant
  /// kt algorithm (E-scheme recombination) with the supplied
  /// value for R.
  ///
  /// If strategy=DumbN3 a very stupid N^3 algorithm is used for the
  /// clustering; otherwise strategy = NlnN* uses cylinders algorithms
  /// with some number of pi coverage. If writeout_combinations=true a
  /// summary of the recombination sequence is written out
  template<class L> ClusterSequence (const std::vector<L> & pseudojets, 
		   const double & R = 1.0,
		   const Strategy & strategy = Best,
		   const bool & writeout_combinations = false);


  /// create a clustersequence starting from the supplied set
  /// of pseudojets and clustering them with jet definition specified
  /// by jet_def (which also specifies the clustering strategy)
  template<class L> ClusterSequence (
			          const std::vector<L> & pseudojets,
				  const JetDefinition & jet_def,
				  const bool & writeout_combinations = false);

  // virtual ClusterSequence destructor, in case any derived class
  // thinks of needing a destructor at some point
  virtual ~ClusterSequence (); //{}

  // NB: in the routines that follow, for extracting lists of jets, a
  //     list structure might be more efficient, if sometimes a little
  //     more awkward to use (at least for old fortran hands).

  /// return a vector of all jets (in the sense of the inclusive
  /// algorithm) with pt >= ptmin. Time taken should be of the order
  /// of the number of jets returned.
  std::vector<PseudoJet> inclusive_jets (const double & ptmin = 0.0) const;

  /// return the number of jets (in the sense of the exclusive
  /// algorithm) that would be obtained when running the algorithm
  /// with the given dcut.
  int n_exclusive_jets (const double & dcut) const;

  /// return a vector of all jets (in the sense of the exclusive
  /// algorithm) that would be obtained when running the algorithm
  /// with the given dcut.
  std::vector<PseudoJet> exclusive_jets (const double & dcut) const;

  /// return a vector of all jets when the event is clustered (in the
  /// exclusive sense) to exactly njets.
  std::vector<PseudoJet> exclusive_jets (const int & njets) const;

  /// return the dmin corresponding to the recombination that went from
  /// n+1 to n jets (sometimes known as d_{n n+1}).
  double exclusive_dmerge (const int & njets) const;

  /// return the maximum of the dmin encountered during all recombinations 
  /// up to the one that led to an n-jet final state; identical to
  /// exclusive_dmerge, except in cases where the dmin do not increase
  /// monotonically.
  double exclusive_dmerge_max (const int & njets) const;

  /// return the ymin corresponding to the recombination that went from
  /// n+1 to n jets (sometimes known as y_{n n+1}).
  double exclusive_ymerge (int njets) const {return exclusive_dmerge(njets) / Q2();}

  /// same as exclusive_dmerge_max, but normalised to squared total energy
  double exclusive_ymerge_max (int njets) const {return exclusive_ymerge_max(njets)/Q2();}

  /// the number of exclusive jets at the given ycut
  int n_exclusive_jets_ycut (double ycut) const {return n_exclusive_jets(ycut*Q2());}

  /// the exclusive jets obtained at the given ycut
  std::vector<PseudoJet> exclusive_jets_ycut (double ycut) const {
    int njets = n_exclusive_jets_ycut(ycut);
    return exclusive_jets(njets);
  }


  //int n_exclusive_jets (const PseudoJet & jet, const double & dcut) const;

  /// return a vector of all subjets of the current jet (in the sense
  /// of the exclusive algorithm) that would be obtained when running
  /// the algorithm with the given dcut. 
  ///
  /// Time taken is O(m ln m), where m is the number of subjets that
  /// are found. If m gets to be of order of the total number of
  /// constituents in the jet, this could be substantially slower than
  /// just getting that list of constituents.
  std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet, 
                                            const double & dcut) const;

  /// return the size of exclusive_subjets(...); still n ln n with same
  /// coefficient, but marginally more efficient than manually taking
  /// exclusive_subjets.size()
  int n_exclusive_subjets(const PseudoJet & jet, 
                          const double & dcut) const;

  /// return the list of subjets obtained by unclustering the supplied
  /// jet down to n subjets (or all constituents if there are fewer
  /// than n).
  ///
  /// requires n ln n time
  std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet, 
                                            int nsub) const;

  /// return the dij that was present in the merging nsub+1 -> nsub 
  /// subjets inside this jet.
  double exclusive_subdmerge(const PseudoJet & jet, int nsub) const;

  /// return the maximum dij that occurred in the whole event at the
  /// stage that the nsub+1 -> nsub merge of subjets occurred inside 
  /// this jet.
  double exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const;

  //std::vector<PseudoJet> exclusive_jets (const PseudoJet & jet, 
  //                                       const int & njets) const;
  //double exclusive_dmerge (const PseudoJet & jet, const int & njets) const;

  /// returns the sum of all energies in the event (relevant mainly for e+e-)
  double Q() const {return _Q;}
  /// return Q()^2
  double Q2() const {return _Q*_Q;}

  /// returns true iff the object is included in the jet. 
  ///
  /// NB: this is only sensible if the object is already registered
  /// within the cluster sequence, so you cannot use it with an input
  /// particle to the CS (since the particle won't have the history
  /// index set properly).
  ///
  /// For nice clustering structures it should run in O(ln(N)) time
  /// but in worst cases (certain cone plugins) it can take O(n) time,
  /// where n is the number of particles in the jet.
  bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;

  /// if the jet has parents in the clustering, it returns true
  /// and sets parent1 and parent2 equal to them.
  ///
  /// if it has no parents it returns false and sets parent1 and
  /// parent2 to zero
  bool has_parents(const PseudoJet & jet, PseudoJet & parent1, 
               PseudoJet & parent2) const;

  /// if the jet has a child then return true and give the child jet
  /// otherwise return false and set the child to zero
  bool has_child(const PseudoJet & jet, PseudoJet & child) const;

  /// Version of has_child that sets a pointer to the child if the child
  /// exists;
  bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;

  /// if this jet has a child (and so a partner) return true
  /// and give the partner, otherwise return false and set the
  /// partner to zero
  bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;

  
  /// return a vector of the particles that make up jet
  std::vector<PseudoJet> constituents (const PseudoJet & jet) const;


  /// output the supplied vector of jets in a format that can be read
  /// by an appropriate root script; the format is:
  /// jet-n jet-px jet-py jet-pz jet-E 
  ///   particle-n particle-rap particle-phi particle-pt
  ///   particle-n particle-rap particle-phi particle-pt
  ///   ...
  /// #END
  /// ... [i.e. above repeated]
  void print_jets_for_root(const std::vector<PseudoJet> & jets, 
                           std::ostream & ostr = std::cout) const;

  /// print jets for root to the file labelled filename, with an
  /// optional comment at the beginning
  void print_jets_for_root(const std::vector<PseudoJet> & jets, 
                           const std::string & filename,
			   const std::string & comment = "") const;

// Not yet. Perhaps in a future release.
//   /// print out all inclusive jets with pt > ptmin
//   virtual void print_jets (const double & ptmin=0.0) const;

  /// add on to subjet_vector the constituents of jet (for internal use mainly)
  void add_constituents (const PseudoJet & jet, 
			 std::vector<PseudoJet> & subjet_vector) const;

  /// return the enum value of the strategy used to cluster the event
  inline Strategy strategy_used () const {return _strategy;}
  std::string strategy_string () const;

  /// return a reference to the jet definition
  const JetDefinition & jet_def() const {return _jet_def;}

  /// returns the scale associated with a jet as required for this
  /// clustering algorithm (kt^2 for the kt-algorithm, 1 for the 
  /// Cambridge algorithm). [May become virtual at some point]
  double jet_scale_for_algorithm(const PseudoJet & jet) const;

  //----- next follow functions designed specifically for plugins, which
  //      may only be called when plugin_activated() returns true

  /// record the fact that there has been a recombination between
  /// jets()[jet_i] and jets()[jet_k], with the specified dij, and
  /// return the index (newjet_k) allocated to the new jet, whose
  /// momentum is assumed to be the 4-vector sum of that of jet_i and
  /// jet_j
  void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, 
				      int & newjet_k) {
    assert(plugin_activated());
    _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
  }

  /// as for the simpler variant of plugin_record_ij_recombination,
  /// except that the new jet is attributed the momentum and
  /// user_index of newjet
  void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, 
				      const PseudoJet & newjet, 
				      int & newjet_k);

  /// record the fact that there has been a recombination between
  /// jets()[jet_i] and the beam, with the specified diB; when looking
  /// for inclusive jets, any iB recombination will returned to the user 
  /// as a jet.
  void plugin_record_iB_recombination(int jet_i, double diB) {
    assert(plugin_activated());
    _do_iB_recombination_step(jet_i, diB);
  }

  /// a class intended to serve as a base in case a plugin needs to
  /// associate extra information with a ClusterSequence (see
  /// SISConePlugin.* for an example).
  class Extras {
  public:
    virtual ~Extras() {}
    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";}
  };

  /// the plugin can associated some extra information with the
  /// ClusterSequence object by calling this function
  inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) {
    _extras = extras_in;
  }

  /// returns true when the plugin is allowed to run the show.
  inline bool plugin_activated() const {return _plugin_activated;}

  /// returns a pointer to the extras object (may be null)
  const Extras * extras() const {return _extras.get();}

  /// allows a plugin to run a templated clustering (nearest-neighbour heuristic)
  ///
  /// This has N^2 behaviour on "good" distance, but a worst case behaviour
  /// of N^3 (and many algs trigger the worst case behaviour)
  ///
  /// 
  /// For more details on how this works, see GenBriefJet below
  template<class GBJ> void plugin_simple_N2_cluster () {
    assert(plugin_activated());
    _simple_N2_cluster<GBJ>();
  }

  //----------------------------------------------------------------------
  /// class to help with a generic clustering sequence
  class GenBriefJet {
  public:
    /// function that initialises the GenBriefJet given a PseudoJet.
    ///
    /// In a derived class, this member has a responsability to call
    ///
    ///  - set_scale_squared
    ///  - set_geom_iB 
    ///
    /// The clustering will be performed by finding the minimum of
    ///
    ///   diB = scale_squared[i] * geom_iB * _invR2
    ///   dij = min(scale_squared[i],scale_squared[j]) * geom_ij * _invR2
    ///
    virtual void init(const PseudoJet & jet) = 0;
    
    /// Returns the "geometric" part of distance between this jet
    /// and jet_j
    virtual double geom_ij(const GenBriefJet * jet_j) const = 0;

    void set_scale_squared(double scale_squared) {kt2 = scale_squared;}
    void set_geom_iB(double diB) {NN_dist = diB; NN = NULL;}
    
  public: // formally public: but users should think of it as private!
    double NN_dist;  // dij
    double kt2;      // squared scale
    GenBriefJet * NN; // pointer to nearest neighbour
    int    _jets_index; // index of this jet
  };


public:
  /// set the default (static) jet finder across all current and future
  /// ClusterSequence objects -- deprecated and obsolescent (i.e. may be
  /// suppressed in a future release).
  static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
  /// same as above for backward compatibility
  static void set_jet_finder (JetAlgorithm jet_algorithm)    {_default_jet_algorithm = jet_algorithm;}


  /// a single element in the clustering history (see vector _history
  /// below).
  struct history_element{
    int parent1; /// index in _history where first parent of this jet
                 /// was created (InexistentParent if this jet is an
                 /// original particle)

    int parent2; /// index in _history where second parent of this jet
                 /// was created (InexistentParent if this jet is an
                 /// original particle); BeamJet if this history entry
                 /// just labels the fact that the jet has recombined
                 /// with the beam)

    int child;   /// index in _history where the current jet is
		 /// recombined with another jet to form its child. It
		 /// is Invalid if this jet does not further
		 /// recombine.

    int jetp_index; /// index in the _jets vector where we will find the
                 /// PseudoJet object corresponding to this jet
                 /// (i.e. the jet created at this entry of the
                 /// history). NB: if this element of the history
                 /// corresponds to a beam recombination, then
                 /// jetp_index=Invalid.

    double dij;  /// the distance corresponding to the recombination
		 /// at this stage of the clustering.

    double max_dij_so_far; /// the largest recombination distance seen
			   /// so far in the clustering history.
  };

  enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};

  /// allow the user to access the internally stored _jets() array,
  /// which contains both the initial particles and the various
  /// intermediate and final stages of recombination.
  ///
  /// The first n_particles() entries are the original particles,
  /// in the order in which they were supplied to the ClusterSequence
  /// constructor. It can be useful to access them for example when
  /// examining whether a given input object is part of a specific
  /// jet, via the objects_in_jet(...) member function (which only takes
  /// PseudoJets that are registered in the ClusterSequence).
  ///
  /// One of the other (internal uses) is related to the fact
  /// because we don't seem to be able to access protected elements of
  /// the class for an object that is not "this" (at least in case where
  /// "this" is of a slightly different kind from the object, both
  /// derived from ClusterSequence).
  const std::vector<PseudoJet> & jets()    const;

  /// allow the user to access the raw internal history.
  ///
  /// This is present (as for jets()) in part so that protected
  /// derived classes can access this information about other
  /// ClusterSequences.
  ///
  /// A user who wishes to follow the details of the ClusterSequence
  /// can also make use of this information (and should consult the
  /// history_element documentation for more information), but should
  /// be aware that these internal structures may evolve in future
  /// FastJet versions.
  const std::vector<history_element> & history() const;

  /// returns the number of particles that were provided to the
  /// clustering algorithm (helps the user find their way around the
  /// history and jets objects if they weren't paying attention
  /// beforehand).
  unsigned int n_particles() const;

  /// returns a vector of size n_particles() which indicates, for 
  /// each of the initial particles (in the order in which they were
  /// supplied), which of the supplied jets it belongs to; if it does
  /// not belong to any of the supplied jets, the index is set to -1;
  std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;

  /// routine that returns an order in which to read the history
  /// such that clusterings that lead to identical jet compositions
  /// but different histories (because of degeneracies in the
  /// clustering order) will have matching constituents for each
  /// matching entry in the unique_history_order.
  ///
  /// The order has the property that an entry's parents will always
  /// appear prior to that entry itself. 
  ///
  /// Roughly speaking the order is such that we first provide all
  /// steps that lead to the final jet containing particle 1; then we
  /// have the steps that lead to reconstruction of the jet containing
  /// the next-lowest-numbered unclustered particle, etc...
  /// [see GPS CCN28-12 for more info -- of course a full explanation
  /// here would be better...]
  std::vector<int> unique_history_order() const;

  /// return the set of particles that have not been clustered. For 
  /// kt and cam/aachen algorithms this should always be null, but for
  /// cone type algorithms it can be non-null;
  std::vector<PseudoJet> unclustered_particles() const;

  /// transfer the sequence contained in other_seq into our own;
  /// any plugin "extras" contained in the from_seq will be lost
  /// from there.
  void transfer_from_sequence(ClusterSequence & from_seq);


protected:
  static JetAlgorithm _default_jet_algorithm;
  JetDefinition _jet_def;

  /// returns true if the jet has a history index contained within
  /// the range of this CS
  bool _potentially_valid(const PseudoJet & jet) const {
    return jet.cluster_hist_index() >= 0 
      && jet.cluster_hist_index() < int(_history.size());
  }

  /// transfer the vector<L> of input jets into our own vector<PseudoJet>
  /// _jets (with some reserved space for future growth).
  template<class L> void _transfer_input_jets(
                                     const std::vector<L> & pseudojets);

  /// This is the routine that will do all the initialisation and
  /// then run the clustering (may be called by various constructors).
  /// It assumes _jets contains the momenta to be clustered.
  void _initialise_and_run (const JetDefinition & jet_def,
			    const bool & writeout_combinations);

  /// This is an alternative routine for initialising and running the
  /// clustering, provided for legacy purposes. The jet finder is that
  /// specified in the static member _default_jet_algorithm.
  void _initialise_and_run (const double & R,
			    const Strategy & strategy,
			    const bool & writeout_combinations);

  /// fills in the various member variables with "decanted" options from
  /// the jet_definition and writeout_combinations variables
  void _decant_options(const JetDefinition & jet_def,
                       const bool & writeout_combinations);

  /// fill out the history (and jet cross refs) related to the initial
  /// set of jets (assumed already to have been "transferred"),
  /// without any clustering
  void _fill_initial_history();

  /// carry out the recombination between the jets numbered jet_i and
  /// jet_j, at distance scale dij; return the index newjet_k of the
  /// result of the recombination of i and j.
  void _do_ij_recombination_step(const int & jet_i, const int & jet_j, 
				 const double & dij, int & newjet_k);

  /// carry out an recombination step in which _jets[jet_i] merges with
  /// the beam, 
  void _do_iB_recombination_step(const int & jet_i, const double & diB);


  /// This contains the physical PseudoJets; for each PseudoJet one
  /// can find the corresponding position in the _history by looking
  /// at _jets[i].cluster_hist_index().
  std::vector<PseudoJet> _jets;


  /// this vector will contain the branching history; for each stage,
  /// _history[i].jetp_index indicates where to look in the _jets
  /// vector to get the physical PseudoJet.
  std::vector<history_element> _history;

  /// set subhist to be a set pointers to history entries corresponding to the
  /// subjets of this jet; one stops going working down through the
  /// subjets either when 
  ///   - there is no further to go
  ///   - one has found maxjet entries
  ///   - max_dij_so_far <= dcut
  /// By setting maxjet=0 one can use just dcut; by setting dcut<0
  /// one can use jet maxjet
  void get_subhist_set(std::set<const history_element*> & subhist,
                       const  PseudoJet & jet, double dcut, int maxjet) const;

  bool _writeout_combinations;
  int  _initial_n;
  double _Rparam, _R2, _invR2;
  double _Q;
  Strategy    _strategy;
  JetAlgorithm  _jet_algorithm;


 private:

  bool _plugin_activated;
  std::auto_ptr<Extras> _extras; // things the plugin might want to add

  void _really_dumb_cluster ();
  void _delaunay_cluster ();
  //void _simple_N2_cluster ();
  template<class BJ> void _simple_N2_cluster ();
  void _tiled_N2_cluster ();
  void _faster_tiled_N2_cluster ();

  //
  void _minheap_faster_tiled_N2_cluster();

  // things needed specifically for Cambridge with Chan's 2D closest
  // pairs method
  void _CP2DChan_cluster();
  void _CP2DChan_cluster_2pi2R ();
  void _CP2DChan_cluster_2piMultD ();
  void _CP2DChan_limited_cluster(double D);
  void _do_Cambridge_inclusive_jets();

  void _add_step_to_history(const int & step_number, const int & parent1, 
			       const int & parent2, const int & jetp_index,
			       const double & dij);

  /// internal routine associated with the construction of the unique
  /// history order (following children in the tree)
  void _extract_tree_children(int pos, std::valarray<bool> &, 
		const std::valarray<int> &, std::vector<int> &) const;

  /// internal routine associated with the construction of the unique
  /// history order (following parents in the tree)
  void _extract_tree_parents (int pos, std::valarray<bool> &, 
                const std::valarray<int> &,  std::vector<int> &) const;


  // these will be useful shorthands in the Voronoi-based code
  typedef std::pair<int,int> TwoVertices;
  typedef std::pair<double,TwoVertices> DijEntry;
  typedef std::multimap<double,TwoVertices> DistMap;

  /// currently used only in the Voronoi based code
  void _add_ktdistance_to_map(const int & ii, 
			      DistMap & DijMap,
  			      const DynamicNearestNeighbours * DNN);

  /// for making sure the user knows what it is they're running...
  void _print_banner();
  /// will be set by default to be true for the first run
  static bool _first_time;

  /// record the number of warnings provided about the exclusive
  /// algorithm -- so that we don't print it out more than a few
  /// times.
  static int _n_exclusive_warnings;

  //----------------------------------------------------------------------
  /// the fundamental structure which contains the minimal info about
  /// a jet, as needed for our plain N^2 algorithm -- the idea is to
  /// put all info that will be accessed N^2 times into an array of
  /// BriefJets...
  struct BriefJet {
    double     eta, phi, kt2, NN_dist;
    BriefJet * NN;
    int        _jets_index;
  };


  /// structure analogous to BriefJet, but with the extra information
  /// needed for dealing with tiles
  class TiledJet {
  public:
    double     eta, phi, kt2, NN_dist;
    TiledJet * NN, *previous, * next; 
    int        _jets_index, tile_index, diJ_posn;
    // routines that are useful in the minheap version of tiled
    // clustering ("misuse" the otherwise unused diJ_posn, so as
    // to indicate whether jets need to have their minheap entries
    // updated).
    inline void label_minheap_update_needed() {diJ_posn = 1;}
    inline void label_minheap_update_done()   {diJ_posn = 0;}
    inline bool minheap_update_needed() const {return diJ_posn==1;}
  };

  //-- some of the functions that follow are templates and will work
  //as well for briefjet and tiled jets

  /// set the kinematic and labelling info for jeta so that it corresponds
  /// to _jets[_jets_index]
  template <class J> void _bj_set_jetinfo( J * const jet, 
						 const int _jets_index) const;

  /// "remove" this jet, which implies updating links of neighbours and
  /// perhaps modifying the tile structure
  void _bj_remove_from_tiles( TiledJet * const jet) const;

  /// return the distance between two BriefJet objects
  template <class J> double _bj_dist(const J * const jeta, 
			const J * const jetb) const;

  // return the diJ (multiplied by _R2) for this jet assuming its NN
  // info is correct
  template <class J> double _bj_diJ(const J * const jeta) const;

  /// for testing purposes only: if in the range head--tail-1 there is a
  /// a jet which corresponds to hist_index in the history, then
  /// return a pointer to that jet; otherwise return tail.
  template <class J> inline J * _bj_of_hindex(
                          const int hist_index, 
			  J * const head, J * const tail) 
    const {
    J * res;
    for(res = head; res<tail; res++) {
      if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
    }
    return res;
  }


  //-- remaining functions are different in various cases, so we
  //   will use templates but are not sure if they're useful...

  /// updates (only towards smaller distances) the NN for jeta without checking
  /// whether in the process jeta itself might be a new NN of one of
  /// the jets being scanned -- span the range head to tail-1 with
  /// assumption that jeta is not contained in that range
  template <class J> void _bj_set_NN_nocross(J * const jeta, 
            J * const head, const J * const tail) const;

  /// reset the NN for jeta and DO check whether in the process jeta
  /// itself might be a new NN of one of the jets being scanned --
  /// span the range head to tail-1 with assumption that jeta is not
  /// contained in that range
  template <class J> void _bj_set_NN_crosscheck(J * const jeta, 
            J * const head, const J * const tail) const;
  


  /// number of neighbours that a tile will have (rectangular geometry
  /// gives 9 neighbours).
  static const int n_tile_neighbours = 9;
  //----------------------------------------------------------------------
  /// The fundamental structures to be used for the tiled N^2 algorithm
  /// (see CCN27-44 for some discussion of pattern of tiling)
  struct Tile {
    /// pointers to neighbouring tiles, including self
    Tile *   begin_tiles[n_tile_neighbours]; 
    /// neighbouring tiles, excluding self
    Tile **  surrounding_tiles; 
    /// half of neighbouring tiles, no self
    Tile **  RH_tiles;  
    /// just beyond end of tiles
    Tile **  end_tiles; 
    /// start of list of BriefJets contained in this tile
    TiledJet * head;    
    /// sometimes useful to be able to tag a tile
    bool     tagged;    
  };
  std::vector<Tile> _tiles;
  double _tiles_eta_min, _tiles_eta_max;
  double _tile_size_eta, _tile_size_phi;
  int    _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;

  // reasonably robust return of tile index given ieta and iphi, in particular
  // it works even if iphi is negative
  inline int _tile_index (int ieta, int iphi) const {
    // note that (-1)%n = -1 so that we have to add _n_tiles_phi
    // before performing modulo operation
    return (ieta-_tiles_ieta_min)*_n_tiles_phi
                  + (iphi+_n_tiles_phi) % _n_tiles_phi;
  }

  // routines for tiled case, including some overloads of the plain
  // BriefJet cases
  int  _tile_index(const double & eta, const double & phi) const;
  void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
  void  _bj_remove_from_tiles(TiledJet * const jet);
  void _initialise_tiles();
  void _print_tiles(TiledJet * briefjets ) const;
  void _add_neighbours_to_tile_union(const int tile_index, 
		 std::vector<int> & tile_union, int & n_near_tiles) const;
  void _add_untagged_neighbours_to_tile_union(const int tile_index, 
		 std::vector<int> & tile_union, int & n_near_tiles);


  //----------------------------------------------------------------------
  /// fundamental structure for e+e- clustering
  struct EEBriefJet {
    double NN_dist;  // obligatorily present
    double kt2;      // obligatorily present == E^2 in general
    EEBriefJet * NN; // must be present too
    int    _jets_index; // must also be present!
    //...........................................................
    double nx, ny, nz;  // our internal storage for fast distance calcs
  };

  /// to help instantiation
  void _dummy_N2_cluster_instantiation();

};


//**********************************************************************
//**************    START   OF   INLINE   MATERIAL    ******************
//**********************************************************************


//----------------------------------------------------------------------
// Transfer the initial jets into our internal structure
template<class L> void ClusterSequence::_transfer_input_jets(
                                       const std::vector<L> & pseudojets) {

  // this will ensure that we can point to jets without difficulties
  // arising.
  _jets.reserve(pseudojets.size()*2);

  // insert initial jets this way so that any type L that can be
  // converted to a pseudojet will work fine (basically PseudoJet
  // and any type that has [] subscript access to the momentum
  // components, such as CLHEP HepLorentzVector).
  for (unsigned int i = 0; i < pseudojets.size(); i++) {
    _jets.push_back(pseudojets[i]);}
  
}

//----------------------------------------------------------------------
// initialise from some generic type... Has to be made available
// here in order for it the template aspect of it to work...
template<class L> ClusterSequence::ClusterSequence (
			          const std::vector<L> & pseudojets,
				  const double & R,
				  const Strategy & strategy,
				  const bool & writeout_combinations) {

  // transfer the initial jets (type L) into our own array
  _transfer_input_jets(pseudojets);

  // run the clustering
  _initialise_and_run(R,strategy,writeout_combinations);
}


//----------------------------------------------------------------------
/// constructor of a jet-clustering sequence from a vector of
/// four-momenta, with the jet definition specified by jet_def
template<class L> ClusterSequence::ClusterSequence (
			          const std::vector<L> & pseudojets,
				  const JetDefinition & jet_def,
				  const bool & writeout_combinations) {

  // transfer the initial jets (type L) into our own array
  _transfer_input_jets(pseudojets);

  // run the clustering
  _initialise_and_run(jet_def,writeout_combinations);
}


inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
  return _jets;
}

inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
  return _history;
}

inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}



//----------------------------------------------------------------------
template <class J> inline void ClusterSequence::_bj_set_jetinfo(
                            J * const jetA, const int _jets_index) const {
    jetA->eta  = _jets[_jets_index].rap();
    jetA->phi  = _jets[_jets_index].phi_02pi();
    jetA->kt2  = jet_scale_for_algorithm(_jets[_jets_index]);
    jetA->_jets_index = _jets_index;
    // initialise NN info as well
    jetA->NN_dist = _R2;
    jetA->NN      = NULL;
}




//----------------------------------------------------------------------
template <class J> inline double ClusterSequence::_bj_dist(
                const J * const jetA, const J * const jetB) const {
  double dphi = std::abs(jetA->phi - jetB->phi);
  double deta = (jetA->eta - jetB->eta);
  if (dphi > pi) {dphi = twopi - dphi;}
  return dphi*dphi + deta*deta;
}

//----------------------------------------------------------------------
template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
  double kt2 = jet->kt2;
  if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
  return jet->NN_dist * kt2;
}


//----------------------------------------------------------------------
// set the NN for jet without checking whether in the process you might
// have discovered a new nearest neighbour for another jet
template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
                 J * const jet, J * const head, const J * const tail) const {
  double NN_dist = _R2;
  J * NN  = NULL;
  if (head < jet) {
    for (J * jetB = head; jetB != jet; jetB++) {
      double dist = _bj_dist(jet,jetB);
      if (dist < NN_dist) {
	NN_dist = dist;
	NN = jetB;
      }
    }
  }
  if (tail > jet) {
    for (J * jetB = jet+1; jetB != tail; jetB++) {
      double dist = _bj_dist(jet,jetB);
      if (dist < NN_dist) {
	NN_dist = dist;
	NN = jetB;
      }
    }
  }
  jet->NN = NN;
  jet->NN_dist = NN_dist;
}


//----------------------------------------------------------------------
template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet, 
		    J * const head, const J * const tail) const {
  double NN_dist = _R2;
  J * NN  = NULL;
  for (J * jetB = head; jetB != tail; jetB++) {
    double dist = _bj_dist(jet,jetB);
    if (dist < NN_dist) {
      NN_dist = dist;
      NN = jetB;
    }
    if (dist < jetB->NN_dist) {
      jetB->NN_dist = dist;
      jetB->NN = jet;
    }
  }
  jet->NN = NN;
  jet->NN_dist = NN_dist;
}



} // fastjet namespace 

#endif // __FASTJET_CLUSTERSEQUENCE_HH__