PseudoJet.hh

//STARTHEADER
// $Id: PseudoJet.hh 1874 2011-01-27 11:37:51Z 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


#ifndef __FASTJET_PSEUDOJET_HH__
#define __FASTJET_PSEUDOJET_HH__

#include<valarray>
#include<vector>
#include<cassert>
#include<cmath>
#include<iostream>
#include "fastjet/internal/numconsts.hh"
#include "fastjet/internal/IsBase.hh"
#include "fastjet/internal/DerivedPseudoJetHelper.hh"
#include "fastjet/SharedPtr.hh"
#include "fastjet/Error.hh"
#include "fastjet/ClusterSequenceWrapper.hh"

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

//using namespace std;

const double MaxRap = 1e5;

const double pseudojet_invalid_phi = -100.0;

// forward definition
class ClusterSequenceAreaBase;

class PseudoJet {

 public:
  //----------------------------------------------------------------------
  //\{
  PseudoJet() {};
  PseudoJet(const double px, const double py, const double pz, const double E);
  template <class L> PseudoJet(const L & some_four_vector) ;

  virtual ~PseudoJet(){};
  //\} ---- end of constructors and destructors --------------------------

  //----------------------------------------------------------------------
  //\{
  //----------------------------------------------------------------------
  inline double E()   const {return _E;}
  inline double e()   const {return _E;} // like CLHEP
  inline double px()  const {return _px;}
  inline double py()  const {return _py;}
  inline double pz()  const {return _pz;}

  inline double phi() const {return phi_02pi();}

  inline double phi_std()  const {
    _ensure_valid_rap_phi();
    return _phi > pi ? _phi-twopi : _phi;}

  inline double phi_02pi() const {
    _ensure_valid_rap_phi();
    return _phi;
  }

  inline double rap() const {
    _ensure_valid_rap_phi();
    return _rap;
  }

  inline double rapidity() const {return rap();} // like CLHEP

  double pseudorapidity() const;
  double eta() const {return pseudorapidity();}

  inline double kt2() const {return _kt2;}
  inline double perp2() const {return _kt2;}  // like CLHEP
  inline double  perp() const {return sqrt(_kt2);}    // like CLHEP
  inline double  m2() const {return (_E+_pz)*(_E-_pz)-_kt2;}    
  inline double mperp2() const {return (_E+_pz)*(_E-_pz);}
  inline double mperp() const {return sqrt(std::abs(mperp2()));}
  inline double  m() const;    
  inline double modp2() const {return _kt2+_pz*_pz;}
  inline double Et() const {return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);}
  inline double Et2() const {return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);}

  double operator () (int i) const ; 
  inline double operator [] (int i) const { return (*this)(i); }; // this too



  double kt_distance(const PseudoJet & other) const;

  double plain_distance(const PseudoJet & other) const;
  inline double squared_distance(const PseudoJet & other) const {
    return plain_distance(other);}

  double delta_phi_to(const PseudoJet & other) const;

  //friend inline double 
  //  kt_distance(const PseudoJet & jet1, const PseudoJet & jet2) { 
  //                                      return jet1.kt_distance(jet2);}

  inline double beam_distance() const {return _kt2;}

  std::valarray<double> four_mom() const;

  //\}  ------- end of kinematic access functions

  // taken from CLHEP
  enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };


  //----------------------------------------------------------------------
  //\{
  //----------------------------------------------------------------------
  PseudoJet & boost(const PseudoJet & prest);
  PseudoJet & unboost(const PseudoJet & prest);

  void operator*=(double);
  void operator/=(double);
  void operator+=(const PseudoJet &);
  void operator-=(const PseudoJet &);

  inline void reset(double px, double py, double pz, double E);
  
  inline void reset(const PseudoJet & psjet) {
    (*this) = psjet;
  }

  template <class L> inline void reset(const L & some_four_vector) {
    reset(some_four_vector[0], some_four_vector[1],
          some_four_vector[2], some_four_vector[3]);
  }

  //\} --- end of kin mod functions ------------------------------------

  //----------------------------------------------------------------------
  //
  //\{

  inline int user_index() const {return _user_index;}
  inline void set_user_index(const int index) {_user_index = index;}

  //\} ----- end of use index functions ---------------------------------

  //----------------------------------------------------------------------
  //\{

  class ExtraInfo{
  public:
    // dummy ctor
    ExtraInfo(){};

    // dummy virtual dtor
    // makes it polymorphic to allow for dynamic_cast
    virtual ~ExtraInfo(){}; 
  };

  class InexistentExtraInfo : public Error {
  public:
    InexistentExtraInfo();
  };

  const ExtraInfo* extra_info() const{
    if (!_extra_info()) return NULL;
    return _extra_info.get();
  }

  void set_extra_info(ExtraInfo * extra_info_in) {
    _extra_info.reset(extra_info_in);
  }

  template<class L>
  const L & extra_info_cast() const{
    if (_extra_info.get() == 0) throw InexistentExtraInfo();
    return dynamic_cast<const L &>(* _extra_info.get());
  }

  SharedPtr<ExtraInfo> & extra_info_shared(){
    return _extra_info;
  }

  const SharedPtr<ExtraInfo> & extra_info_shared() const{
    return _extra_info;
  }

  // \} --- end of extra info functions ---------------------------------

  //-------------------------------------------------------------
  //\{
  //-------------------------------------------------------------
  bool has_associated_cluster_sequence() const;

  const ClusterSequence* associated_cluster_sequence() const;
  //\}

  //-------------------------------------------------------------
  //-------------------------------------------------------------
  //\{

  virtual bool has_partner(PseudoJet &partner) const;

  virtual bool has_child(PseudoJet &child) const;

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

  virtual bool contains(const PseudoJet &constituent) const;

  virtual bool is_inside(const PseudoJet &jet) const;

  virtual std::vector<PseudoJet> constituents() const;

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

  int n_exclusive_subjets(const double & dcut) const;

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

  double exclusive_subdmerge(int nsub) const;

  double exclusive_subdmerge_max(int nsub) const;


  // the following ones require a computation of the area in the
  // parent ClusterSequence (See ClusterSequenceAreaBase for details)
  //------------------------------------------------------------------

  virtual bool has_area() const;

  virtual double area() const;

  virtual double area_error() const;

  virtual PseudoJet area_4vector() const;

  virtual bool is_pure_ghost() const;

  //\} --- end of jet structure -------------------------------------



  //----------------------------------------------------------------------
  //----------------------------------------------------------------------
  //\{
  inline int cluster_hist_index() const {return _cluster_hist_index;}
  inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;}

  inline int cluster_sequence_history_index() const {
    return cluster_hist_index();}
  inline void set_cluster_sequence_history_index(const int index) {
    set_cluster_hist_index(index);}

  void set_associated_csw(const SharedPtr<ClusterSequenceWrapper> &csw){
    _associated_csw = csw;
  }

  const SharedPtr<ClusterSequenceWrapper> & associated_csw() const {
    return _associated_csw;
  }
  
  const ClusterSequence * validated_cs() const;

  const ClusterSequenceAreaBase * validated_csab() const;

  //\} ---- end of internal use functions ---------------------------
  
 private: 
  // NB: following order must be kept for things to behave sensibly...
  double _px,_py,_pz,_E;
  mutable double _phi, _rap;
  double _kt2; 
  int    _cluster_hist_index, _user_index;

  SharedPtr<ClusterSequenceWrapper> _associated_csw;
  SharedPtr<ExtraInfo> _extra_info;

  void _finish_init();
  void _reset_indices();

  inline void _ensure_valid_rap_phi() const {
    if (_phi == pseudojet_invalid_phi) _set_rap_phi();
  }

  void _set_rap_phi() const;
};


