Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:57:18

0001 //FJSTARTHEADER
0002 // $Id$
0003 //
0004 // Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
0005 //
0006 //----------------------------------------------------------------------
0007 // This file is part of FastJet.
0008 //
0009 //  FastJet is free software; you can redistribute it and/or modify
0010 //  it under the terms of the GNU General Public License as published by
0011 //  the Free Software Foundation; either version 2 of the License, or
0012 //  (at your option) any later version.
0013 //
0014 //  The algorithms that underlie FastJet have required considerable
0015 //  development. They are described in the original FastJet paper,
0016 //  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
0017 //  FastJet as part of work towards a scientific publication, please
0018 //  quote the version you use and include a citation to the manual and
0019 //  optionally also to hep-ph/0512210.
0020 //
0021 //  FastJet is distributed in the hope that it will be useful,
0022 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
0023 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0024 //  GNU General Public License for more details.
0025 //
0026 //  You should have received a copy of the GNU General Public License
0027 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
0028 //----------------------------------------------------------------------
0029 //FJENDHEADER
0030 
0031 
0032 #ifndef __FASTJET_PSEUDOJET_HH__
0033 #define __FASTJET_PSEUDOJET_HH__
0034 
0035 #include<valarray>
0036 #include<vector>
0037 #include<cassert>
0038 #include<cmath>
0039 #include<iostream>
0040 #include "fastjet/config.h"
0041 #include "fastjet/internal/numconsts.hh"
0042 #include "fastjet/internal/IsBase.hh"
0043 #include "fastjet/SharedPtr.hh"
0044 #include "fastjet/Error.hh"
0045 #include "fastjet/PseudoJetStructureBase.hh"
0046 
0047 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0048 
0049 //using namespace std;
0050 
0051 /// Used to protect against parton-level events where pt can be zero
0052 /// for some partons, giving rapidity=infinity. KtJet fails in those cases.
0053 const double MaxRap = 1e5;
0054 
0055 /// default value for phi, meaning it (and rapidity) have yet to be calculated) 
0056 const double pseudojet_invalid_phi = -100.0;
0057 const double pseudojet_invalid_rap = -1e200;
0058 
0059 #ifndef __FJCORE__
0060 // forward definition
0061 class ClusterSequenceAreaBase;
0062 #endif  // __FJCORE__
0063 
0064 /// @ingroup basic_classes
0065 /// \class PseudoJet
0066 /// Class to contain pseudojets, including minimal information of use to
0067 /// jet-clustering routines.
0068 class PseudoJet {
0069 
0070  public:
0071   //----------------------------------------------------------------------
0072   /// @name Constructors and destructor
0073   //\{
0074   /// default constructor, which as of FJ3.0 provides an object for
0075   /// which all operations are now valid and which has zero momentum
0076   ///
0077   // (cf. this is actually OK from a timing point of view and in some
0078   // cases better than just having the default constructor for the
0079   // internal shared pointer: see PJtiming.cc and the notes therein)
0080   //
0081   // note: no reset of shared pointers needed
0082   PseudoJet() : _px(0), _py(0), _pz(0), _E(0) {_finish_init(); _reset_indices();}
0083   //PseudoJet() : _px(0), _py(0), _pz(0), _E(0), _phi(pseudojet_invalid_phi), _rap(pseudojet_invalid_rap), _kt2(0) {_reset_indices();}
0084   /// construct a pseudojet from explicit components
0085   PseudoJet(const double px, const double py, const double pz, const double E);
0086 
0087   /// constructor from any object that has px,py,pz,E = some_four_vector[0--3],
0088   #ifndef SWIG
0089   template <class L> PseudoJet(const L & some_four_vector);
0090   #endif
0091   
0092   // Constructor that performs minimal initialisation (only that of
0093   // the shared pointers), of use in certain speed-critical contexts
0094   //
0095   // NB: "dummy" is commented to avoid unused-variable compiler warnings
0096   PseudoJet(bool /* dummy */) {}
0097 
0098 #ifdef FASTJET_HAVE_THREAD_SAFETY
0099   PseudoJet(const PseudoJet &other){ (*this)=other; }
0100   PseudoJet& operator=(const PseudoJet& other);
0101 #endif
0102   
0103   /// default (virtual) destructor
0104   virtual ~PseudoJet(){}
0105 //std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
0106 //std::shared_ptr:     ;
0107 //std::shared_ptr: #else
0108 //std::shared_ptr:     {}
0109 //std::shared_ptr: #endif
0110   //\} ---- end of constructors and destructors --------------------------
0111 
0112   //----------------------------------------------------------------------
0113   /// @name Kinematic access functions
0114   //\{
0115   //----------------------------------------------------------------------
0116   inline double E()   const {return _E;}
0117   inline double e()   const {return _E;} // like CLHEP
0118   inline double px()  const {return _px;}
0119   inline double py()  const {return _py;}
0120   inline double pz()  const {return _pz;}
0121 
0122   /// returns phi (in the range 0..2pi)
0123   inline double phi() const {return phi_02pi();}
0124 
0125   /// returns phi in the range -pi..pi
0126   inline double phi_std()  const {
0127     _ensure_valid_rap_phi();
0128     return _phi > pi ? _phi-twopi : _phi;}
0129 
0130   /// returns phi in the range 0..2pi
0131   inline double phi_02pi() const {
0132     _ensure_valid_rap_phi();
0133     return _phi;
0134   }
0135 
0136   /// returns the rapidity or some large value when the rapidity
0137   /// is infinite
0138   inline double rap() const {
0139     _ensure_valid_rap_phi();
0140     return _rap;
0141   }
0142 
0143   /// the same as rap()
0144   inline double rapidity() const {return rap();} // like CLHEP
0145 
0146   /// returns the pseudo-rapidity or some large value when the
0147   /// rapidity is infinite
0148   double pseudorapidity() const;
0149   double eta() const {return pseudorapidity();}
0150 
0151   /// returns the squared transverse momentum
0152   inline double pt2() const {return _kt2;}
0153   /// returns the scalar transverse momentum
0154   inline double  pt() const {return sqrt(_kt2);} 
0155   /// returns the squared transverse momentum
0156   inline double perp2() const {return _kt2;}  // like CLHEP
0157   /// returns the scalar transverse momentum
0158   inline double  perp() const {return sqrt(_kt2);}    // like CLHEP
0159   /// returns the squared transverse momentum
0160   inline double kt2() const {return _kt2;} // for bkwds compatibility
0161 
0162   /// returns the squared invariant mass // like CLHEP
0163   inline double  m2() const {return (_E+_pz)*(_E-_pz)-_kt2;}    
0164   /// returns the invariant mass 
0165   /// (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP)
0166   inline double  m() const;    
0167 
0168   /// returns the squared transverse mass = kt^2+m^2
0169   inline double mperp2() const {return (_E+_pz)*(_E-_pz);}
0170   /// returns the transverse mass = sqrt(kt^2+m^2)
0171   inline double mperp() const {return sqrt(std::abs(mperp2()));}
0172   /// returns the squared transverse mass = kt^2+m^2
0173   inline double mt2() const {return (_E+_pz)*(_E-_pz);}
0174   /// returns the transverse mass = sqrt(kt^2+m^2)
0175   inline double mt() const {return sqrt(std::abs(mperp2()));}
0176 
0177   /// return the squared 3-vector modulus = px^2+py^2+pz^2
0178   inline double modp2() const {return _kt2+_pz*_pz;}
0179   /// return the 3-vector modulus = sqrt(px^2+py^2+pz^2)
0180   inline double modp() const {return sqrt(_kt2+_pz*_pz);}
0181 
0182   /// return the transverse energy
0183   inline double Et() const {return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);}
0184   /// return the transverse energy squared
0185   inline double Et2() const {return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);}
0186 
0187   /// cos of the polar angle
0188   /// should we have: min(1.0,max(-1.0,_pz/sqrt(modp2()))); 
0189   inline double cos_theta() const {
0190     return std::min(1.0, std::max(-1.0, _pz/sqrt(modp2())));
0191   }
0192   /// polar angle
0193   inline double theta() const { return acos(cos_theta()); }
0194 
0195   /// returns component i, where X==0, Y==1, Z==2, E==3
0196   double operator () (int i) const ; 
0197   /// returns component i, where X==0, Y==1, Z==2, E==3
0198   inline double operator [] (int i) const { return (*this)(i); }; // this too
0199 
0200 
0201 
0202   /// returns kt distance (R=1) between this jet and another
0203   double kt_distance(const PseudoJet & other) const;
0204 
0205   /// returns squared cylinder (rap-phi) distance between this jet and another
0206   double plain_distance(const PseudoJet & other) const;
0207   /// returns squared cylinder (rap-phi) distance between this jet and
0208   /// another
0209   inline double squared_distance(const PseudoJet & other) const {
0210     return plain_distance(other);}
0211 
0212   /// return the cylinder (rap-phi) distance between this jet and another,
0213   /// \f$\Delta_R = \sqrt{\Delta y^2 + \Delta \phi^2}\f$.
0214   inline double delta_R(const PseudoJet & other) const {
0215     return sqrt(squared_distance(other));
0216   }
0217 
0218   /// returns other.phi() - this.phi(), constrained to be in 
0219   /// range -pi .. pi
0220   double delta_phi_to(const PseudoJet & other) const;
0221 
0222   //// this seemed to compile except if it was used
0223   //friend inline double 
0224   //  kt_distance(const PseudoJet & jet1, const PseudoJet & jet2) { 
0225   //                                      return jet1.kt_distance(jet2);}
0226 
0227   /// returns distance between this jet and the beam
0228   inline double beam_distance() const {return _kt2;}
0229 
0230   /// return a valarray containing the four-momentum (components 0-2
0231   /// are 3-mom, component 3 is energy).
0232   std::valarray<double> four_mom() const;
0233 
0234   //\}  ------- end of kinematic access functions
0235 
0236   // taken from CLHEP
0237   enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
0238 
0239 
0240   //----------------------------------------------------------------------
0241   /// @name Kinematic modification functions
0242   //\{
0243   //----------------------------------------------------------------------
0244   /// transform this jet (given in the rest frame of prest) into a jet
0245   /// in the lab frame
0246   PseudoJet & boost(const PseudoJet & prest);
0247   /// transform this jet (given in lab) into a jet in the rest
0248   /// frame of prest
0249   PseudoJet & unboost(const PseudoJet & prest);
0250 
0251   PseudoJet & operator*=(double);
0252   PseudoJet & operator/=(double);
0253   PseudoJet & operator+=(const PseudoJet &);
0254   PseudoJet & operator-=(const PseudoJet &);
0255 
0256 //std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
0257 //std::shared_ptr:   /// overload the assignment through the = operator
0258 //std::shared_ptr:   ///
0259 //std::shared_ptr:   /// this is needed to make sure that release_from_cs is called
0260 //std::shared_ptr:   /// before copying the structure shared pointer!
0261 //std::shared_ptr:   void operator=(const PseudoJet &other){
0262 //std::shared_ptr:     _release_jet_from_cs();
0263 //std::shared_ptr:     _structure = other._structure;
0264 //std::shared_ptr:     _user_info = other._user_info;
0265 //std::shared_ptr: 
0266 //std::shared_ptr:     _px = other._px;
0267 //std::shared_ptr:     _py = other._py;
0268 //std::shared_ptr:     _pz = other._pz;
0269 //std::shared_ptr:     _E  = other._E;
0270 //std::shared_ptr: 
0271 //std::shared_ptr:     _phi = other._phi;
0272 //std::shared_ptr:     _rap = other._rap;
0273 //std::shared_ptr:     _kt2  = other._kt2;
0274 //std::shared_ptr:     
0275 //std::shared_ptr:     _cluster_hist_index = other._cluster_hist_index;
0276 //std::shared_ptr:     _user_index = other._user_index;
0277 //std::shared_ptr:   }
0278 //std::shared_ptr: 
0279 //std::shared_ptr:   /// force resetting the structure pointer to an empty structure
0280 //std::shared_ptr:   ///
0281 //std::shared_ptr:   /// THIS IS STRICTLY MEANT FOR INTERNAL USAGE IN SELF-DELETING
0282 //std::shared_ptr:   /// ClusterSequences. IT SHOULD NOT BE USED BY END-USERS.
0283 //std::shared_ptr:   void force_reset_structure(){
0284 //std::shared_ptr:     _structure.reset();
0285 //std::shared_ptr:   }  
0286 //std::shared_ptr: #endif
0287 
0288   /// reset the 4-momentum according to the supplied components and
0289   /// put the user and history indices back to their default values
0290   inline void reset(double px, double py, double pz, double E);
0291 
0292   /// reset the PseudoJet to be equal to psjet (including its
0293   /// indices); NB if the argument is derived from a PseudoJet then
0294   /// the "reset" used will be the templated version
0295   ///
0296   /// Note: this is included on top of the templated version because
0297   /// PseudoJet is not "derived" from PseudoJet, so the templated
0298   /// reset would not handle this case properly.
0299   inline void reset(const PseudoJet & psjet) {
0300     (*this) = psjet;
0301   }
0302 
0303   /// reset the 4-momentum according to the supplied generic 4-vector
0304   /// (accessible via indexing, [0]==px,...[3]==E) and put the user
0305   /// and history indices back to their default values.
0306 #ifndef SWIG
0307   template <class L> inline void reset(const L & some_four_vector) {
0308     // check if some_four_vector can be cast to a PseudoJet
0309     //
0310     // Note that a regular dynamic_cast would not work here because
0311     // there is no guarantee that L is polymorphic. We use a more
0312     // complex construct here that works also in such a case. As for
0313     // dynamic_cast, NULL is returned if L is not derived from
0314     // PseudoJet
0315     //
0316     // Note the explicit request for fastjet::cast_if_derived; when
0317     // combining fastjet and fjcore, this avoids ambiguity in which of
0318     // the two cast_if_derived calls to use.
0319     const PseudoJet * pj = fastjet::cast_if_derived<const PseudoJet>(&some_four_vector);
0320 
0321     if (pj){
0322       (*this) = *pj;
0323     } else {
0324       reset(some_four_vector[0], some_four_vector[1],
0325         some_four_vector[2], some_four_vector[3]);
0326     }
0327   }
0328 #endif // SWIG
0329 
0330   /// reset the PseudoJet according to the specified pt, rapidity,
0331   /// azimuth and mass (also resetting indices, etc.)
0332   /// (phi should satisfy -2pi<phi<4pi)
0333   inline void reset_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in=0.0) {
0334     reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in);
0335     _reset_indices();
0336     //std::shared_ptr: _reset_shared_pointers();
0337   }
0338 
0339   /// reset the 4-momentum according to the supplied components 
0340   /// but leave all other information (indices, user info, etc.)
0341   /// untouched
0342   inline void reset_momentum(double px, double py, double pz, double E);
0343 
0344   /// reset the 4-momentum according to the components of the supplied
0345   /// PseudoJet, including cached components; note that the template
0346   /// version (below) will be called for classes derived from PJ.
0347   inline void reset_momentum(const PseudoJet & pj);
0348 
0349   /// reset the 4-momentum according to the specified pt, rapidity,
0350   /// azimuth and mass (phi should satisfy -2pi<phi<4pi)
0351   void reset_momentum_PtYPhiM(double pt, double y, double phi, double m=0.0);
0352 
0353   /// reset the 4-momentum according to the supplied generic 4-vector
0354   /// (accessible via indexing, [0]==px,...[3]==E), but leave all
0355   /// other information (indices, user info, etc.)  untouched
0356   template <class L> inline void reset_momentum(const L & some_four_vector) {
0357     reset_momentum(some_four_vector[0], some_four_vector[1],
0358            some_four_vector[2], some_four_vector[3]);
0359   }
0360 
0361   /// in some cases when setting a 4-momentum, the user/program knows
0362   /// what rapidity and azimuth are associated with that 4-momentum;
0363   /// by calling this routine the user can provide the information
0364   /// directly to the PseudoJet and avoid expensive rap-phi
0365   /// recalculations.
0366   ///
0367   /// - \param rap  rapidity
0368   /// - \param phi  (in range -twopi...4*pi)
0369   ///
0370   /// USE WITH CAUTION: there are no checks that the rapidity and
0371   /// azimuth supplied are sensible, nor does this reset the
0372   /// 4-momentum components if things don't match.
0373   void set_cached_rap_phi(double rap, double phi);
0374 
0375 
0376   //\} --- end of kin mod functions ------------------------------------
0377 
0378   //----------------------------------------------------------------------
0379   /// @name User index functions
0380   ///
0381   /// To allow the user to set and access an integer index which can
0382   /// be exploited by the user to associate extra information with a
0383   /// particle/jet (for example pdg id, or an indication of a
0384   /// particle's origin within the user's analysis)
0385   //
0386   //\{
0387 
0388   /// return the user_index, 
0389   inline int user_index() const {return _user_index;}
0390   /// set the user_index, intended to allow the user to add simple
0391   /// identifying information to a particle/jet
0392   inline void set_user_index(const int index) {_user_index = index;}
0393 
0394   //\} ----- end of use index functions ---------------------------------
0395 
0396   //----------------------------------------------------------------------
0397   /// @name User information types and functions
0398   ///
0399   /// Allows PseudoJet to carry extra user info (as an object derived from
0400   /// UserInfoBase).
0401   //\{
0402 
0403   /// @ingroup user_info
0404   /// \class UserInfoBase
0405   /// a base class to hold extra user information in a PseudoJet
0406   ///
0407   /// This is a base class to help associate extra user information
0408   /// with a jet. The user should store their information in a class
0409   /// derived from this. This allows information of arbitrary
0410   /// complexity to be easily associated with a PseudoJet (in contrast
0411   /// to the user index). For example, in a Monte Carlo simulation,
0412   /// the user information might include the PDG ID, and the position
0413   /// of the production vertex for the particle.
0414   ///
0415   /// The PseudoJet is able to store a shared pointer to any object
0416   /// derived from UserInfo. The use of a shared pointer frees the
0417   /// user of the need to handle the memory management associated with
0418   /// the information.
0419   ///
0420   /// Having the user information derive from a common base class also
0421   /// facilitates dynamic casting, etc.
0422   ///
0423   class UserInfoBase{
0424   public:
0425     // dummy ctor
0426     UserInfoBase(){};
0427 
0428     // dummy virtual dtor
0429     // makes it polymorphic to allow for dynamic_cast
0430     virtual ~UserInfoBase(){}; 
0431   };
0432 
0433   /// error class to be thrown if accessing user info when it doesn't
0434   /// exist
0435   class InexistentUserInfo : public Error {
0436   public:
0437     InexistentUserInfo();
0438   };
0439 
0440   /// sets the internal shared pointer to the user information.
0441   ///
0442   /// Note that the PseudoJet will now _own_ the pointer, and delete
0443   /// the corresponding object when it (the jet, and any copies of the jet)
0444   /// goes out of scope. 
0445   void set_user_info(UserInfoBase * user_info_in) {
0446     _user_info.reset(user_info_in);
0447   }
0448 
0449   /// returns a reference to the dynamic cast conversion of user_info
0450   /// to type L.
0451   ///
0452   /// Usage: suppose you have previously set the user info with a pointer
0453   /// to an object of type MyInfo, 
0454   ///
0455   ///   class MyInfo: public PseudoJet::UserInfoBase {
0456   ///      MyInfo(int id) : _pdg_id(id);
0457   ///      int pdg_id() const {return _pdg_id;}
0458   ///      int _pdg_id;
0459   ///   };
0460   ///
0461   ///   PseudoJet particle(...);
0462   ///   particle.set_user_info(new MyInfo(its_pdg_id));
0463   ///
0464   /// Then you would access that pdg_id() as
0465   ///
0466   ///   particle.user_info<MyInfo>().pdg_id();
0467   ///
0468   /// It's overkill for just a single integer, but scales easily to
0469   /// more extensive information.
0470   ///
0471   /// Note that user_info() throws an InexistentUserInfo() error if
0472   /// there is no user info; throws a std::bad_cast if the conversion
0473   /// doesn't work
0474   ///
0475   /// If this behaviour does not fit your needs, use instead the the
0476   /// user_info_ptr() or user_info_shared_ptr() member functions.
0477   template<class L>
0478   const L & user_info() const{
0479     if (_user_info.get() == 0) throw InexistentUserInfo();
0480     return dynamic_cast<const L &>(* _user_info.get());
0481   }
0482 
0483   /// returns true if the PseudoJet has user information
0484   bool has_user_info() const{
0485     return _user_info.get();
0486   }
0487 
0488   /// returns true if the PseudoJet has user information than can be
0489   /// cast to the template argument type.
0490   template<class L>
0491   bool has_user_info() const{
0492     return _user_info.get() && dynamic_cast<const L *>(_user_info.get());
0493   }
0494 
0495   /// retrieve a pointer to the (const) user information
0496   const UserInfoBase * user_info_ptr() const{
0497     // the line below is not needed since the next line would anyway
0498     // return NULL in that case
0499     //if (!_user_info) return NULL;
0500     return _user_info.get();
0501   }
0502 
0503 
0504   /// retrieve a (const) shared pointer to the user information
0505   const SharedPtr<UserInfoBase> & user_info_shared_ptr() const{
0506     return _user_info;
0507   }
0508 
0509   /// retrieve a (non-const) shared pointer to the user information;
0510   /// you can use this, for example, to set the shared pointer, eg
0511   ///
0512   /// \code
0513   ///   p2.user_info_shared_ptr() = p1.user_info_shared_ptr();
0514   /// \endcode
0515   ///
0516   /// or 
0517   ///
0518   /// \code
0519   ///   SharedPtr<PseudoJet::UserInfoBase> info_shared(new MyInfo(...));
0520   ///   p2.user_info_shared_ptr() = info_shared;
0521   /// \endcode
0522   SharedPtr<UserInfoBase> & user_info_shared_ptr(){
0523     return _user_info;
0524   }
0525 
0526   // \} --- end of extra info functions ---------------------------------
0527 
0528   //----------------------------------------------------------------------
0529   /// @name Description
0530   ///
0531   /// Since a PseudoJet can have a structure that contains a variety
0532   /// of information, we provide a description that allows one to check
0533   /// exactly what kind of PseudoJet we are dealing with
0534   //
0535   //\{
0536 
0537   /// return a string describing what kind of PseudoJet we are dealing with 
0538   std::string description() const;
0539 
0540   //\} ----- end of description functions ---------------------------------
0541 
0542   //-------------------------------------------------------------
0543   /// @name Access to the associated ClusterSequence object.
0544   ///
0545   /// In addition to having kinematic information, jets may contain a
0546   /// reference to an associated ClusterSequence (this is the case,
0547   /// for example, if the jet has been returned by a ClusterSequence
0548   /// member function).
0549   //\{
0550   //-------------------------------------------------------------
0551   /// returns true if this PseudoJet has an associated ClusterSequence.
0552   bool has_associated_cluster_sequence() const;
0553   /// shorthand for has_associated_cluster_sequence()
0554   bool has_associated_cs() const {return has_associated_cluster_sequence();}
0555 
0556   /// returns true if this PseudoJet has an associated and still
0557   /// valid(ated) ClusterSequence.
0558   bool has_valid_cluster_sequence() const;
0559   /// shorthand for has_valid_cluster_sequence()
0560   bool has_valid_cs() const {return has_valid_cluster_sequence();}
0561 
0562   /// get a (const) pointer to the parent ClusterSequence (NULL if
0563   /// inexistent)
0564   const ClusterSequence* associated_cluster_sequence() const;
0565   // shorthand for associated_cluster_sequence()
0566   const ClusterSequence* associated_cs() const {return associated_cluster_sequence();}
0567 
0568   /// if the jet has a valid associated cluster sequence then return a
0569   /// pointer to it; otherwise throw an error
0570   inline const ClusterSequence * validated_cluster_sequence() const {
0571     return validated_cs();
0572   }
0573   /// shorthand for validated_cluster_sequence()
0574   const ClusterSequence * validated_cs() const;
0575 
0576 #ifndef __FJCORE__
0577   /// if the jet has valid area information then return a pointer to
0578   /// the associated ClusterSequenceAreaBase object; otherwise throw an error
0579   inline const ClusterSequenceAreaBase * validated_cluster_sequence_area_base() const {
0580     return validated_csab();
0581   }
0582 
0583   /// shorthand for validated_cluster_sequence_area_base()
0584   const ClusterSequenceAreaBase * validated_csab() const;
0585 #endif  //  __FJCORE__
0586 
0587   //\}
0588 
0589   //-------------------------------------------------------------
0590   /// @name Access to the associated PseudoJetStructureBase object.
0591   ///
0592   /// In addition to having kinematic information, jets may contain a
0593   /// reference to an associated ClusterSequence (this is the case,
0594   /// for example, if the jet has been returned by a ClusterSequence
0595   /// member function).
0596   //\{
0597   //-------------------------------------------------------------
0598 
0599   /// set the associated structure
0600   void set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure_in);
0601 
0602   /// return true if there is some structure associated with this PseudoJet
0603   bool has_structure() const;
0604 
0605   /// return a pointer to the structure (of type
0606   /// PseudoJetStructureBase*) associated with this PseudoJet.
0607   ///
0608   /// return NULL if there is no associated structure
0609   const PseudoJetStructureBase* structure_ptr() const;
0610   
0611   /// return a non-const pointer to the structure (of type
0612   /// PseudoJetStructureBase*) associated with this PseudoJet.
0613   ///
0614   /// return NULL if there is no associated structure
0615   ///
0616   /// Only use this if you know what you are doing. In any case,
0617   /// prefer the 'structure_ptr()' (the const version) to this method,
0618   /// unless you really need a write access to the PseudoJet's
0619   /// underlying structure.
0620   PseudoJetStructureBase* structure_non_const_ptr();
0621   
0622   /// return a pointer to the structure (of type
0623   /// PseudoJetStructureBase*) associated with this PseudoJet.
0624   ///
0625   /// throw an error if there is no associated structure
0626   const PseudoJetStructureBase* validated_structure_ptr() const;
0627   
0628   /// return a reference to the shared pointer to the
0629   /// PseudoJetStructureBase associated with this PseudoJet
0630   const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const;
0631 
0632   /// returns a reference to the structure casted to the requested
0633   /// structure type
0634   ///
0635   /// If there is no structure associated, an Error is thrown.
0636   /// If the type is not met, a std::bad_cast error is thrown.
0637   template<typename StructureType>
0638   const StructureType & structure() const;
0639 
0640   /// check if the PseudoJet has the structure resulting from a Transformer 
0641   /// (that is, its structure is compatible with a Transformer::StructureType).
0642   /// If there is no structure, false is returned.
0643   template<typename TransformerType>
0644   bool has_structure_of() const;
0645 
0646   /// this is a helper to access any structure created by a Transformer 
0647   /// (that is, of type Transformer::StructureType).
0648   ///
0649   /// If there is no structure, or if the structure is not compatible
0650   /// with TransformerType, an error is thrown.
0651   template<typename TransformerType>
0652   const typename TransformerType::StructureType & structure_of() const;
0653 
0654   //\}
0655 
0656   //-------------------------------------------------------------
0657   /// @name Methods for access to information about jet structure
0658   ///
0659   /// These allow access to jet constituents, and other jet
0660   /// subtructure information. They only work if the jet is associated
0661   /// with a ClusterSequence.
0662   //-------------------------------------------------------------
0663   //\{
0664 
0665   /// check if it has been recombined with another PseudoJet in which
0666   /// case, return its partner through the argument. Otherwise,
0667   /// 'partner' is set to 0.
0668   ///
0669   /// an Error is thrown if this PseudoJet has no currently valid
0670   /// associated ClusterSequence
0671   virtual bool has_partner(PseudoJet &partner) const;
0672 
0673   /// check if it has been recombined with another PseudoJet in which
0674   /// case, return its child through the argument. Otherwise, 'child'
0675   /// is set to 0.
0676   /// 
0677   /// an Error is thrown if this PseudoJet has no currently valid
0678   /// associated ClusterSequence
0679   virtual bool has_child(PseudoJet &child) const;
0680 
0681   /// check if it is the product of a recombination, in which case
0682   /// return the 2 parents through the 'parent1' and 'parent2'
0683   /// arguments. Otherwise, set these to 0.
0684   ///
0685   /// an Error is thrown if this PseudoJet has no currently valid
0686   /// associated ClusterSequence
0687   virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const;
0688 
0689   /// check if the current PseudoJet contains the one passed as
0690   /// argument.
0691   ///
0692   /// an Error is thrown if this PseudoJet has no currently valid
0693   /// associated ClusterSequence
0694   virtual bool contains(const PseudoJet &constituent) const;
0695 
0696   /// check if the current PseudoJet is contained the one passed as
0697   /// argument.
0698   ///
0699   /// an Error is thrown if this PseudoJet has no currently valid
0700   /// associated ClusterSequence
0701   virtual bool is_inside(const PseudoJet &jet) const;
0702 
0703 
0704   /// returns true if the PseudoJet has constituents
0705   virtual bool has_constituents() const;
0706 
0707   /// retrieve the constituents. 
0708   ///
0709   /// an Error is thrown if this PseudoJet has no currently valid
0710   /// associated ClusterSequence or other substructure information
0711   virtual std::vector<PseudoJet> constituents() const;
0712 
0713 
0714   /// returns true if the PseudoJet has support for exclusive subjets
0715   virtual bool has_exclusive_subjets() const;
0716 
0717   /// return a vector of all subjets of the current jet (in the sense
0718   /// of the exclusive algorithm) that would be obtained when running
0719   /// the algorithm with the given dcut. 
0720   ///
0721   /// Time taken is O(m ln m), where m is the number of subjets that
0722   /// are found. If m gets to be of order of the total number of
0723   /// constituents in the jet, this could be substantially slower than
0724   /// just getting that list of constituents.
0725   ///
0726   /// an Error is thrown if this PseudoJet has no currently valid
0727   /// associated ClusterSequence
0728   std::vector<PseudoJet> exclusive_subjets (const double dcut) const;
0729 
0730   /// return the size of exclusive_subjets(...); still n ln n with same
0731   /// coefficient, but marginally more efficient than manually taking
0732   /// exclusive_subjets.size()
0733   ///
0734   /// an Error is thrown if this PseudoJet has no currently valid
0735   /// associated ClusterSequence
0736   int n_exclusive_subjets(const double dcut) const;
0737 
0738   /// return the list of subjets obtained by unclustering the supplied
0739   /// jet down to nsub subjets. Throws an error if there are fewer than
0740   /// nsub particles in the jet.
0741   ///
0742   /// For ClusterSequence type jets, requires nsub ln nsub time
0743   ///
0744   /// An Error is thrown if this PseudoJet has no currently valid
0745   /// associated ClusterSequence
0746   std::vector<PseudoJet> exclusive_subjets (int nsub) const;
0747 
0748   /// return the list of subjets obtained by unclustering the supplied
0749   /// jet down to nsub subjets (or all constituents if there are fewer
0750   /// than nsub).
0751   ///
0752   /// For ClusterSequence type jets, requires nsub ln nsub time
0753   ///
0754   /// An Error is thrown if this PseudoJet has no currently valid
0755   /// associated ClusterSequence
0756   std::vector<PseudoJet> exclusive_subjets_up_to (int nsub) const;
0757 
0758   /// Returns the dij that was present in the merging nsub+1 -> nsub 
0759   /// subjets inside this jet.
0760   ///
0761   /// Returns 0 if there were nsub or fewer constituents in the jet.
0762   ///
0763   /// an Error is thrown if this PseudoJet has no currently valid
0764   /// associated ClusterSequence
0765   double exclusive_subdmerge(int nsub) const;
0766 
0767   /// Returns the maximum dij that occurred in the whole event at the
0768   /// stage that the nsub+1 -> nsub merge of subjets occurred inside 
0769   /// this jet.
0770   ///
0771   /// Returns 0 if there were nsub or fewer constituents in the jet.
0772   ///
0773   /// an Error is thrown if this PseudoJet has no currently valid
0774   /// associated ClusterSequence
0775   double exclusive_subdmerge_max(int nsub) const;
0776 
0777 
0778   /// returns true if a jet has pieces
0779   ///
0780   /// By default a single particle or a jet coming from a
0781   /// ClusterSequence have no pieces and this methos will return false.
0782   ///
0783   /// In practice, this is equivalent to have an structure of type
0784   /// CompositeJetStructure.
0785   virtual bool has_pieces() const;
0786 
0787 
0788   /// retrieve the pieces that make up the jet. 
0789   ///
0790   /// If the jet does not support pieces, an error is throw
0791   virtual std::vector<PseudoJet> pieces() const;
0792 
0793 
0794   // the following ones require a computation of the area in the
0795   // parent ClusterSequence (See ClusterSequenceAreaBase for details)
0796   //------------------------------------------------------------------
0797 #ifndef __FJCORE__
0798 
0799   /// check if it has a defined area
0800   virtual bool has_area() const;
0801 
0802   /// return the jet (scalar) area.
0803   /// throws an Error if there is no support for area in the parent CS
0804   virtual double area() const;
0805 
0806   /// return the error (uncertainty) associated with the determination
0807   /// of the area of this jet.
0808   /// throws an Error if there is no support for area in the parent CS
0809   virtual double area_error() const;
0810 
0811   /// return the jet 4-vector area.
0812   /// throws an Error if there is no support for area in the parent CS
0813   virtual PseudoJet area_4vector() const;
0814 
0815   /// true if this jet is made exclusively of ghosts.
0816   /// throws an Error if there is no support for area in the parent CS
0817   virtual bool is_pure_ghost() const;
0818 
0819 #endif  // __FJCORE__
0820   //\} --- end of jet structure -------------------------------------
0821 
0822 
0823 
0824   //----------------------------------------------------------------------
0825   /// @name Members mainly intended for internal use
0826   //----------------------------------------------------------------------
0827   //\{
0828   /// return the cluster_hist_index, intended to be used by clustering
0829   /// routines.
0830   inline int cluster_hist_index() const {return _cluster_hist_index;}
0831   /// set the cluster_hist_index, intended to be used by clustering routines.
0832   inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;}
0833 
0834   /// alternative name for cluster_hist_index() [perhaps more meaningful]
0835   inline int cluster_sequence_history_index() const {
0836     return cluster_hist_index();}
0837   /// alternative name for set_cluster_hist_index(...) [perhaps more
0838   /// meaningful]
0839   inline void set_cluster_sequence_history_index(const int index) {
0840     set_cluster_hist_index(index);}
0841 
0842   //\} ---- end of internal use functions ---------------------------
0843 
0844  protected:  
0845 
0846   SharedPtr<PseudoJetStructureBase> _structure;
0847   SharedPtr<UserInfoBase> _user_info;
0848 
0849 
0850  private: 
0851   // NB: following order must be kept for things to behave sensibly...
0852   double _px,_py,_pz,_E;
0853   mutable double _phi, _rap;
0854   double _kt2; 
0855   int    _cluster_hist_index, _user_index;
0856 
0857 #ifdef FASTJET_HAVE_THREAD_SAFETY
0858   enum {
0859     Init_Done=1,
0860     Init_NotDone=0,
0861     Init_InProgress=-1
0862   };
0863 
0864   mutable std::atomic<int> _init_status;
0865 #endif
0866   
0867   /// calculate phi, rap, kt2 based on the 4-momentum components
0868   void _finish_init();
0869   /// set the indices to default values
0870   void _reset_indices();
0871   
0872 //std::shared_ptr: /// reset the shared pointers to empty ones
0873 //std::shared_ptr: void _reset_shared_pointers();
0874 //std::shared_ptr: 
0875 //std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
0876 //std::shared_ptr:   /// For jets associated with a ClusterSequence, this will "free" the
0877 //std::shared_ptr:   /// jet from the ClusterSequence. This means that, for self-deleting
0878 //std::shared_ptr:   /// cluster sequences, it will reset the structure pointer and check
0879 //std::shared_ptr:   /// if the CS needs to be deleted.
0880 //std::shared_ptr:   ///
0881 //std::shared_ptr:   /// This is the replacement mechanism for self-deleting cs (with
0882 //std::shared_ptr:   /// set_count not working with std shared_ptr) and...
0883 //std::shared_ptr:   ///
0884 //std::shared_ptr:   /// IT HAS TO BE CALLED BEFORE ANY CHANGE OF THE JET STRUCTURE
0885 //std::shared_ptr:   /// POINTER
0886 //std::shared_ptr:   void _release_jet_from_cs();
0887 //std::shared_ptr: #endif //FASTJET_HAVE_THREAD_SAFETY
0888 
0889   /// ensure that the internal values for rapidity and phi 
0890   /// correspond to 4-momentum structure
0891 #ifdef FASTJET_HAVE_THREAD_SAFETY
0892   void _ensure_valid_rap_phi() const;
0893 #else
0894   inline void _ensure_valid_rap_phi() const{
0895     if (_phi == pseudojet_invalid_phi) _set_rap_phi();
0896   }
0897 #endif
0898   
0899   /// set cached rapidity and phi values
0900   void _set_rap_phi() const;
0901 
0902   // needed for operator* to have access to _ensure_valid_rap_phi()
0903   friend PseudoJet operator*(double, const PseudoJet &);
0904 };
0905 
0906 
0907 //----------------------------------------------------------------------
0908 // routines for basic binary operations
0909 
0910 PseudoJet operator+(const PseudoJet &, const PseudoJet &);
0911 PseudoJet operator-(const PseudoJet &, const PseudoJet &);
0912 PseudoJet operator*(double, const PseudoJet &);
0913 PseudoJet operator*(const PseudoJet &, double);
0914 PseudoJet operator/(const PseudoJet &, double);
0915 
0916 /// returns true if the 4 momentum components of the two PseudoJets
0917 /// are identical and all the internal indices (user, cluster_history)
0918 /// + structure and user-info shared pointers are too
0919 bool operator==(const PseudoJet &, const PseudoJet &);
0920 
0921 /// inequality test which is exact opposite of operator==
0922 inline bool operator!=(const PseudoJet & a, const PseudoJet & b) {return !(a==b);}
0923 
0924 /// Can only be used with val=0 and tests whether all four
0925 /// momentum components are equal to val (=0.0)
0926 bool operator==(const PseudoJet & jet, const double val);
0927 inline bool operator==(const double val, const PseudoJet & jet) {return jet == val;}
0928 
0929 /// Can only be used with val=0 and tests whether at least one of the
0930 /// four momentum components is different from val (=0.0)
0931 inline bool operator!=(const PseudoJet & a, const double val)  {return !(a==val);}
0932 inline bool operator!=( const double val, const PseudoJet & a) {return !(a==val);}
0933 
0934 /// returns the 4-vector dot product of a and b
0935 inline double dot_product(const PseudoJet & a, const PseudoJet & b) {
0936   return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
0937 }
0938 
0939 /// returns the cosine of the angle between a and b
0940 inline double cos_theta(const PseudoJet & a, const PseudoJet & b) {
0941   double dot_3d = a.px()*b.px() + a.py()*b.py() + a.pz()*b.pz();
0942   return std::min(1.0, std::max(-1.0, dot_3d/sqrt(a.modp2()*b.modp2())));
0943 }
0944 
0945 /// returns the angle between a and b
0946 inline double theta(const PseudoJet & a, const PseudoJet & b) {
0947   return acos(cos_theta(a,b));
0948 }
0949 
0950 /// returns true if the momenta of the two input jets are identical
0951 bool have_same_momentum(const PseudoJet &, const PseudoJet &);
0952 
0953 /// return a pseudojet with the given pt, y, phi and mass
0954 /// (phi should satisfy -2pi<phi<4pi)
0955 PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0);
0956 
0957 //----------------------------------------------------------------------
0958 // Routines to do with providing sorted arrays of vectors.
0959 
0960 /// return a vector of jets sorted into decreasing transverse momentum
0961 std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets);
0962 
0963 /// return a vector of jets sorted into increasing rapidity
0964 std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets);
0965 
0966 /// return a vector of jets sorted into decreasing energy
0967 std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets);
0968 
0969 /// return a vector of jets sorted into increasing pz
0970 std::vector<PseudoJet> sorted_by_pz(const std::vector<PseudoJet> & jets);
0971 
0972 //----------------------------------------------------------------------
0973 // some code to help sorting
0974 
0975 /// sort the indices so that values[indices[0->n-1]] is sorted
0976 /// into increasing order 
0977 void sort_indices(std::vector<int> & indices, 
0978           const std::vector<double> & values);
0979 
0980 /// given a vector of values with a one-to-one correspondence with the
0981 /// vector of objects, sort objects into an order such that the
0982 /// associated values would be in increasing order (but don't actually
0983 /// touch the values vector in the process).
0984 template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects, 
0985                           const std::vector<double> & values) {
0986   //assert(objects.size() == values.size());
0987   if (objects.size() != values.size()){
0988     throw Error("fastjet::objects_sorted_by_values(...): the size of the 'objects' vector must match the size of the 'values' vector");
0989   }
0990   
0991   // get a vector of indices
0992   std::vector<int> indices(values.size());
0993   for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
0994   
0995   // sort the indices
0996   sort_indices(indices, values);
0997   
0998   // copy the objects 
0999   std::vector<T> objects_sorted(objects.size());
1000   
1001   // place the objects in the correct order
1002   for (size_t i = 0; i < indices.size(); i++) {
1003     objects_sorted[i] = objects[indices[i]];
1004   }
1005 
1006   return objects_sorted;
1007 }
1008 
1009 /// \if internal_doc
1010 /// @ingroup internal
1011 /// \class IndexedSortHelper
1012 /// a class that helps us carry out indexed sorting.
1013 /// \endif
1014 class IndexedSortHelper {
1015 public:
1016   inline IndexedSortHelper (const std::vector<double> * reference_values) {
1017     _ref_values = reference_values;
1018   };
1019   inline int operator() (const int i1, const int i2) const {
1020     return  (*_ref_values)[i1] < (*_ref_values)[i2];
1021   };
1022 private:
1023   const std::vector<double> * _ref_values;
1024 };
1025 
1026 
1027 //----------------------------------------------------------------------
1028 /// constructor from any object that has px,py,pz,E = some_four_vector[0--3],
1029 // NB: do not know if it really needs to be inline, but when it wasn't
1030 //     linking failed with g++ (who knows what was wrong...)
1031 #ifndef SWIG
1032 template <class L> inline  PseudoJet::PseudoJet(const L & some_four_vector) {
1033   reset(some_four_vector);
1034 }
1035 #endif
1036 
1037 //----------------------------------------------------------------------
1038 inline void PseudoJet::_reset_indices() { 
1039   set_cluster_hist_index(-1);
1040   set_user_index(-1);
1041 
1042   _structure.reset();
1043   _user_info.reset();
1044 }
1045 
1046 //std::shared_ptr: inline void PseudoJet::_reset_shared_pointers() {
1047 //std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
1048 //std::shared_ptr:   // if the jet currently belongs to a cs, we need to release it before any chenge
1049 //std::shared_ptr:   _release_jet_from_cs();
1050 //std::shared_ptr: #endif // FASTJET_HAVE_THREAD_SAFETY
1051 //std::shared_ptr:   _structure.reset();
1052 //std::shared_ptr:   _user_info.reset();
1053 //std::shared_ptr:   // and do not forget to remove it in reset_indices above
1054 //std::shared_ptr: }
1055 
1056 
1057 
1058 
1059 // taken literally from CLHEP
1060 inline double PseudoJet::m() const {
1061   double mm = m2();
1062   return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
1063 }
1064 
1065 
1066 inline void PseudoJet::reset(double px_in, double py_in, double pz_in, double E_in) {
1067   _px = px_in;
1068   _py = py_in;
1069   _pz = pz_in;
1070   _E  = E_in;
1071   _finish_init();
1072   _reset_indices();
1073   //std::shared_ptr: _reset_shared_pointers();
1074 }
1075 
1076 inline void PseudoJet::reset_momentum(double px_in, double py_in, double pz_in, double E_in) {
1077   _px = px_in;
1078   _py = py_in;
1079   _pz = pz_in;
1080   _E  = E_in;
1081   _finish_init();
1082 }
1083 
1084 inline void PseudoJet::reset_momentum(const PseudoJet & pj) {
1085   _px  = pj._px ;
1086   _py  = pj._py ;
1087   _pz  = pj._pz ;
1088   _E   = pj._E  ;
1089   _phi = pj._phi;
1090   _rap = pj._rap;
1091   _kt2 = pj._kt2;
1092 }
1093 
1094 //-------------------------------------------------------------------------------
1095 // implementation of the templated accesses to the underlying structyre
1096 //-------------------------------------------------------------------------------
1097 
1098 // returns a reference to the structure casted to the requested
1099 // structure type
1100 //
1101 // If there is no sructure associated, an Error is thrown.
1102 // If the type is not met, a std::bad_cast error is thrown.
1103 template<typename StructureType>
1104 const StructureType & PseudoJet::structure() const{
1105   return dynamic_cast<const StructureType &>(* validated_structure_ptr());
1106   
1107 }
1108 
1109 // check if the PseudoJet has the structure resulting from a Transformer 
1110 // (that is, its structure is compatible with a Transformer::StructureType)
1111 template<typename TransformerType>
1112 bool PseudoJet::has_structure_of() const{
1113   if (!_structure) return false;
1114 
1115   return dynamic_cast<const typename TransformerType::StructureType *>(_structure.get()) != 0;
1116 }
1117 
1118 // this is a helper to access a structure created by a Transformer 
1119 // (that is, of type Transformer::StructureType)
1120 // NULL is returned if the corresponding type is not met
1121 template<typename TransformerType>
1122 const typename TransformerType::StructureType & PseudoJet::structure_of() const{
1123   if (!_structure) 
1124     throw Error("Trying to access the structure of a PseudoJet without an associated structure");
1125 
1126   return dynamic_cast<const typename TransformerType::StructureType &>(*_structure);
1127 }
1128 
1129 
1130 
1131 //-------------------------------------------------------------------------------
1132 // helper functions to build a jet made of pieces
1133 //
1134 // Note that there are more complete versions of these functions, with
1135 // an additional argument for a recombination scheme, in
1136 // JetDefinition.hh
1137 // -------------------------------------------------------------------------------
1138 
1139 /// build a "CompositeJet" from the vector of its pieces
1140 ///
1141 /// In this case, E-scheme recombination is assumed to compute the
1142 /// total momentum
1143 PseudoJet join(const std::vector<PseudoJet> & pieces);
1144 
1145 /// build a MergedJet from a single PseudoJet
1146 PseudoJet join(const PseudoJet & j1);
1147 
1148 /// build a MergedJet from 2 PseudoJet
1149 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2);
1150 
1151 /// build a MergedJet from 3 PseudoJet
1152 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3);
1153 
1154 /// build a MergedJet from 4 PseudoJet
1155 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4);
1156 
1157 
1158 
1159 FASTJET_END_NAMESPACE
1160 
1161 #endif // __FASTJET_PSEUDOJET_HH__