ClusterSequence.hh

//STARTHEADER
// $Id: ClusterSequence.hh 1881 2011-01-27 14:25:24Z soyez $
//
// 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


#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"
#include "fastjet/SharedPtr.hh"
#include "fastjet/internal/LimitedWarning.hh"


FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh

//----------------------------------------------------------------------
// here's where we put the main page for fastjet (as explained in the
// Doxygen faq)
// We put in inside te fastjet namespace to have the links without
// having to specify (fastjet::)
//......................................................................
//----------------------------------------------------------------------

// forward declaration
//class ClusterSequence;

class ClusterSequence {


 public: 

  ClusterSequence () {}

  template<class L> ClusterSequence (const std::vector<L> & pseudojets, 
                   const double & R = 1.0,
                   const Strategy & strategy = Best,
                   const bool & writeout_combinations = false);


  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).

  std::vector<PseudoJet> inclusive_jets (const double & ptmin = 0.0) const;

  int n_exclusive_jets (const double & dcut) const;

  std::vector<PseudoJet> exclusive_jets (const double & dcut) const;

  std::vector<PseudoJet> exclusive_jets (const int & njets) const;

  double exclusive_dmerge (const int & njets) const;

  double exclusive_dmerge_max (const int & njets) const;

  double exclusive_ymerge (int njets) const {return exclusive_dmerge(njets) / Q2();}

  double exclusive_ymerge_max (int njets) const {return exclusive_dmerge_max(njets)/Q2();}

  int n_exclusive_jets_ycut (double ycut) const {return n_exclusive_jets(ycut*Q2());}

  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;

  std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet, 
                                            const double & dcut) const;

  int n_exclusive_subjets(const PseudoJet & jet, 
                          const double & dcut) const;

  std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet, 
                                            int nsub) const;

  double exclusive_subdmerge(const PseudoJet & jet, int nsub) const;

  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;

  double Q() const {return _Qtot;}
  double Q2() const {return _Qtot*_Qtot;}

  bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;

  bool has_parents(const PseudoJet & jet, PseudoJet & parent1, 
               PseudoJet & parent2) const;

  bool has_child(const PseudoJet & jet, PseudoJet & child) const;

  bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;

  bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;

  
  std::vector<PseudoJet> constituents (const PseudoJet & jet) const;


  void print_jets_for_root(const std::vector<PseudoJet> & jets, 
                           std::ostream & ostr = std::cout) const;

  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;

  void add_constituents (const PseudoJet & jet, 
                         std::vector<PseudoJet> & subjet_vector) const;

  inline Strategy strategy_used () const {return _strategy;}

  std::string strategy_string () const {return strategy_string(_strategy);}

  std::string strategy_string (Strategy strategy_in) const;


  const JetDefinition & jet_def() const {return _jet_def;}

  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

  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);
  }

  void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, 
                                      const PseudoJet & newjet, 
                                      int & newjet_k);

  void plugin_record_iB_recombination(int jet_i, double diB) {
    assert(plugin_activated());
    _do_iB_recombination_step(jet_i, diB);
  }

  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";}
  };

  inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) {
    _extras = extras_in;
  }

  inline bool plugin_activated() const {return _plugin_activated;}

  const Extras * extras() const {return _extras.get();}

  template<class GBJ> void plugin_simple_N2_cluster () {
    assert(plugin_activated());
    _simple_N2_cluster<GBJ>();
  }


public:
  static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
  static void set_jet_finder (JetAlgorithm jet_algorithm)    {_default_jet_algorithm = jet_algorithm;}


  struct history_element{
    int parent1; 



    int parent2; 





    int child;   




    int jetp_index; 






    double dij;  


    double max_dij_so_far; 

  };

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

  const std::vector<PseudoJet> & jets()    const;

  const std::vector<history_element> & history() const;

  unsigned int n_particles() const;

  std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;

  std::vector<int> unique_history_order() const;

  std::vector<PseudoJet> unclustered_particles() const;

  void transfer_from_sequence(ClusterSequence & from_seq, bool transfer_ownership=true);


protected:
  static JetAlgorithm _default_jet_algorithm;
  JetDefinition _jet_def;

  bool _potentially_valid(const PseudoJet & jet) const {
    return jet.cluster_hist_index() >= 0 
      && jet.cluster_hist_index() < int(_history.size());
  }

  template<class L> void _transfer_input_jets(
                                     const std::vector<L> & pseudojets);

  void _initialise_and_run (const JetDefinition & jet_def,
                            const bool & writeout_combinations);

  void _initialise_and_run (const double & R,
                            const Strategy & strategy,
                            const bool & writeout_combinations);

  void _decant_options(const JetDefinition & jet_def,
                       const bool & writeout_combinations);

  void _fill_initial_history();

  void _do_ij_recombination_step(const int & jet_i, const int & jet_j, 
                                 const double & dij, int & newjet_k);

  void _do_iB_recombination_step(const int & jet_i, const double & diB);


  std::vector<PseudoJet> _jets;


  std::vector<history_element> _history;

  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 _Qtot;
  Strategy    _strategy;
  JetAlgorithm  _jet_algorithm;

  SharedPtr<ClusterSequenceWrapper> _wrapper_to_this;

 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();

  // NSqrtN method for C/A
  void _fast_NsqrtN_cluster();

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

  void _extract_tree_children(int pos, std::valarray<bool> &, 
                const std::valarray<int> &, std::vector<int> &) const;

  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;

  void _add_ktdistance_to_map(const int & ii, 
                              DistMap & DijMap,
                              const DynamicNearestNeighbours * DNN);

  void _print_banner();
  static bool _first_time;

  static int _n_exclusive_warnings;

  static LimitedWarning _changed_strategy_warning;

  //----------------------------------------------------------------------
  struct BriefJet {
    double     eta, phi, kt2, NN_dist;
    BriefJet * NN;
    int        _jets_index;
  };


  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

  template <class J> void _bj_set_jetinfo( J * const jet, 
                                                 const int _jets_index) const;

  void _bj_remove_from_tiles( TiledJet * const jet) const;

  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;

  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...

  template <class J> void _bj_set_NN_nocross(J * const jeta, 
            J * const head, const J * const tail) const;

  template <class J> void _bj_set_NN_crosscheck(J * const jeta, 
            J * const head, const J * const tail) const;
  


  static const int n_tile_neighbours = 9;
  //----------------------------------------------------------------------
  struct Tile {
    Tile *   begin_tiles[n_tile_neighbours]; 
    Tile **  surrounding_tiles; 
    Tile **  RH_tiles;  
    Tile **  end_tiles; 
    TiledJet * head;    
    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);


  //----------------------------------------------------------------------
  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
  };

  //void _dummy_N2_cluster_instantiation();


  void _simple_N2_cluster_BriefJet();
  void _simple_N2_cluster_EEBriefJet();
};


//**********************************************************************
//**************    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);
}


//----------------------------------------------------------------------
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_END_NAMESPACE

#endif // __FASTJET_CLUSTERSEQUENCE_HH__