32#ifndef __FASTJET_PSEUDOJET_HH__
33#define __FASTJET_PSEUDOJET_HH__
40#include "fastjet/config.h"
41#include "fastjet/internal/numconsts.hh"
42#include "fastjet/internal/IsBase.hh"
43#include "fastjet/SharedPtr.hh"
44#include "fastjet/Error.hh"
45#include "fastjet/PseudoJetStructureBase.hh"
47FASTJET_BEGIN_NAMESPACE
57const double pseudojet_invalid_rap = -1e200;
82 PseudoJet() : _px(0), _py(0), _pz(0), _E(0) {_finish_init(); _reset_indices();}
85 PseudoJet(
const double px,
const double py,
const double pz,
const double E);
89 template <
class L>
PseudoJet(
const L & some_four_vector);
98#ifdef FASTJET_HAVE_THREAD_SAFETY
99 PseudoJet(
const PseudoJet &other){ (*this)=other; }
100 PseudoJet& operator=(
const PseudoJet& other);
116 inline double E()
const {
return _E;}
117 inline double e()
const {
return _E;}
118 inline double px()
const {
return _px;}
119 inline double py()
const {
return _py;}
120 inline double pz()
const {
return _pz;}
123 inline double phi()
const {
return phi_02pi();}
127 _ensure_valid_rap_phi();
128 return _phi > pi ? _phi-twopi : _phi;}
132 _ensure_valid_rap_phi();
138 inline double rap()
const {
139 _ensure_valid_rap_phi();
148 double pseudorapidity()
const;
149 double eta()
const {
return pseudorapidity();}
152 inline double pt2()
const {
return _kt2;}
154 inline double pt()
const {
return sqrt(_kt2);}
156 inline double perp2()
const {
return _kt2;}
158 inline double perp()
const {
return sqrt(_kt2);}
160 inline double kt2()
const {
return _kt2;}
163 inline double m2()
const {
return (_E+_pz)*(_E-_pz)-_kt2;}
166 inline double m()
const;
169 inline double mperp2()
const {
return (_E+_pz)*(_E-_pz);}
171 inline double mperp()
const {
return sqrt(std::abs(mperp2()));}
173 inline double mt2()
const {
return (_E+_pz)*(_E-_pz);}
175 inline double mt()
const {
return sqrt(std::abs(mperp2()));}
178 inline double modp2()
const {
return _kt2+_pz*_pz;}
180 inline double modp()
const {
return sqrt(_kt2+_pz*_pz);}
183 inline double Et()
const {
return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);}
185 inline double Et2()
const {
return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);}
190 return std::min(1.0, std::max(-1.0, _pz/sqrt(modp2())));
196 double operator () (
int i)
const ;
198 inline double operator [] (
int i)
const {
return (*
this)(i); };
203 double kt_distance(
const PseudoJet & other)
const;
206 double plain_distance(
const PseudoJet & other)
const;
210 return plain_distance(other);}
215 return sqrt(squared_distance(other));
220 double delta_phi_to(
const PseudoJet & other)
const;
232 std::valarray<double> four_mom()
const;
237 enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
246 PseudoJet & boost(
const PseudoJet & prest);
249 PseudoJet & unboost(
const PseudoJet & prest);
251 PseudoJet & operator*=(
double);
252 PseudoJet & operator/=(
double);
253 PseudoJet & operator+=(
const PseudoJet &);
254 PseudoJet & operator-=(
const PseudoJet &);
290 inline void reset(
double px,
double py,
double pz,
double E);
307 template <
class L>
inline void reset(
const L & some_four_vector) {
319 const PseudoJet * pj = fastjet::cast_if_derived<const PseudoJet>(&some_four_vector);
324 reset(some_four_vector[0], some_four_vector[1],
325 some_four_vector[2], some_four_vector[3]);
333 inline void reset_PtYPhiM(
double pt_in,
double y_in,
double phi_in,
double m_in=0.0) {
334 reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in);
342 inline void reset_momentum(
double px,
double py,
double pz,
double E);
347 inline void reset_momentum(
const PseudoJet & pj);
351 void reset_momentum_PtYPhiM(
double pt,
double y,
double phi,
double m=0.0);
357 reset_momentum(some_four_vector[0], some_four_vector[1],
358 some_four_vector[2], some_four_vector[3]);
373 void set_cached_rap_phi(
double rap,
double phi);
446 _user_info.reset(user_info_in);
480 return dynamic_cast<const L &
>(* _user_info.get());
485 return _user_info.get();
492 return _user_info.get() &&
dynamic_cast<const L *
>(_user_info.get());
500 return _user_info.get();
527 _user_info = user_info_in;
542 std::string description()
const;
556 bool has_associated_cluster_sequence()
const;
562 bool has_valid_cluster_sequence()
const;
570 const ClusterSequence* associated_cs()
const {
return associated_cluster_sequence();}
575 return validated_cs();
584 return validated_csab();
607 bool has_structure()
const;
641 template<
typename StructureType>
642 const StructureType & structure()
const;
647 template<
typename TransformerType>
648 bool has_structure_of()
const;
655 template<
typename TransformerType>
656 const typename TransformerType::StructureType & structure_of()
const;
675 virtual bool has_partner(
PseudoJet &partner)
const;
683 virtual bool has_child(
PseudoJet &child)
const;
698 virtual bool contains(
const PseudoJet &constituent)
const;
705 virtual bool is_inside(
const PseudoJet &jet)
const;
709 virtual bool has_constituents()
const;
715 virtual std::vector<PseudoJet> constituents()
const;
719 virtual bool has_exclusive_subjets()
const;
732 std::vector<PseudoJet> exclusive_subjets (
const double dcut)
const;
740 int n_exclusive_subjets(
const double dcut)
const;
750 std::vector<PseudoJet> exclusive_subjets (
int nsub)
const;
760 std::vector<PseudoJet> exclusive_subjets_up_to (
int nsub)
const;
769 double exclusive_subdmerge(
int nsub)
const;
779 double exclusive_subdmerge_max(
int nsub)
const;
789 virtual bool has_pieces()
const;
795 virtual std::vector<PseudoJet> pieces()
const;
804 virtual bool has_area()
const;
808 virtual double area()
const;
813 virtual double area_error()
const;
821 virtual bool is_pure_ghost()
const;
840 return cluster_hist_index();}
844 set_cluster_hist_index(index);}
856 double _px,_py,_pz,_E;
857 mutable double _phi, _rap;
859 int _cluster_hist_index, _user_index;
861#ifdef FASTJET_HAVE_THREAD_SAFETY
868 mutable std::atomic<int> _init_status;
874 void _reset_indices();
895#ifdef FASTJET_HAVE_THREAD_SAFETY
896 void _ensure_valid_rap_phi()
const;
898 inline void _ensure_valid_rap_phi()
const{
899 if (_phi == pseudojet_invalid_phi) _set_rap_phi();
904 void _set_rap_phi()
const;
907 friend PseudoJet operator*(
double,
const PseudoJet &);
914PseudoJet operator+(
const PseudoJet &,
const PseudoJet &);
915PseudoJet operator-(
const PseudoJet &,
const PseudoJet &);
916PseudoJet operator*(
double,
const PseudoJet &);
917PseudoJet operator*(
const PseudoJet &,
double);
918PseudoJet operator/(
const PseudoJet &,
double);
923bool operator==(
const PseudoJet &,
const PseudoJet &);
930bool operator==(
const PseudoJet & jet,
const double val);
931inline bool operator==(
const double val,
const PseudoJet & jet) {
return jet == val;}
936inline bool operator!=(
const double val,
const PseudoJet & a) {
return !(a==val);}
940 return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
945 double dot_3d = a.px()*b.px() + a.py()*b.py() + a.pz()*b.pz();
946 return std::min(1.0, std::max(-1.0, dot_3d/sqrt(a.
modp2()*b.
modp2())));
955bool have_same_momentum(
const PseudoJet &,
const PseudoJet &);
959PseudoJet PtYPhiM(
double pt,
double y,
double phi,
double m = 0.0);
965std::vector<PseudoJet>
sorted_by_pt(
const std::vector<PseudoJet> & jets);
971std::vector<PseudoJet>
sorted_by_E(
const std::vector<PseudoJet> & jets);
974std::vector<PseudoJet>
sorted_by_pz(
const std::vector<PseudoJet> & jets);
981void sort_indices(std::vector<int> & indices,
982 const std::vector<double> & values);
989 const std::vector<double> & values) {
991 if (objects.size() != values.size()){
992 throw Error(
"fastjet::objects_sorted_by_values(...): the size of the 'objects' vector must match the size of the 'values' vector");
996 std::vector<int> indices(values.size());
997 for (
size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
1000 sort_indices(indices, values);
1003 std::vector<T> objects_sorted(objects.size());
1006 for (
size_t i = 0; i < indices.size(); i++) {
1007 objects_sorted[i] = objects[indices[i]];
1010 return objects_sorted;
1021 _ref_values = reference_values;
1023 inline int operator() (
const int i1,
const int i2)
const {
1024 return (*_ref_values)[i1] < (*_ref_values)[i2];
1027 const std::vector<double> * _ref_values;
1036template <
class L>
inline PseudoJet::PseudoJet(
const L & some_four_vector) {
1037 reset(some_four_vector);
1042inline void PseudoJet::_reset_indices() {
1043 set_cluster_hist_index(-1);
1064inline double PseudoJet::m()
const {
1066 return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
1070inline void PseudoJet::reset(
double px_in,
double py_in,
double pz_in,
double E_in) {
1080inline void PseudoJet::reset_momentum(
double px_in,
double py_in,
double pz_in,
double E_in) {
1088inline void PseudoJet::reset_momentum(
const PseudoJet & pj) {
1093#ifdef FASTJET_HAVE_THREAD_SAFETY
1102 int expected = Init_Done;
1103 if (pj._init_status.compare_exchange_weak(expected, Init_Done)) {
1104 _init_status = Init_Done;
1127template<
typename StructureType>
1128const StructureType & PseudoJet::structure()
const{
1129 return dynamic_cast<const StructureType &
>(* validated_structure_ptr());
1135template<
typename TransformerType>
1136bool PseudoJet::has_structure_of()
const{
1137 if (!_structure)
return false;
1139 return dynamic_cast<const typename TransformerType::StructureType *
>(_structure.
get()) != 0;
1145template<
typename TransformerType>
1146const typename TransformerType::StructureType & PseudoJet::structure_of()
const{
1148 throw Error(
"Trying to access the structure of a PseudoJet without an associated structure");
1150 return dynamic_cast<const typename TransformerType::StructureType &
>(*_structure);
1183FASTJET_END_NAMESPACE
base class that sets interface for extensions of ClusterSequence that provide information about the a...
base class corresponding to errors that can be thrown by FastJet
Contains any information related to the clustering that should be directly accessible to PseudoJet.
error class to be thrown if accessing user info when it doesn't exist
a base class to hold extra user information in a PseudoJet
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
PseudoJet()
default constructor, which as of FJ3.0 provides an object for which all operations are now valid and ...
double cos_theta() const
cos of the polar angle should we have: min(1.0,max(-1.0,_pz/sqrt(modp2())));
const ClusterSequenceAreaBase * validated_cluster_sequence_area_base() const
if the jet has valid area information then return a pointer to the associated ClusterSequenceAreaBase...
double mt() const
returns the transverse mass = sqrt(kt^2+m^2)
int cluster_sequence_history_index() const
alternative name for cluster_hist_index() [perhaps more meaningful]
double rap() const
returns the rapidity or some large value when the rapidity is infinite
double Et() const
return the transverse energy
double delta_R(const PseudoJet &other) const
return the cylinder (rap-phi) distance between this jet and another, .
const L & user_info() const
returns a reference to the dynamic cast conversion of user_info to type L.
double mperp() const
returns the transverse mass = sqrt(kt^2+m^2)
double modp() const
return the 3-vector modulus = sqrt(px^2+py^2+pz^2)
void set_user_index(const int index)
set the user_index, intended to allow the user to add simple identifying information to a particle/je...
const UserInfoBase * user_info_ptr() const
retrieve a pointer to the (const) user information
void reset_momentum(const L &some_four_vector)
reset the 4-momentum according to the supplied generic 4-vector (accessible via indexing,...
bool has_user_info() const
returns true if the PseudoJet has user information
double squared_distance(const PseudoJet &other) const
returns squared cylinder (rap-phi) distance between this jet and another
bool has_valid_cs() const
shorthand for has_valid_cluster_sequence()
double phi() const
returns phi (in the range 0..2pi)
virtual ~PseudoJet()
default (virtual) destructor
double perp() const
returns the scalar transverse momentum
double perp2() const
returns the squared transverse momentum
double m2() const
returns the squared invariant mass // like CLHEP
void set_cluster_hist_index(const int index)
set the cluster_hist_index, intended to be used by clustering routines.
double kt2() const
returns the squared transverse momentum
bool has_associated_cs() const
shorthand for has_associated_cluster_sequence()
bool has_user_info() const
returns true if the PseudoJet has user information than can be cast to the template argument type.
double theta() const
polar angle
double phi_std() const
returns phi in the range -pi..pi
double phi_02pi() const
returns phi in the range 0..2pi
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
void reset_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in=0.0)
reset the PseudoJet according to the specified pt, rapidity, azimuth and mass (also resetting indices...
SharedPtr< UserInfoBase > & user_info_shared_ptr()
retrieve a (non-const) shared pointer to the user information; you can use this, for example,...
double pt() const
returns the scalar transverse momentum
void set_cluster_sequence_history_index(const int index)
alternative name for set_cluster_hist_index(...) [perhaps more meaningful]
double mperp2() const
returns the squared transverse mass = kt^2+m^2
void reset(const L &some_four_vector)
reset the 4-momentum according to the supplied generic 4-vector (accessible via indexing,...
int user_index() const
return the user_index,
const SharedPtr< UserInfoBase > & user_info_shared_ptr() const
retrieve a (const) shared pointer to the user information
const ClusterSequence * validated_cluster_sequence() const
if the jet has a valid associated cluster sequence then return a pointer to it; otherwise throw an er...
double modp2() const
return the squared 3-vector modulus = px^2+py^2+pz^2
double Et2() const
return the transverse energy squared
double rapidity() const
the same as rap()
double beam_distance() const
returns distance between this jet and the beam
double pt2() const
returns the squared transverse momentum
void set_user_info(UserInfoBase *user_info_in)
sets the internal shared pointer to the user information.
void reset(const PseudoJet &psjet)
reset the PseudoJet to be equal to psjet (including its indices); NB if the argument is derived from ...
double mt2() const
returns the squared transverse mass = kt^2+m^2
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
T * get() const
get the stored pointer
void reset()
reset the pointer to default value (NULL)
bool operator!=(const PseudoJet &a, const PseudoJet &b)
inequality test which is exact opposite of operator==
vector< PseudoJet > sorted_by_rapidity(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing rapidity
const double MaxRap
Used to protect against parton-level events where pt can be zero for some partons,...
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
bool operator==(const PseudoJet &a, const PseudoJet &b)
returns true if the 4 momentum components of the two PseudoJets are identical and all the internal in...
std::vector< T > objects_sorted_by_values(const std::vector< T > &objects, const std::vector< double > &values)
given a vector of values with a one-to-one correspondence with the vector of objects,...
double dot_product(const PseudoJet &a, const PseudoJet &b)
returns the 4-vector dot product of a and b
const double pseudojet_invalid_phi
default value for phi, meaning it (and rapidity) have yet to be calculated)
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
double theta(const PseudoJet &a, const PseudoJet &b)
returns the angle between a and b
double cos_theta(const PseudoJet &a, const PseudoJet &b)
returns the cosine of the angle between a and b
vector< PseudoJet > sorted_by_pz(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing pz