PseudoJet.hh

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: PseudoJet.hh 958 2007-11-28 17:43:48Z cacciari $
00003 //
00004 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet; if not, write to the Free Software
00026 //  Foundation, Inc.:
00027 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00028 //----------------------------------------------------------------------
00029 //ENDHEADER
00030 
00031 
00032 #ifndef __FASTJET_PSEUDOJET_HH__
00033 #define __FASTJET_PSEUDOJET_HH__
00034 
00035 #include<valarray>
00036 #include<vector>
00037 #include<cassert>
00038 #include<cmath>
00039 #include<iostream>
00040 #include "fastjet/internal/numconsts.hh"
00041 
00042 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00043 
00044 //using namespace std;
00045 
00048 const double MaxRap = 1e5;
00049 
00052 class PseudoJet {
00053 
00054  public:
00055   PseudoJet() {};
00057   PseudoJet(const double px, const double py, const double pz, const double E);
00059   template <class L> PseudoJet(const L & some_four_vector) ;
00060 
00061   // first "const double &" says that result is a reference to the
00062   // stored value and that we will not change that stored value.
00063   //
00064   // second "const" says that "this" will not be modified by these
00065   // functions.
00066   inline double E()   const {return _E;}
00067   inline double e()   const {return _E;} // like CLHEP
00068   inline double px()  const {return _px;}
00069   inline double py()  const {return _py;}
00070   inline double pz()  const {return _pz;}
00071 
00073   inline const double phi() const {return phi_02pi();}
00074 
00076   inline const double phi_std()  const {
00077     return _phi > pi ? _phi-twopi : _phi;}
00078 
00080   inline const double phi_02pi() const {return _phi;}
00081 
00084   inline double rap() const {return _rap;}
00085 
00087   inline double rapidity() const {return _rap;} // like CLHEP
00088 
00091   double pseudorapidity() const;
00092   double eta() const {return pseudorapidity();}
00093 
00095   inline double kt2() const {return _kt2;}
00097   inline double perp2() const {return _kt2;}  // like CLHEP
00099   inline double  perp() const {return sqrt(_kt2);}    // like CLHEP
00101   inline double  m2() const {return (_E+_pz)*(_E-_pz)-_kt2;}    
00103   inline double mperp2() const {return (_E+_pz)*(_E-_pz);}
00105   inline double mperp() const {return sqrt(std::abs(mperp2()));}
00108   inline double  m() const;    
00110   double operator () (int i) const ; 
00112   inline double operator [] (int i) const { return (*this)(i); }; // this too
00113 
00114   // taken from CLHEP
00115   enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
00116 
00117 
00120   PseudoJet & boost(const PseudoJet & prest);
00123   PseudoJet & unboost(const PseudoJet & prest);
00124 
00127   inline const int & cluster_hist_index() const {return _cluster_hist_index;}
00129   inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;}
00130 
00132   inline const int cluster_sequence_history_index() const {
00133     return cluster_hist_index();}
00136   inline void set_cluster_sequence_history_index(const int index) {
00137     set_cluster_hist_index(index);}
00138 
00139 
00141   inline const int & user_index() const {return _user_index;}
00143   inline void set_user_index(const int index) {_user_index = index;}
00144 
00147   std::valarray<double> four_mom() const;
00148 
00150   double kt_distance(const PseudoJet & other) const;
00151 
00153   double plain_distance(const PseudoJet & other) const;
00156   inline double squared_distance(const PseudoJet & other) const {
00157     return plain_distance(other);}
00158 
00160   //friend inline double 
00161   //  kt_distance(const PseudoJet & jet1, const PseudoJet & jet2) { 
00162   //                                      return jet1.kt_distance(jet2);}
00163 
00165   inline double beam_distance() const {return _kt2;}
00166 
00167 
00168   void operator*=(double);
00169   void operator/=(double);
00170   void operator+=(const PseudoJet &);
00171   void operator-=(const PseudoJet &);
00172 
00175   inline void reset(double px, double py, double pz, double E);
00176   
00181   inline void reset(const PseudoJet & psjet) {
00182     (*this) = psjet;
00183   }
00184 
00188   template <class L> inline void reset(const L & some_four_vector) {
00189     reset(some_four_vector[0], some_four_vector[1],
00190           some_four_vector[2], some_four_vector[3]);
00191   }
00192 
00193  private: 
00194   // NB: following order must be kept for things to behave sensibly...
00195   double _px,_py,_pz,_E;
00196   double _phi, _rap, _kt2; 
00197   int    _cluster_hist_index, _user_index;
00199   void _finish_init();
00201   void _reset_indices();
00202 
00203   //vertex_type * vertex0, vertex1;
00204 };
00205 
00206 
00207 //----------------------------------------------------------------------
00208 // routines for basic binary operations
00209 
00210 PseudoJet operator+(const PseudoJet &, const PseudoJet &);
00211 PseudoJet operator-(const PseudoJet &, const PseudoJet &);
00212 PseudoJet operator*(double, const PseudoJet &);
00213 PseudoJet operator*(const PseudoJet &, double);
00214 PseudoJet operator/(const PseudoJet &, double);
00215 
00217 bool have_same_momentum(const PseudoJet &, const PseudoJet &);
00218 
00220 PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0);
00221 
00222 //----------------------------------------------------------------------
00223 // Routines to do with providing sorted arrays of vectors.
00224 
00226 std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets);
00227 
00229 std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets);
00230 
00232 std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets);
00233 
00234 //----------------------------------------------------------------------
00235 // some code to help sorting
00236 
00239 void sort_indices(std::vector<int> & indices, 
00240                   const std::vector<double> & values);
00241 
00246 template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects, 
00247                                               const std::vector<double> & values);
00248 
00250 class IndexedSortHelper {
00251 public:
00252   inline IndexedSortHelper (const std::vector<double> * reference_values) {
00253     _ref_values = reference_values;
00254   };
00255   inline int operator() (const int & i1, const int & i2) const {
00256     return  (*_ref_values)[i1] < (*_ref_values)[i2];
00257   };
00258 private:
00259   const std::vector<double> * _ref_values;
00260 };
00261 
00262 
00263 //----------------------------------------------------------------------
00265 // NB: do not know if it really needs to be inline, but when it wasn't
00266 //     linking failed with g++ (who knows what was wrong...)
00267 template <class L> inline  PseudoJet::PseudoJet(const L & some_four_vector) {
00268 
00269   _px = some_four_vector[0];
00270   _py = some_four_vector[1];
00271   _pz = some_four_vector[2];
00272   _E  = some_four_vector[3];
00273   _finish_init();
00274   // some default values for these two indices
00275   _reset_indices();
00276 }
00277 
00278 
00279 //----------------------------------------------------------------------
00280 inline void PseudoJet::_reset_indices() { 
00281   set_cluster_hist_index(-1);
00282   set_user_index(-1);
00283 }
00284 
00285 //----------------------------------------------------------------------
00289 // template<> inline void PseudoJet::reset<PseudoJet>(const PseudoJet & psjet) {
00290 //   (*this) = psjet;
00291 // }
00292 
00300 
00301 // taken literally from CLHEP
00302 inline double PseudoJet::m() const {
00303   double mm = m2();
00304   return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
00305 }
00306 
00307 
00308 inline void PseudoJet::reset(double px, double py, double pz, double E) {
00309   _px = px;
00310   _py = py;
00311   _pz = pz;
00312   _E  = E;
00313   _finish_init();
00314   _reset_indices();
00315 }
00316 
00317 
00318 FASTJET_END_NAMESPACE
00319 
00320 #endif // __FASTJET_PSEUDOJET_HH__

Generated on Tue Dec 18 17:05:02 2007 for fastjet by  doxygen 1.5.2