00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
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
00043
00044
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
00062
00063
00064
00065
00066 inline double E() const {return _E;}
00067 inline double e() const {return _E;}
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;}
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;}
00099 inline double perp() const {return sqrt(_kt2);}
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); };
00113
00114
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
00161
00162
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
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
00204 };
00205
00206
00207
00208
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
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
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
00266
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
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
00290
00291
00292
00300
00301
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__