30 #ifndef __FASTJET_PSEUDOJET_HH__
31 #define __FASTJET_PSEUDOJET_HH__
38 #include "fastjet/internal/numconsts.hh"
39 #include "fastjet/internal/IsBase.hh"
40 #include "fastjet/SharedPtr.hh"
41 #include "fastjet/Error.hh"
42 #include "fastjet/PseudoJetStructureBase.hh"
44 FASTJET_BEGIN_NAMESPACE
74 PseudoJet() : _px(0), _py(0), _pz(0), _E(0) {_finish_init(); _reset_indices();}
76 PseudoJet(
const double px,
const double py,
const double pz,
const double E);
79 template <
class L>
PseudoJet(
const L & some_four_vector);
95 inline double E()
const {
return _E;}
96 inline double e()
const {
return _E;}
97 inline double px()
const {
return _px;}
98 inline double py()
const {
return _py;}
99 inline double pz()
const {
return _pz;}
102 inline double phi()
const {
return phi_02pi();}
106 _ensure_valid_rap_phi();
107 return _phi > pi ? _phi-twopi : _phi;}
111 _ensure_valid_rap_phi();
117 inline double rap()
const {
118 _ensure_valid_rap_phi();
127 double pseudorapidity()
const;
128 double eta()
const {
return pseudorapidity();}
131 inline double pt2()
const {
return _kt2;}
133 inline double pt()
const {
return sqrt(_kt2);}
135 inline double perp2()
const {
return _kt2;}
137 inline double perp()
const {
return sqrt(_kt2);}
139 inline double kt2()
const {
return _kt2;}
142 inline double m2()
const {
return (_E+_pz)*(_E-_pz)-_kt2;}
145 inline double m()
const;
148 inline double mperp2()
const {
return (_E+_pz)*(_E-_pz);}
150 inline double mperp()
const {
return sqrt(std::abs(mperp2()));}
152 inline double mt2()
const {
return (_E+_pz)*(_E-_pz);}
154 inline double mt()
const {
return sqrt(std::abs(mperp2()));}
157 inline double modp2()
const {
return _kt2+_pz*_pz;}
159 inline double modp()
const {
return sqrt(_kt2+_pz*_pz);}
162 inline double Et()
const {
return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);}
164 inline double Et2()
const {
return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);}
167 double operator () (
int i)
const ;
169 inline double operator [] (
int i)
const {
return (*
this)(i); };
174 double kt_distance(
const PseudoJet & other)
const;
177 double plain_distance(
const PseudoJet & other)
const;
181 return plain_distance(other);}
186 return sqrt(squared_distance(other));
191 double delta_phi_to(
const PseudoJet & other)
const;
203 std::valarray<double> four_mom()
const;
208 enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
217 PseudoJet & boost(
const PseudoJet & prest);
220 PseudoJet & unboost(
const PseudoJet & prest);
222 void operator*=(
double);
223 void operator/=(
double);
224 void operator+=(
const PseudoJet &);
225 void operator-=(
const PseudoJet &);
229 inline void reset(
double px,
double py,
double pz,
double E);
245 template <
class L>
inline void reset(
const L & some_four_vector) {
253 const PseudoJet * pj = cast_if_derived<const PseudoJet>(&some_four_vector);
258 reset(some_four_vector[0], some_four_vector[1],
259 some_four_vector[2], some_four_vector[3]);
266 inline void reset_PtYPhiM(
double pt_in,
double y_in,
double phi_in,
double m_in=0.0) {
267 reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in);
274 inline void reset_momentum(
double px,
double py,
double pz,
double E);
279 inline void reset_momentum(
const PseudoJet & pj);
283 void reset_momentum_PtYPhiM(
double pt,
double y,
double phi,
double m=0.0);
289 reset_momentum(some_four_vector[0], some_four_vector[1],
290 some_four_vector[2], some_four_vector[3]);
305 void set_cached_rap_phi(
double rap,
double phi);
378 _user_info.reset(user_info_in);
412 return dynamic_cast<const L &
>(* _user_info.get());
417 return _user_info.get();
424 return _user_info.get() &&
dynamic_cast<const L *
>(_user_info.get());
429 if (!_user_info())
return NULL;
430 return _user_info.get();
468 std::string description()
const;
482 bool has_associated_cluster_sequence()
const;
488 bool has_valid_cluster_sequence()
const;
496 const ClusterSequence* associated_cs()
const {
return associated_cluster_sequence();}
501 return validated_cs();
509 return validated_csab();
530 bool has_structure()
const;
564 template<
typename StructureType>
565 const StructureType & structure()
const;
570 template<
typename TransformerType>
571 bool has_structure_of()
const;
578 template<
typename TransformerType>
579 const typename TransformerType::StructureType & structure_of()
const;
598 virtual bool has_partner(
PseudoJet &partner)
const;
606 virtual bool has_child(
PseudoJet &child)
const;
621 virtual bool contains(
const PseudoJet &constituent)
const;
628 virtual bool is_inside(
const PseudoJet &jet)
const;
632 virtual bool has_constituents()
const;
638 virtual std::vector<PseudoJet> constituents()
const;
642 virtual bool has_exclusive_subjets()
const;
655 std::vector<PseudoJet> exclusive_subjets (
const double & dcut)
const;
663 int n_exclusive_subjets(
const double & dcut)
const;
673 std::vector<PseudoJet> exclusive_subjets (
int nsub)
const;
683 std::vector<PseudoJet> exclusive_subjets_up_to (
int nsub)
const;
690 double exclusive_subdmerge(
int nsub)
const;
698 double exclusive_subdmerge_max(
int nsub)
const;
708 virtual bool has_pieces()
const;
714 virtual std::vector<PseudoJet> pieces()
const;
722 virtual bool has_area()
const;
726 virtual double area()
const;
731 virtual double area_error()
const;
739 virtual bool is_pure_ghost()
const;
757 return cluster_hist_index();}
761 set_cluster_hist_index(index);}
773 double _px,_py,_pz,_E;
774 mutable double _phi, _rap;
776 int _cluster_hist_index, _user_index;
781 void _reset_indices();
785 inline void _ensure_valid_rap_phi()
const {
786 if (_phi == pseudojet_invalid_phi) _set_rap_phi();
790 void _set_rap_phi()
const;
797 PseudoJet operator+(
const PseudoJet &,
const PseudoJet &);
798 PseudoJet operator-(
const PseudoJet &,
const PseudoJet &);
799 PseudoJet
operator*(
double,
const PseudoJet &);
800 PseudoJet
operator*(
const PseudoJet &,
double);
801 PseudoJet operator/(
const PseudoJet &,
double);
806 bool operator==(
const PseudoJet &,
const PseudoJet &);
813 bool operator==(
const PseudoJet & jet,
const double val);
819 inline double dot_product(
const PseudoJet & a,
const PseudoJet & b) {
820 return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
828 PseudoJet
PtYPhiM(
double pt,
double y,
double phi,
double m = 0.0);
834 std::vector<PseudoJet>
sorted_by_pt(
const std::vector<PseudoJet> & jets);
840 std::vector<PseudoJet>
sorted_by_E(
const std::vector<PseudoJet> & jets);
843 std::vector<PseudoJet>
sorted_by_pz(
const std::vector<PseudoJet> & jets);
850 void sort_indices(std::vector<int> & indices,
851 const std::vector<double> & values);
858 const std::vector<double> & values);
865 class IndexedSortHelper {
867 inline IndexedSortHelper (
const std::vector<double> * reference_values) {
868 _ref_values = reference_values;
870 inline int operator() (
const int & i1,
const int & i2)
const {
871 return (*_ref_values)[i1] < (*_ref_values)[i2];
874 const std::vector<double> * _ref_values;
882 template <
class L>
inline PseudoJet::PseudoJet(
const L & some_four_vector) {
883 reset(some_four_vector);
887 inline void PseudoJet::_reset_indices() {
888 set_cluster_hist_index(-1);
896 inline double PseudoJet::m()
const {
898 return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
902 inline void PseudoJet::reset(
double px_in,
double py_in,
double pz_in,
double E_in) {
911 inline void PseudoJet::reset_momentum(
double px_in,
double py_in,
double pz_in,
double E_in) {
919 inline void PseudoJet::reset_momentum(
const PseudoJet & pj) {
938 template<
typename StructureType>
939 const StructureType & PseudoJet::structure()
const{
940 return dynamic_cast<const StructureType &
>(* validated_structure_ptr());
946 template<
typename TransformerType>
947 bool PseudoJet::has_structure_of()
const{
948 if (!_structure())
return false;
950 return dynamic_cast<const typename TransformerType::StructureType *
>(_structure.get()) != 0;
956 template<
typename TransformerType>
957 const typename TransformerType::StructureType & PseudoJet::structure_of()
const{
959 throw Error(
"Trying to access the structure of a PseudoJet without an associated structure");
961 return dynamic_cast<const typename TransformerType::StructureType &
>(*_structure);
978 PseudoJet join(
const std::vector<PseudoJet> & pieces);
994 FASTJET_END_NAMESPACE
996 #endif // __FASTJET_PSEUDOJET_HH__