//----------------------------------------------------------------------
// routines for basic binary operations

PseudoJet operator+(const PseudoJet &, const PseudoJet &);
PseudoJet operator-(const PseudoJet &, const PseudoJet &);
PseudoJet operator*(double, const PseudoJet &);
PseudoJet operator*(const PseudoJet &, double);
PseudoJet operator/(const PseudoJet &, double);

inline double dot_product(const PseudoJet & a, const PseudoJet & b) {
  return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
}

bool have_same_momentum(const PseudoJet &, const PseudoJet &);

PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0);

//----------------------------------------------------------------------
// Routines to do with providing sorted arrays of vectors.

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

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

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

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

//----------------------------------------------------------------------
// some code to help sorting

void sort_indices(std::vector<int> & indices, 
                  const std::vector<double> & values);

template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects, 
                                              const std::vector<double> & values);

class IndexedSortHelper {
public:
  inline IndexedSortHelper (const std::vector<double> * reference_values) {
    _ref_values = reference_values;
  };
  inline int operator() (const int & i1, const int & i2) const {
    return  (*_ref_values)[i1] < (*_ref_values)[i2];
  };
private:
  const std::vector<double> * _ref_values;
};


//----------------------------------------------------------------------
// NB: do not know if it really needs to be inline, but when it wasn't
//     linking failed with g++ (who knows what was wrong...)
template <class L> inline  PseudoJet::PseudoJet(const L & some_four_vector) {
  // now check whether L is simply a class that implements
  // some_fuor_vector[0--3] or actually is derived from PseudoJet and
  // has extra information
  DerivedPseudoJetHelper<L, IsBaseAndDerived<PseudoJet,L>::value> dpj_helper(some_four_vector);

  if (dpj_helper() != NULL){
    const PseudoJet *pj = dpj_helper();
    reset(*pj);
  } else {
    reset(some_four_vector);
  }
}


//----------------------------------------------------------------------
inline void PseudoJet::_reset_indices() { 
  set_cluster_hist_index(-1);
  set_user_index(-1);
  _associated_csw.reset();
}

//----------------------------------------------------------------------
// template<> inline void PseudoJet::reset<PseudoJet>(const PseudoJet & psjet) {
//   (*this) = psjet;
// }


// taken literally from CLHEP
inline double PseudoJet::m() const {
  double mm = m2();
  return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
}


inline void PseudoJet::reset(double px, double py, double pz, double E) {
  _px = px;
  _py = py;
  _pz = pz;
  _E  = E;
  _finish_init();
  _reset_indices();
}


FASTJET_END_NAMESPACE

#endif // __FASTJET_PSEUDOJET_HH__