32 #ifndef __FASTJET_PSEUDOJET_HH__
33 #define __FASTJET_PSEUDOJET_HH__
40 #include "fastjet/internal/numconsts.hh"
41 #include "fastjet/internal/IsBase.hh"
42 #include "fastjet/SharedPtr.hh"
43 #include "fastjet/Error.hh"
44 #include "fastjet/PseudoJetStructureBase.hh"
46 FASTJET_BEGIN_NAMESPACE
56 const double pseudojet_invalid_rap = -1e200;
79 PseudoJet() : _px(0), _py(0), _pz(0), _E(0) {_finish_init(); _reset_indices();}
82 PseudoJet(
const double px,
const double py,
const double pz,
const double E);
85 template <
class L>
PseudoJet(
const L & some_four_vector);
101 inline double E()
const {
return _E;}
102 inline double e()
const {
return _E;}
103 inline double px()
const {
return _px;}
104 inline double py()
const {
return _py;}
105 inline double pz()
const {
return _pz;}
108 inline double phi()
const {
return phi_02pi();}
112 _ensure_valid_rap_phi();
113 return _phi > pi ? _phi-twopi : _phi;}
117 _ensure_valid_rap_phi();
123 inline double rap()
const {
124 _ensure_valid_rap_phi();
133 double pseudorapidity()
const;
134 double eta()
const {
return pseudorapidity();}
137 inline double pt2()
const {
return _kt2;}
139 inline double pt()
const {
return sqrt(_kt2);}
141 inline double perp2()
const {
return _kt2;}
143 inline double perp()
const {
return sqrt(_kt2);}
145 inline double kt2()
const {
return _kt2;}
148 inline double m2()
const {
return (_E+_pz)*(_E-_pz)-_kt2;}
151 inline double m()
const;
154 inline double mperp2()
const {
return (_E+_pz)*(_E-_pz);}
156 inline double mperp()
const {
return sqrt(std::abs(mperp2()));}
158 inline double mt2()
const {
return (_E+_pz)*(_E-_pz);}
160 inline double mt()
const {
return sqrt(std::abs(mperp2()));}
163 inline double modp2()
const {
return _kt2+_pz*_pz;}
165 inline double modp()
const {
return sqrt(_kt2+_pz*_pz);}
168 inline double Et()
const {
return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);}
170 inline double Et2()
const {
return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);}
173 double operator () (
int i)
const ;
175 inline double operator [] (
int i)
const {
return (*
this)(i); };
180 double kt_distance(
const PseudoJet & other)
const;
183 double plain_distance(
const PseudoJet & other)
const;
187 return plain_distance(other);}
192 return sqrt(squared_distance(other));
197 double delta_phi_to(
const PseudoJet & other)
const;
209 std::valarray<double> four_mom()
const;
214 enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
223 PseudoJet & boost(
const PseudoJet & prest);
226 PseudoJet & unboost(
const PseudoJet & prest);
228 void operator*=(
double);
229 void operator/=(
double);
230 void operator+=(
const PseudoJet &);
231 void operator-=(
const PseudoJet &);
235 inline void reset(
double px,
double py,
double pz,
double E);
251 template <
class L>
inline void reset(
const L & some_four_vector) {
263 const PseudoJet * pj = fastjet::cast_if_derived<const PseudoJet>(&some_four_vector);
268 reset(some_four_vector[0], some_four_vector[1],
269 some_four_vector[2], some_four_vector[3]);
276 inline void reset_PtYPhiM(
double pt_in,
double y_in,
double phi_in,
double m_in=0.0) {
277 reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in);
284 inline void reset_momentum(
double px,
double py,
double pz,
double E);
289 inline void reset_momentum(
const PseudoJet & pj);
293 void reset_momentum_PtYPhiM(
double pt,
double y,
double phi,
double m=0.0);
299 reset_momentum(some_four_vector[0], some_four_vector[1],
300 some_four_vector[2], some_four_vector[3]);
315 void set_cached_rap_phi(
double rap,
double phi);
388 _user_info.reset(user_info_in);
422 return dynamic_cast<const L &
>(* _user_info.get());
427 return _user_info.get();
434 return _user_info.get() &&
dynamic_cast<const L *
>(_user_info.get());
439 if (!_user_info())
return NULL;
440 return _user_info.get();
478 std::string description()
const;
492 bool has_associated_cluster_sequence()
const;
498 bool has_valid_cluster_sequence()
const;
506 const ClusterSequence* associated_cs()
const {
return associated_cluster_sequence();}
511 return validated_cs();
520 return validated_csab();
543 bool has_structure()
const;
577 template<
typename StructureType>
578 const StructureType & structure()
const;
583 template<
typename TransformerType>
584 bool has_structure_of()
const;
591 template<
typename TransformerType>
592 const typename TransformerType::StructureType & structure_of()
const;
611 virtual bool has_partner(
PseudoJet &partner)
const;
619 virtual bool has_child(
PseudoJet &child)
const;
634 virtual bool contains(
const PseudoJet &constituent)
const;
641 virtual bool is_inside(
const PseudoJet &jet)
const;
645 virtual bool has_constituents()
const;
651 virtual std::vector<PseudoJet> constituents()
const;
655 virtual bool has_exclusive_subjets()
const;
668 std::vector<PseudoJet> exclusive_subjets (
const double dcut)
const;
676 int n_exclusive_subjets(
const double dcut)
const;
686 std::vector<PseudoJet> exclusive_subjets (
int nsub)
const;
696 std::vector<PseudoJet> exclusive_subjets_up_to (
int nsub)
const;
705 double exclusive_subdmerge(
int nsub)
const;
715 double exclusive_subdmerge_max(
int nsub)
const;
725 virtual bool has_pieces()
const;
731 virtual std::vector<PseudoJet> pieces()
const;
740 virtual bool has_area()
const;
744 virtual double area()
const;
749 virtual double area_error()
const;
757 virtual bool is_pure_ghost()
const;
776 return cluster_hist_index();}
780 set_cluster_hist_index(index);}
792 double _px,_py,_pz,_E;
793 mutable double _phi, _rap;
795 int _cluster_hist_index, _user_index;
800 void _reset_indices();
804 inline void _ensure_valid_rap_phi()
const {
805 if (_phi == pseudojet_invalid_phi) _set_rap_phi();
809 void _set_rap_phi()
const;
812 friend PseudoJet
operator*(
double,
const PseudoJet &);
819 PseudoJet operator+(
const PseudoJet &,
const PseudoJet &);
820 PseudoJet operator-(
const PseudoJet &,
const PseudoJet &);
821 PseudoJet
operator*(
double,
const PseudoJet &);
822 PseudoJet
operator*(
const PseudoJet &,
double);
823 PseudoJet operator/(
const PseudoJet &,
double);
828 bool operator==(
const PseudoJet &,
const PseudoJet &);
835 bool operator==(
const PseudoJet & jet,
const double val);
836 inline bool operator==(
const double val,
const PseudoJet & jet) {
return jet == val;}
841 inline bool operator!=(
const double val,
const PseudoJet & a) {
return !(a==val);}
843 inline double dot_product(
const PseudoJet & a,
const PseudoJet & b) {
844 return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
852 PseudoJet
PtYPhiM(
double pt,
double y,
double phi,
double m = 0.0);
858 std::vector<PseudoJet>
sorted_by_pt(
const std::vector<PseudoJet> & jets);
864 std::vector<PseudoJet>
sorted_by_E(
const std::vector<PseudoJet> & jets);
867 std::vector<PseudoJet>
sorted_by_pz(
const std::vector<PseudoJet> & jets);
874 void sort_indices(std::vector<int> & indices,
875 const std::vector<double> & values);
882 const std::vector<double> & values);
889 class IndexedSortHelper {
891 inline IndexedSortHelper (
const std::vector<double> * reference_values) {
892 _ref_values = reference_values;
894 inline int operator() (
const int i1,
const int i2)
const {
895 return (*_ref_values)[i1] < (*_ref_values)[i2];
898 const std::vector<double> * _ref_values;
906 template <
class L>
inline PseudoJet::PseudoJet(
const L & some_four_vector) {
907 reset(some_four_vector);
911 inline void PseudoJet::_reset_indices() {
912 set_cluster_hist_index(-1);
920 inline double PseudoJet::m()
const {
922 return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
926 inline void PseudoJet::reset(
double px_in,
double py_in,
double pz_in,
double E_in) {
935 inline void PseudoJet::reset_momentum(
double px_in,
double py_in,
double pz_in,
double E_in) {
943 inline void PseudoJet::reset_momentum(
const PseudoJet & pj) {
962 template<
typename StructureType>
963 const StructureType & PseudoJet::structure()
const{
964 return dynamic_cast<const StructureType &
>(* validated_structure_ptr());
970 template<
typename TransformerType>
971 bool PseudoJet::has_structure_of()
const{
972 if (!_structure())
return false;
974 return dynamic_cast<const typename TransformerType::StructureType *
>(_structure.get()) != 0;
980 template<
typename TransformerType>
981 const typename TransformerType::StructureType & PseudoJet::structure_of()
const{
983 throw Error(
"Trying to access the structure of a PseudoJet without an associated structure");
985 return dynamic_cast<const typename TransformerType::StructureType &
>(*_structure);
1002 PseudoJet join(
const std::vector<PseudoJet> & pieces);
1018 FASTJET_END_NAMESPACE
1020 #endif // __FASTJET_PSEUDOJET_HH__