PseudoJet.hh

Go to the documentation of this file.

//STARTHEADER
// $Id: PseudoJet.hh 492 2007-03-02 08:12:18Z 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


// NB: to do at some point
//     - add += and *= operators

#ifndef __FASTJET_PSEUDOJET_HH__
#define __FASTJET_PSEUDOJET_HH__

#include<valarray>
#include<vector>
#include<cassert>
#include <cmath>
#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 const double phi() const {return phi_02pi();};

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

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

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

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

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


  inline const int & cluster_hist_index() const {return _cluster_hist_index;};
  inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;};

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

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

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

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

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

//----------------------------------------------------------------------
// 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];
  this->_finish_init();
  // some default values for these two indices
  set_cluster_hist_index(-1);
  set_user_index(-1);
}



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


FASTJET_END_NAMESPACE

#endif // __FASTJET_PSEUDOJET_HH__

---