EIC code displayed by LXR master/​include/​fastjet/​PseudoJet.hh	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-09-17 08:57:14

 //FJSTARTHEADER
 // $Id$
 //
 // Copyright (c) 2005-2025, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 //
 //----------------------------------------------------------------------
 // This file is part of FastJet.
 //
 //  FastJet is free software; you can redistribute it and/or modify
 //  it under the terms of the GNU General Public License as published by
 //  the Free Software Foundation; either version 2 of the License, or
 //  (at your option) any later version.
 //
 //  The algorithms that underlie FastJet have required considerable
 //  development. They are described in the original FastJet paper,
 //  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
 //  FastJet as part of work towards a scientific publication, please
 //  quote the version you use and include a citation to the manual and
 //  optionally also to hep-ph/0512210.
 //
 //  FastJet is distributed in the hope that it will be useful,
 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 //  GNU General Public License for more details.
 //
 //  You should have received a copy of the GNU General Public License
 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
 //----------------------------------------------------------------------
 //FJENDHEADER
 
 
 #ifndef __FASTJET_PSEUDOJET_HH__
 #define __FASTJET_PSEUDOJET_HH__
 
 #include<valarray>
 #include<vector>
 #include<cassert>
 #include<cmath>
 #include<iostream>
 #include "fastjet/config.h"
 #include "fastjet/internal/numconsts.hh"
 #include "fastjet/internal/IsBase.hh"
 #include "fastjet/SharedPtr.hh"
 #include "fastjet/Error.hh"
 #include "fastjet/PseudoJetStructureBase.hh"
 
 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
 //using namespace std;
 
 /// Used to protect against parton-level events where pt can be zero
 /// for some partons, giving rapidity=infinity. KtJet fails in those cases.
 const double MaxRap = 1e5;
 
 /// default value for phi, meaning it (and rapidity) have yet to be calculated) 
 const double pseudojet_invalid_phi = -100.0;
 const double pseudojet_invalid_rap = -1e200;
 
 #ifndef __FJCORE__
 // forward definition
 class ClusterSequenceAreaBase;
 #endif  // __FJCORE__
 
 /// @ingroup basic_classes
 /// \class PseudoJet
 /// Class to contain pseudojets, including minimal information of use to
 /// jet-clustering routines.
 class PseudoJet {
 
  public:
   //----------------------------------------------------------------------
   /// @name Constructors and destructor
   //\{
   /// default constructor, which as of FJ3.0 provides an object for
   /// which all operations are now valid and which has zero momentum
   ///
   // (cf. this is actually OK from a timing point of view and in some
   // cases better than just having the default constructor for the
   // internal shared pointer: see PJtiming.cc and the notes therein)
   //
   // note: no reset of shared pointers needed
   PseudoJet() : _px(0), _py(0), _pz(0), _E(0) {_finish_init(); _reset_indices();}
   //PseudoJet() : _px(0), _py(0), _pz(0), _E(0), _phi(pseudojet_invalid_phi), _rap(pseudojet_invalid_rap), _kt2(0) {_reset_indices();}
   /// construct a pseudojet from explicit components
   PseudoJet(const double px, const double py, const double pz, const double E);
 
   /// constructor from any object that has px,py,pz,E = some_four_vector[0--3],
   #ifndef SWIG
   template <class L> PseudoJet(const L & some_four_vector);
   #endif
   
   // Constructor that performs minimal initialisation (only that of
   // the shared pointers), of use in certain speed-critical contexts
   //
   // NB: "dummy" is commented to avoid unused-variable compiler warnings
   PseudoJet(bool /* dummy */) {}
 
 #ifdef FASTJET_HAVE_THREAD_SAFETY
   PseudoJet(const PseudoJet &other){ (*this)=other; }
   PseudoJet& operator=(const PseudoJet& other);
 #endif
   
   /// default (virtual) destructor
   virtual ~PseudoJet(){}
 //std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
 //std::shared_ptr:     ;
 //std::shared_ptr: #else
 //std::shared_ptr:     {}
 //std::shared_ptr: #endif
   //\} ---- end of constructors and destructors --------------------------
 
   //----------------------------------------------------------------------
   /// @name Kinematic access functions
   //\{
   //----------------------------------------------------------------------
   inline double E()   const {return _E;}
   inline double e()   const {return _E;} // like CLHEP
   inline double px()  const {return _px;}
   inline double py()  const {return _py;}
   inline double pz()  const {return _pz;}
 
   /// returns phi (in the range 0..2pi)
   inline double phi() const {return phi_02pi();}
 
   /// returns phi in the range -pi..pi
   inline double phi_std()  const {
     _ensure_valid_rap_phi();
     return _phi > pi ? _phi-twopi : _phi;}
 
   /// returns phi in the range 0..2pi
   inline double phi_02pi() const {
     _ensure_valid_rap_phi();
     return _phi;
   }
 
   /// returns the rapidity or some large value when the rapidity
   /// is infinite
   inline double rap() const {
     _ensure_valid_rap_phi();
     return _rap;
   }
 
   /// the same as rap()
   inline double rapidity() const {return rap();} // like CLHEP
 
   /// returns the pseudo-rapidity or some large value when the
   /// rapidity is infinite
   double pseudorapidity() const;
   double eta() const {return pseudorapidity();}
 
   /// returns the squared transverse momentum
   inline double pt2() const {return _kt2;}
   /// returns the scalar transverse momentum
   inline double  pt() const {return sqrt(_kt2);} 
   /// returns the squared transverse momentum
   inline double perp2() const {return _kt2;}  // like CLHEP
   /// returns the scalar transverse momentum
   inline double  perp() const {return sqrt(_kt2);}    // like CLHEP
   /// returns the squared transverse momentum
   inline double kt2() const {return _kt2;} // for bkwds compatibility
 
   /// returns the squared invariant mass // like CLHEP
   inline double  m2() const {return (_E+_pz)*(_E-_pz)-_kt2;}    
   /// returns the invariant mass 
   /// (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP)
   inline double  m() const;    
 
   /// returns the squared transverse mass = kt^2+m^2
   inline double mperp2() const {return (_E+_pz)*(_E-_pz);}
   /// returns the transverse mass = sqrt(kt^2+m^2)
   inline double mperp() const {return sqrt(std::abs(mperp2()));}
   /// returns the squared transverse mass = kt^2+m^2
   inline double mt2() const {return (_E+_pz)*(_E-_pz);}
   /// returns the transverse mass = sqrt(kt^2+m^2)
   inline double mt() const {return sqrt(std::abs(mperp2()));}
 
   /// return the squared 3-vector modulus = px^2+py^2+pz^2
   inline double modp2() const {return _kt2+_pz*_pz;}
   /// return the 3-vector modulus = sqrt(px^2+py^2+pz^2)
   inline double modp() const {return sqrt(_kt2+_pz*_pz);}
 
   /// return the transverse energy
   inline double Et() const {return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);}
   /// return the transverse energy squared
   inline double Et2() const {return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);}
 
   /// cos of the polar angle
   /// should we have: min(1.0,max(-1.0,_pz/sqrt(modp2()))); 
   inline double cos_theta() const {
     return std::min(1.0, std::max(-1.0, _pz/sqrt(modp2())));
   }
   /// polar angle
   inline double theta() const { return acos(cos_theta()); }
 
   /// returns component i, where X==0, Y==1, Z==2, E==3
   double operator () (int i) const ; 
   /// returns component i, where X==0, Y==1, Z==2, E==3
   inline double operator [] (int i) const { return (*this)(i); }; // this too
 
 
 
   /// returns kt distance (R=1) between this jet and another
   double kt_distance(const PseudoJet & other) const;
 
   /// returns squared cylinder (rap-phi) distance between this jet and another
   double plain_distance(const PseudoJet & other) const;
   /// returns squared cylinder (rap-phi) distance between this jet and
   /// another
   inline double squared_distance(const PseudoJet & other) const {
     return plain_distance(other);}
 
   /// return the cylinder (rap-phi) distance between this jet and another,
   /// \f$\Delta_R = \sqrt{\Delta y^2 + \Delta \phi^2}\f$.
   inline double delta_R(const PseudoJet & other) const {
     return sqrt(squared_distance(other));
   }
 
   /// returns other.phi() - this.phi(), constrained to be in 
   /// range -pi .. pi
   double delta_phi_to(const PseudoJet & other) const;
 
   //// this seemed to compile except if it was used
   //friend inline double 
   //  kt_distance(const PseudoJet & jet1, const PseudoJet & jet2) { 
   //                                      return jet1.kt_distance(jet2);}
 
   /// returns distance between this jet and the beam
   inline double beam_distance() const {return _kt2;}
 
   /// return a valarray containing the four-momentum (components 0-2
   /// are 3-mom, component 3 is energy).
   std::valarray<double> four_mom() const;
 
   //\}  ------- end of kinematic access functions
 
   // taken from CLHEP
   enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
 
 
   //----------------------------------------------------------------------
   /// @name Kinematic modification functions
   //\{
   //----------------------------------------------------------------------
   /// transform this jet (given in the rest frame of prest) into a jet
   /// in the lab frame
   PseudoJet & boost(const PseudoJet & prest);
   /// transform this jet (given in lab) into a jet in the rest
   /// frame of prest
   PseudoJet & unboost(const PseudoJet & prest);
 
   PseudoJet & operator*=(double);
   PseudoJet & operator/=(double);
   PseudoJet & operator+=(const PseudoJet &);
   PseudoJet & operator-=(const PseudoJet &);
 
 //std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
 //std::shared_ptr:   /// overload the assignment through the = operator
 //std::shared_ptr:   ///
 //std::shared_ptr:   /// this is needed to make sure that release_from_cs is called
 //std::shared_ptr:   /// before copying the structure shared pointer!
 //std::shared_ptr:   void operator=(const PseudoJet &other){
 //std::shared_ptr:     _release_jet_from_cs();
 //std::shared_ptr:     _structure = other._structure;
 //std::shared_ptr:     _user_info = other._user_info;
 //std::shared_ptr: 
 //std::shared_ptr:     _px = other._px;
 //std::shared_ptr:     _py = other._py;
 //std::shared_ptr:     _pz = other._pz;
 //std::shared_ptr:     _E  = other._E;
 //std::shared_ptr: 
 //std::shared_ptr:     _phi = other._phi;
 //std::shared_ptr:     _rap = other._rap;
 //std::shared_ptr:     _kt2  = other._kt2;
 //std::shared_ptr:     
 //std::shared_ptr:     _cluster_hist_index = other._cluster_hist_index;
 //std::shared_ptr:     _user_index = other._user_index;
 //std::shared_ptr:   }
 //std::shared_ptr: 
 //std::shared_ptr:   /// force resetting the structure pointer to an empty structure
 //std::shared_ptr:   ///
 //std::shared_ptr:   /// THIS IS STRICTLY MEANT FOR INTERNAL USAGE IN SELF-DELETING
 //std::shared_ptr:   /// ClusterSequences. IT SHOULD NOT BE USED BY END-USERS.
 //std::shared_ptr:   void force_reset_structure(){
 //std::shared_ptr:     _structure.reset();
 //std::shared_ptr:   }  
 //std::shared_ptr: #endif
 
   /// reset the 4-momentum according to the supplied components and
   /// put the user and history indices back to their default values
   inline void reset(double px, double py, double pz, double E);
 
   /// reset the PseudoJet to be equal to psjet (including its
   /// indices); NB if the argument is derived from a PseudoJet then
   /// the "reset" used will be the templated version
   ///
   /// Note: this is included on top of the templated version because
   /// PseudoJet is not "derived" from PseudoJet, so the templated
   /// reset would not handle this case properly.
   inline void reset(const PseudoJet & psjet) {
     (*this) = psjet;
   }
 
   /// reset the 4-momentum according to the supplied generic 4-vector
   /// (accessible via indexing, [0]==px,...[3]==E) and put the user
   /// and history indices back to their default values.
 #ifndef SWIG
   template <class L> inline void reset(const L & some_four_vector) {
     // check if some_four_vector can be cast to a PseudoJet
     //
     // Note that a regular dynamic_cast would not work here because
     // there is no guarantee that L is polymorphic. We use a more
     // complex construct here that works also in such a case. As for
     // dynamic_cast, NULL is returned if L is not derived from
     // PseudoJet
     //
     // Note the explicit request for fastjet::cast_if_derived; when
     // combining fastjet and fjcore, this avoids ambiguity in which of
     // the two cast_if_derived calls to use.
     const PseudoJet * pj = fastjet::cast_if_derived<const PseudoJet>(&some_four_vector);
 
     if (pj){
       (*this) = *pj;
     } else {
       reset(some_four_vector[0], some_four_vector[1],
         some_four_vector[2], some_four_vector[3]);
     }
   }
 #endif // SWIG
 
   /// reset the PseudoJet according to the specified pt, rapidity,
   /// azimuth and mass (also resetting indices, etc.)
   /// (phi should satisfy -2pi<phi<4pi)
   inline void reset_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in=0.0) {
     reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in);
     _reset_indices();
     //std::shared_ptr: _reset_shared_pointers();
   }
 
   /// reset the 4-momentum according to the supplied components 
   /// but leave all other information (indices, user info, etc.)
   /// untouched
   inline void reset_momentum(double px, double py, double pz, double E);
 
   /// reset the 4-momentum according to the components of the supplied
   /// PseudoJet, including cached components; note that the template
   /// version (below) will be called for classes derived from PJ.
   inline void reset_momentum(const PseudoJet & pj);
 
   /// reset the 4-momentum according to the specified pt, rapidity,
   /// azimuth and mass (phi should satisfy -2pi<phi<4pi)
   void reset_momentum_PtYPhiM(double pt, double y, double phi, double m=0.0);
 
   /// reset the 4-momentum according to the supplied generic 4-vector
   /// (accessible via indexing, [0]==px,...[3]==E), but leave all
   /// other information (indices, user info, etc.)  untouched
   template <class L> inline void reset_momentum(const L & some_four_vector) {
     reset_momentum(some_four_vector[0], some_four_vector[1],
            some_four_vector[2], some_four_vector[3]);
   }
 
   /// in some cases when setting a 4-momentum, the user/program knows
   /// what rapidity and azimuth are associated with that 4-momentum;
   /// by calling this routine the user can provide the information
   /// directly to the PseudoJet and avoid expensive rap-phi
   /// recalculations.
   ///
   /// - \param rap  rapidity
   /// - \param phi  (in range -twopi...4*pi)
   ///
   /// USE WITH CAUTION: there are no checks that the rapidity and
   /// azimuth supplied are sensible, nor does this reset the
   /// 4-momentum components if things don't match.
   void set_cached_rap_phi(double rap, double phi);
 
 
   //\} --- end of kin mod functions ------------------------------------
 
   //----------------------------------------------------------------------
   /// @name User index functions
   ///
   /// To allow the user to set and access an integer index which can
   /// be exploited by the user to associate extra information with a
   /// particle/jet (for example pdg id, or an indication of a
   /// particle's origin within the user's analysis)
   //
   //\{
 
   /// return the user_index, 
   inline int user_index() const {return _user_index;}
   /// set the user_index, intended to allow the user to add simple
   /// identifying information to a particle/jet
   inline void set_user_index(const int index) {_user_index = index;}
 
   //\} ----- end of use index functions ---------------------------------
 
   //----------------------------------------------------------------------
   /// @name User information types and functions
   ///
   /// Allows PseudoJet to carry extra user info (as an object derived from
   /// UserInfoBase).
   //\{
 
   /// @ingroup user_info
   /// \class UserInfoBase
   /// a base class to hold extra user information in a PseudoJet
   ///
   /// This is a base class to help associate extra user information
   /// with a jet. The user should store their information in a class
   /// derived from this. This allows information of arbitrary
   /// complexity to be easily associated with a PseudoJet (in contrast
   /// to the user index). For example, in a Monte Carlo simulation,
   /// the user information might include the PDG ID, and the position
   /// of the production vertex for the particle.
   ///
   /// The PseudoJet is able to store a shared pointer to any object
   /// derived from UserInfo. The use of a shared pointer frees the
   /// user of the need to handle the memory management associated with
   /// the information.
   ///
   /// Having the user information derive from a common base class also
   /// facilitates dynamic casting, etc.
   ///
   class UserInfoBase{
   public:
     // dummy ctor
     UserInfoBase(){};
 
     // dummy virtual dtor
     // makes it polymorphic to allow for dynamic_cast
     virtual ~UserInfoBase(){}; 
   };
 
   /// error class to be thrown if accessing user info when it doesn't
   /// exist
   class InexistentUserInfo : public Error {
   public:
     InexistentUserInfo();
   };
 
   /// sets the internal shared pointer to the user information.
   ///
   /// Note that the PseudoJet will now _own_ the pointer, and delete
   /// the corresponding object when it (the jet, and any copies of the jet)
   /// goes out of scope. 
   void set_user_info(UserInfoBase * user_info_in) {
     _user_info.reset(user_info_in);
   }
 
   /// returns a reference to the dynamic cast conversion of user_info
   /// to type L.
   ///
   /// Usage: suppose you have previously set the user info with a pointer
   /// to an object of type MyInfo, 
   ///
   ///   class MyInfo: public PseudoJet::UserInfoBase {
   ///      MyInfo(int id) : _pdg_id(id);
   ///      int pdg_id() const {return _pdg_id;}
   ///      int _pdg_id;
   ///   };
   ///
   ///   PseudoJet particle(...);
   ///   particle.set_user_info(new MyInfo(its_pdg_id));
   ///
   /// Then you would access that pdg_id() as
   ///
   ///   particle.user_info<MyInfo>().pdg_id();
   ///
   /// It's overkill for just a single integer, but scales easily to
   /// more extensive information.
   ///
   /// Note that user_info() throws an InexistentUserInfo() error if
   /// there is no user info; throws a std::bad_cast if the conversion
   /// doesn't work
   ///
   /// If this behaviour does not fit your needs, use instead the the
   /// user_info_ptr() or user_info_shared_ptr() member functions.
   template<class L>
   const L & user_info() const{
     if (_user_info.get() == 0) throw InexistentUserInfo();
     return dynamic_cast<const L &>(* _user_info.get());
   }
 
   /// returns true if the PseudoJet has user information
   bool has_user_info() const{
     return _user_info.get();
   }
 
   /// returns true if the PseudoJet has user information than can be
   /// cast to the template argument type.
   template<class L>
   bool has_user_info() const{
     return _user_info.get() && dynamic_cast<const L *>(_user_info.get());
   }
 
   /// retrieve a pointer to the (const) user information
   const UserInfoBase * user_info_ptr() const{
     // the line below is not needed since the next line would anyway
     // return NULL in that case
     //if (!_user_info) return NULL;
     return _user_info.get();
   }
 
 
   /// retrieve a (const) shared pointer to the user information
   const SharedPtr<UserInfoBase> & user_info_shared_ptr() const{
     return _user_info;
   }
 
   /// retrieve a (non-const) shared pointer to the user information;
   /// you can use this, for example, to set the shared pointer, eg
   ///
   /// \code
   ///   p2.user_info_shared_ptr() = p1.user_info_shared_ptr();
   /// \endcode
   ///
   /// or 
   ///
   /// \code
   ///   SharedPtr<PseudoJet::UserInfoBase> info_shared(new MyInfo(...));
   ///   p2.user_info_shared_ptr() = info_shared;
   /// \endcode
   SharedPtr<UserInfoBase> & user_info_shared_ptr(){
     return _user_info;
   }
 
   void set_user_info_shared_ptr(const SharedPtr<UserInfoBase> & user_info_in) {
     _user_info = user_info_in;
   }
 
   // \} --- end of extra info functions ---------------------------------
 
   //----------------------------------------------------------------------
   /// @name Description
   ///
   /// Since a PseudoJet can have a structure that contains a variety
   /// of information, we provide a description that allows one to check
   /// exactly what kind of PseudoJet we are dealing with
   //
   //\{
 
   /// return a string describing what kind of PseudoJet we are dealing with 
   std::string description() const;
 
   //\} ----- end of description functions ---------------------------------
 
   //-------------------------------------------------------------
   /// @name Access to the associated ClusterSequence object.
   ///
   /// In addition to having kinematic information, jets may contain a
   /// reference to an associated ClusterSequence (this is the case,
   /// for example, if the jet has been returned by a ClusterSequence
   /// member function).
   //\{
   //-------------------------------------------------------------
   /// returns true if this PseudoJet has an associated ClusterSequence.
   bool has_associated_cluster_sequence() const;
   /// shorthand for has_associated_cluster_sequence()
   bool has_associated_cs() const {return has_associated_cluster_sequence();}
 
   /// returns true if this PseudoJet has an associated and still
   /// valid(ated) ClusterSequence.
   bool has_valid_cluster_sequence() const;
   /// shorthand for has_valid_cluster_sequence()
   bool has_valid_cs() const {return has_valid_cluster_sequence();}
 
   /// get a (const) pointer to the parent ClusterSequence (NULL if
   /// inexistent)
   const ClusterSequence* associated_cluster_sequence() const;
   // shorthand for associated_cluster_sequence()
   const ClusterSequence* associated_cs() const {return associated_cluster_sequence();}
 
   /// if the jet has a valid associated cluster sequence then return a
   /// pointer to it; otherwise throw an error
   inline const ClusterSequence * validated_cluster_sequence() const {
     return validated_cs();
   }
   /// shorthand for validated_cluster_sequence()
   const ClusterSequence * validated_cs() const;
 
 #ifndef __FJCORE__
   /// if the jet has valid area information then return a pointer to
   /// the associated ClusterSequenceAreaBase object; otherwise throw an error
   inline const ClusterSequenceAreaBase * validated_cluster_sequence_area_base() const {
     return validated_csab();
   }
 
   /// shorthand for validated_cluster_sequence_area_base()
   const ClusterSequenceAreaBase * validated_csab() const;
 #endif  //  __FJCORE__
 
   //\}
 
   //-------------------------------------------------------------
   /// @name Access to the associated PseudoJetStructureBase object.
   ///
   /// In addition to having kinematic information, jets may contain a
   /// reference to an associated ClusterSequence (this is the case,
   /// for example, if the jet has been returned by a ClusterSequence
   /// member function).
   //\{
   //-------------------------------------------------------------
 
   /// set the associated structure
   void set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure_in);
 
   /// return true if there is some structure associated with this PseudoJet
   bool has_structure() const;
 
   /// return a pointer to the structure (of type
   /// PseudoJetStructureBase*) associated with this PseudoJet.
   ///
   /// return NULL if there is no associated structure
   const PseudoJetStructureBase* structure_ptr() const;
   
   /// return a non-const pointer to the structure (of type
   /// PseudoJetStructureBase*) associated with this PseudoJet.
   ///
   /// return NULL if there is no associated structure
   ///
   /// Only use this if you know what you are doing. In any case,
   /// prefer the 'structure_ptr()' (the const version) to this method,
   /// unless you really need a write access to the PseudoJet's
   /// underlying structure.
   PseudoJetStructureBase* structure_non_const_ptr();
   
   /// return a pointer to the structure (of type
   /// PseudoJetStructureBase*) associated with this PseudoJet.
   ///
   /// throw an error if there is no associated structure
   const PseudoJetStructureBase* validated_structure_ptr() const;
   
   /// return a reference to the shared pointer to the
   /// PseudoJetStructureBase associated with this PseudoJet
   const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const;
 
   /// returns a reference to the structure casted to the requested
   /// structure type
   ///
   /// If there is no structure associated, an Error is thrown.
   /// If the type is not met, a std::bad_cast error is thrown.
   template<typename StructureType>
   const StructureType & structure() const;
 
   /// check if the PseudoJet has the structure resulting from a Transformer 
   /// (that is, its structure is compatible with a Transformer::StructureType).
   /// If there is no structure, false is returned.
   template<typename TransformerType>
   bool has_structure_of() const;
 
   /// this is a helper to access any structure created by a Transformer 
   /// (that is, of type Transformer::StructureType).
   ///
   /// If there is no structure, or if the structure is not compatible
   /// with TransformerType, an error is thrown.
   template<typename TransformerType>
   const typename TransformerType::StructureType & structure_of() const;
 
   //\}
 
   //-------------------------------------------------------------
   /// @name Methods for access to information about jet structure
   ///
   /// These allow access to jet constituents, and other jet
   /// subtructure information. They only work if the jet is associated
   /// with a ClusterSequence.
   //-------------------------------------------------------------
   //\{
 
   /// check if it has been recombined with another PseudoJet in which
   /// case, return its partner through the argument. Otherwise,
   /// 'partner' is set to 0.
   ///
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   virtual bool has_partner(PseudoJet &partner) const;
 
   /// check if it has been recombined with another PseudoJet in which
   /// case, return its child through the argument. Otherwise, 'child'
   /// is set to 0.
   /// 
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   virtual bool has_child(PseudoJet &child) const;
 
   /// check if it is the product of a recombination, in which case
   /// return the 2 parents through the 'parent1' and 'parent2'
   /// arguments. Otherwise, set these to 0.
   ///
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const;
 
   /// check if the current PseudoJet contains the one passed as
   /// argument.
   ///
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   virtual bool contains(const PseudoJet &constituent) const;
 
   /// check if the current PseudoJet is contained in the one passed as
   /// argument.
   ///
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   virtual bool is_inside(const PseudoJet &jet) const;
 
 
   /// returns true if the PseudoJet has constituents
   virtual bool has_constituents() const;
 
   /// retrieve the constituents. 
   ///
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence or other substructure information
   virtual std::vector<PseudoJet> constituents() const;
 
 
   /// returns true if the PseudoJet has support for exclusive subjets
   virtual bool has_exclusive_subjets() const;
 
   /// return a vector of all subjets of the current jet (in the sense
   /// of the exclusive algorithm) that would be obtained when running
   /// the algorithm with the given dcut. 
   ///
   /// Time taken is O(m ln m), where m is the number of subjets that
   /// are found. If m gets to be of order of the total number of
   /// constituents in the jet, this could be substantially slower than
   /// just getting that list of constituents.
   ///
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   std::vector<PseudoJet> exclusive_subjets (const double dcut) const;
 
   /// return the size of exclusive_subjets(...); still n ln n with same
   /// coefficient, but marginally more efficient than manually taking
   /// exclusive_subjets.size()
   ///
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   int n_exclusive_subjets(const double dcut) const;
 
   /// return the list of subjets obtained by unclustering the supplied
   /// jet down to nsub subjets. Throws an error if there are fewer than
   /// nsub particles in the jet.
   ///
   /// For ClusterSequence type jets, requires nsub ln nsub time
   ///
   /// An Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   std::vector<PseudoJet> exclusive_subjets (int nsub) const;
 
   /// return the list of subjets obtained by unclustering the supplied
   /// jet down to nsub subjets (or all constituents if there are fewer
   /// than nsub).
   ///
   /// For ClusterSequence type jets, requires nsub ln nsub time
   ///
   /// An Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   std::vector<PseudoJet> exclusive_subjets_up_to (int nsub) const;
 
   /// Returns the dij that was present in the merging nsub+1 -> nsub 
   /// subjets inside this jet.
   ///
   /// Returns 0 if there were nsub or fewer constituents in the jet.
   ///
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   double exclusive_subdmerge(int nsub) const;
 
   /// Returns the maximum dij that occurred in the whole event at the
   /// stage that the nsub+1 -> nsub merge of subjets occurred inside 
   /// this jet.
   ///
   /// Returns 0 if there were nsub or fewer constituents in the jet.
   ///
   /// an Error is thrown if this PseudoJet has no currently valid
   /// associated ClusterSequence
   double exclusive_subdmerge_max(int nsub) const;
 
 
   /// returns true if a jet has pieces
   ///
   /// By default a single particle or a jet coming from a
   /// ClusterSequence have no pieces and this methos will return false.
   ///
   /// In practice, this is equivalent to have an structure of type
   /// CompositeJetStructure.
   virtual bool has_pieces() const;
 
 
   /// retrieve the pieces that make up the jet. 
   ///
   /// If the jet does not support pieces, an error is throw
   virtual std::vector<PseudoJet> pieces() const;
 
 
   // the following ones require a computation of the area in the
   // parent ClusterSequence (See ClusterSequenceAreaBase for details)
   //------------------------------------------------------------------
 #ifndef __FJCORE__
 
   /// check if it has a defined area
   virtual bool has_area() const;
 
   /// return the jet (scalar) area.
   /// throws an Error if there is no support for area in the parent CS
   virtual double area() const;
 
   /// return the error (uncertainty) associated with the determination
   /// of the area of this jet.
   /// throws an Error if there is no support for area in the parent CS
   virtual double area_error() const;
 
   /// return the jet 4-vector area.
   /// throws an Error if there is no support for area in the parent CS
   virtual PseudoJet area_4vector() const;
 
   /// true if this jet is made exclusively of ghosts.
   /// throws an Error if there is no support for area in the parent CS
   virtual bool is_pure_ghost() const;
 
 #endif  // __FJCORE__
   //\} --- end of jet structure -------------------------------------
 
 
 
   //----------------------------------------------------------------------
   /// @name Members mainly intended for internal use
   //----------------------------------------------------------------------
   //\{
   /// return the cluster_hist_index, intended to be used by clustering
   /// routines.
   inline int cluster_hist_index() const {return _cluster_hist_index;}
   /// set the cluster_hist_index, intended to be used by clustering routines.
   inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;}
 
   /// alternative name for cluster_hist_index() [perhaps more meaningful]
   inline int cluster_sequence_history_index() const {
     return cluster_hist_index();}
   /// alternative name for set_cluster_hist_index(...) [perhaps more
   /// meaningful]
   inline void set_cluster_sequence_history_index(const int index) {
     set_cluster_hist_index(index);}
 
   //\} ---- end of internal use functions ---------------------------
 
  protected:  
 
   SharedPtr<PseudoJetStructureBase> _structure;
   SharedPtr<UserInfoBase> _user_info;
 
 
  private: 
   // NB: following order must be kept for things to behave sensibly...
   double _px,_py,_pz,_E;
   mutable double _phi, _rap;
   double _kt2; 
   int    _cluster_hist_index, _user_index;
 
 #ifdef FASTJET_HAVE_THREAD_SAFETY
   enum {
     Init_Done=1,
     Init_NotDone=0,
     Init_InProgress=-1
   };
 
   mutable std::atomic<int> _init_status;
 #endif
   
   /// calculate phi, rap, kt2 based on the 4-momentum components
   void _finish_init();
   /// set the indices to default values
   void _reset_indices();
   
 //std::shared_ptr: /// reset the shared pointers to empty ones
 //std::shared_ptr: void _reset_shared_pointers();
 //std::shared_ptr: 
 //std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
 //std::shared_ptr:   /// For jets associated with a ClusterSequence, this will "free" the
 //std::shared_ptr:   /// jet from the ClusterSequence. This means that, for self-deleting
 //std::shared_ptr:   /// cluster sequences, it will reset the structure pointer and check
 //std::shared_ptr:   /// if the CS needs to be deleted.
 //std::shared_ptr:   ///
 //std::shared_ptr:   /// This is the replacement mechanism for self-deleting cs (with
 //std::shared_ptr:   /// set_count not working with std shared_ptr) and...
 //std::shared_ptr:   ///
 //std::shared_ptr:   /// IT HAS TO BE CALLED BEFORE ANY CHANGE OF THE JET STRUCTURE
 //std::shared_ptr:   /// POINTER
 //std::shared_ptr:   void _release_jet_from_cs();
 //std::shared_ptr: #endif //FASTJET_HAVE_THREAD_SAFETY
 
   /// ensure that the internal values for rapidity and phi 
   /// correspond to 4-momentum structure
 #ifdef FASTJET_HAVE_THREAD_SAFETY
   void _ensure_valid_rap_phi() const;
 #else
   inline void _ensure_valid_rap_phi() const{
     if (_phi == pseudojet_invalid_phi) _set_rap_phi();
   }
 #endif
   
   /// set cached rapidity and phi values
   void _set_rap_phi() const;
 
   // needed for operator* to have access to _ensure_valid_rap_phi()
   friend PseudoJet operator*(double, const PseudoJet &);
 };
 
 
 //----------------------------------------------------------------------
 // routines for basic binary operations
 
 PseudoJet operator+(const PseudoJet &, const PseudoJet &);
 PseudoJet operator-(const PseudoJet &, const PseudoJet &);
 PseudoJet operator*(double, const PseudoJet &);
 PseudoJet operator*(const PseudoJet &, double);
 PseudoJet operator/(const PseudoJet &, double);
 
 /// returns true if the 4 momentum components of the two PseudoJets
 /// are identical and all the internal indices (user, cluster_history)
 /// + structure and user-info shared pointers are too
 bool operator==(const PseudoJet &, const PseudoJet &);
 
 /// inequality test which is exact opposite of operator==
 inline bool operator!=(const PseudoJet & a, const PseudoJet & b) {return !(a==b);}
 
 /// Can only be used with val=0 and tests whether all four
 /// momentum components are equal to val (=0.0)
 bool operator==(const PseudoJet & jet, const double val);
 inline bool operator==(const double val, const PseudoJet & jet) {return jet == val;}
 
 /// Can only be used with val=0 and tests whether at least one of the
 /// four momentum components is different from val (=0.0)
 inline bool operator!=(const PseudoJet & a, const double val)  {return !(a==val);}
 inline bool operator!=( const double val, const PseudoJet & a) {return !(a==val);}
 
 /// returns the 4-vector dot product of a and b
 inline double dot_product(const PseudoJet & a, const PseudoJet & b) {
   return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
 }
 
 /// returns the cosine of the angle between a and b
 inline double cos_theta(const PseudoJet & a, const PseudoJet & b) {
   double dot_3d = a.px()*b.px() + a.py()*b.py() + a.pz()*b.pz();
   return std::min(1.0, std::max(-1.0, dot_3d/sqrt(a.modp2()*b.modp2())));
 }
 
 /// returns the angle between a and b
 inline double theta(const PseudoJet & a, const PseudoJet & b) {
   return acos(cos_theta(a,b));
 }
 
 /// returns true if the momenta of the two input jets are identical
 bool have_same_momentum(const PseudoJet &, const PseudoJet &);
 
 /// return a pseudojet with the given pt, y, phi and mass
 /// (phi should satisfy -2pi<phi<4pi)
 PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0);
 
 //----------------------------------------------------------------------
 // Routines to do with providing sorted arrays of vectors.
 
 /// return a vector of jets sorted into decreasing transverse momentum
 std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets);
 
 /// return a vector of jets sorted into increasing rapidity
 std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets);
 
 /// return a vector of jets sorted into decreasing energy
 std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets);
 
 /// return a vector of jets sorted into increasing pz
 std::vector<PseudoJet> sorted_by_pz(const std::vector<PseudoJet> & jets);
 
 //----------------------------------------------------------------------
 // some code to help sorting
 
 /// sort the indices so that values[indices[0->n-1]] is sorted
 /// into increasing order 
 void sort_indices(std::vector<int> & indices, 
           const std::vector<double> & values);
 
 /// given a vector of values with a one-to-one correspondence with the
 /// vector of objects, sort objects into an order such that the
 /// associated values would be in increasing order (but don't actually
 /// touch the values vector in the process).
 template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects, 
                           const std::vector<double> & values) {
   //assert(objects.size() == values.size());
   if (objects.size() != values.size()){
     throw Error("fastjet::objects_sorted_by_values(...): the size of the 'objects' vector must match the size of the 'values' vector");
   }
   
   // get a vector of indices
   std::vector<int> indices(values.size());
   for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
   
   // sort the indices
   sort_indices(indices, values);
   
   // copy the objects 
   std::vector<T> objects_sorted(objects.size());
   
   // place the objects in the correct order
   for (size_t i = 0; i < indices.size(); i++) {
     objects_sorted[i] = objects[indices[i]];
   }
 
   return objects_sorted;
 }
 
 /// \if internal_doc
 /// @ingroup internal
 /// \class IndexedSortHelper
 /// a class that helps us carry out indexed sorting.
 /// \endif
 class IndexedSortHelper {
 public:
   inline IndexedSortHelper (const std::vector<double> * reference_values) {
     _ref_values = reference_values;
   };
   inline int operator() (const int i1, const int i2) const {
     return  (*_ref_values)[i1] < (*_ref_values)[i2];
   };
 private:
   const std::vector<double> * _ref_values;
 };
 
 
 //----------------------------------------------------------------------
 /// constructor from any object that has px,py,pz,E = some_four_vector[0--3],
 // NB: do not know if it really needs to be inline, but when it wasn't
 //     linking failed with g++ (who knows what was wrong...)
 #ifndef SWIG
 template <class L> inline  PseudoJet::PseudoJet(const L & some_four_vector) {
   reset(some_four_vector);
 }
 #endif
 
 //----------------------------------------------------------------------
 inline void PseudoJet::_reset_indices() { 
   set_cluster_hist_index(-1);
   set_user_index(-1);
 
   _structure.reset();
   _user_info.reset();
 }
 
 //std::shared_ptr: inline void PseudoJet::_reset_shared_pointers() {
 //std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
 //std::shared_ptr:   // if the jet currently belongs to a cs, we need to release it before any chenge
 //std::shared_ptr:   _release_jet_from_cs();
 //std::shared_ptr: #endif // FASTJET_HAVE_THREAD_SAFETY
 //std::shared_ptr:   _structure.reset();
 //std::shared_ptr:   _user_info.reset();
 //std::shared_ptr:   // and do not forget to remove it in reset_indices above
 //std::shared_ptr: }
 
 
 
 
 // taken literally from CLHEP
 inline double PseudoJet::m() const {
   double mm = m2();
   return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
 }
 
 
 inline void PseudoJet::reset(double px_in, double py_in, double pz_in, double E_in) {
   _px = px_in;
   _py = py_in;
   _pz = pz_in;
   _E  = E_in;
   _finish_init();
   _reset_indices();
   //std::shared_ptr: _reset_shared_pointers();
 }
 
 inline void PseudoJet::reset_momentum(double px_in, double py_in, double pz_in, double E_in) {
   _px = px_in;
   _py = py_in;
   _pz = pz_in;
   _E  = E_in;
   _finish_init();
 }
 
 inline void PseudoJet::reset_momentum(const PseudoJet & pj) {
   _px  = pj._px ;
   _py  = pj._py ;
   _pz  = pj._pz ;
   _E   = pj._E  ;
 #ifdef FASTJET_HAVE_THREAD_SAFETY
   // the following lines are here to address the situation where
   // the current form of the jet has Init_Done, but the other jet does
   // not (or is in the process of being initialized, so that copying
   // is dangerous). 
   // We take the approach of verifying if the pj has been initialised
   // and if it has we take its cached phi & rapidity, otherwise
   // we call _finish_init(), which sets Init_NotDone and puts
   // rapidity & phi to invalid values. 
   int expected = Init_Done;
   if (pj._init_status.compare_exchange_weak(expected, Init_Done)) {
     _init_status = Init_Done;
     _phi = pj._phi;
     _rap = pj._rap;
     _kt2 = pj._kt2;
   } else {
     _finish_init();
   }
 #else 
   _phi = pj._phi;
   _rap = pj._rap;
   _kt2 = pj._kt2;
 #endif
 }
 
 //-------------------------------------------------------------------------------
 // implementation of the templated accesses to the underlying structyre
 //-------------------------------------------------------------------------------
 
 // returns a reference to the structure casted to the requested
 // structure type
 //
 // If there is no sructure associated, an Error is thrown.
 // If the type is not met, a std::bad_cast error is thrown.
 template<typename StructureType>
 const StructureType & PseudoJet::structure() const{
   return dynamic_cast<const StructureType &>(* validated_structure_ptr());
   
 }
 
 // check if the PseudoJet has the structure resulting from a Transformer 
 // (that is, its structure is compatible with a Transformer::StructureType)
 template<typename TransformerType>
 bool PseudoJet::has_structure_of() const{
   if (!_structure) return false;
 
   return dynamic_cast<const typename TransformerType::StructureType *>(_structure.get()) != 0;
 }
 
 // this is a helper to access a structure created by a Transformer 
 // (that is, of type Transformer::StructureType)
 // NULL is returned if the corresponding type is not met
 template<typename TransformerType>
 const typename TransformerType::StructureType & PseudoJet::structure_of() const{
   if (!_structure) 
     throw Error("Trying to access the structure of a PseudoJet without an associated structure");
 
   return dynamic_cast<const typename TransformerType::StructureType &>(*_structure);
 }
 
 
 
 //-------------------------------------------------------------------------------
 // helper functions to build a jet made of pieces
 //
 // Note that there are more complete versions of these functions, with
 // an additional argument for a recombination scheme, in
 // JetDefinition.hh
 // -------------------------------------------------------------------------------
 
 /// build a "CompositeJet" from the vector of its pieces
 ///
 /// In this case, E-scheme recombination is assumed to compute the
 /// total momentum
 PseudoJet join(const std::vector<PseudoJet> & pieces);
 
 /// build a MergedJet from a single PseudoJet
 PseudoJet join(const PseudoJet & j1);
 
 /// build a MergedJet from 2 PseudoJet
 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2);
 
 /// build a MergedJet from 3 PseudoJet
 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3);
 
 /// build a MergedJet from 4 PseudoJet
 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4);
 
 
 
 FASTJET_END_NAMESPACE
 
 #endif // __FASTJET_PSEUDOJET_HH__

This page was automatically generated by the 2.3.7 LXR engine. The LXR team