PseudoJet.hh

Go to the documentation of this file.

//STARTHEADER
// $Id: PseudoJet.hh 1510 2009-04-13 08:48:41Z 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"

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

//using namespace std;

const double MaxRap = 1e5;

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

  // first "const double &" says that result is a reference to the
  // stored value and that we will not change that stored value.
  //
  // second "const" says that "this" will not be modified by these
  // functions.
  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 {
    return _phi > pi ? _phi-twopi : _phi;}

  inline double phi_02pi() const {return _phi;}

  inline double rap() const {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


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

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


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

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

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


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

 private: 
  // NB: following order must be kept for things to behave sensibly...
  double _px,_py,_pz,_E;
  double _phi, _rap, _kt2; 
  int    _cluster_hist_index, _user_index;
  void _finish_init();
  void _reset_indices();

  //vertex_type * vertex0, vertex1;
};


//----------------------------------------------------------------------
// 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) {

  _px = some_four_vector[0];
  _py = some_four_vector[1];
  _pz = some_four_vector[2];
  _E  = some_four_vector[3];
  _finish_init();
  // some default values for these two indices
  _reset_indices();
}


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

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