FastJet 3.0.4
PseudoJet.hh
00001 //STARTHEADER
00002 // $Id: PseudoJet.hh 3111 2013-05-04 08:17:27Z salam $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 
00030 #ifndef __FASTJET_PSEUDOJET_HH__
00031 #define __FASTJET_PSEUDOJET_HH__
00032 
00033 #include<valarray>
00034 #include<vector>
00035 #include<cassert>
00036 #include<cmath>
00037 #include<iostream>
00038 #include "fastjet/internal/numconsts.hh"
00039 #include "fastjet/internal/IsBase.hh"
00040 #include "fastjet/SharedPtr.hh"
00041 #include "fastjet/Error.hh"
00042 #include "fastjet/PseudoJetStructureBase.hh"
00043 
00044 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00045 
00046 //using namespace std;
00047 
00048 /// Used to protect against parton-level events where pt can be zero
00049 /// for some partons, giving rapidity=infinity. KtJet fails in those cases.
00050 const double MaxRap = 1e5;
00051 
00052 /// default value for phi, meaning it (and rapidity) have yet to be calculated) 
00053 const double pseudojet_invalid_phi = -100.0;
00054 const double pseudojet_invalid_rap = -1e200;
00055 
00056 #ifndef __FJCORE__
00057 // forward definition
00058 class ClusterSequenceAreaBase;
00059 #endif  // __FJCORE__
00060 
00061 /// @ingroup basic_classes
00062 /// \class PseudoJet
00063 /// Class to contain pseudojets, including minimal information of use to
00064 /// jet-clustering routines.
00065 class PseudoJet {
00066 
00067  public:
00068   //----------------------------------------------------------------------
00069   /// @name Constructors and destructor
00070   //\{
00071   /// default constructor, which as of FJ3.0 provides an object for
00072   /// which all operations are now valid and which has zero momentum
00073   ///
00074   // (cf. this is actually OK from a timing point of view and in some
00075   // cases better than just having the default constructor for the
00076   // internal shared pointer: see PJtiming.cc and the notes therein)
00077   PseudoJet() : _px(0), _py(0), _pz(0), _E(0) {_finish_init(); _reset_indices();}
00078   //PseudoJet() : _px(0), _py(0), _pz(0), _E(0), _phi(pseudojet_invalid_phi), _rap(pseudojet_invalid_rap), _kt2(0) {_reset_indices();}
00079   /// construct a pseudojet from explicit components
00080   PseudoJet(const double px, const double py, const double pz, const double E);
00081 
00082   /// constructor from any object that has px,py,pz,E = some_four_vector[0--3],
00083   template <class L> PseudoJet(const L & some_four_vector);
00084 
00085   // Constructor that performs minimal initialisation (only that of
00086   // the shared pointers), of use in certain speed-critical contexts
00087   //
00088   // NB: "dummy" is commented to avoid unused-variable compiler warnings
00089   PseudoJet(bool /* dummy */) {}
00090 
00091   /// default (virtual) destructor
00092   virtual ~PseudoJet(){};
00093   //\} ---- end of constructors and destructors --------------------------
00094 
00095   //----------------------------------------------------------------------
00096   /// @name Kinematic access functions
00097   //\{
00098   //----------------------------------------------------------------------
00099   inline double E()   const {return _E;}
00100   inline double e()   const {return _E;} // like CLHEP
00101   inline double px()  const {return _px;}
00102   inline double py()  const {return _py;}
00103   inline double pz()  const {return _pz;}
00104 
00105   /// returns phi (in the range 0..2pi)
00106   inline double phi() const {return phi_02pi();}
00107 
00108   /// returns phi in the range -pi..pi
00109   inline double phi_std()  const {
00110     _ensure_valid_rap_phi();
00111     return _phi > pi ? _phi-twopi : _phi;}
00112 
00113   /// returns phi in the range 0..2pi
00114   inline double phi_02pi() const {
00115     _ensure_valid_rap_phi();
00116     return _phi;
00117   }
00118 
00119   /// returns the rapidity or some large value when the rapidity
00120   /// is infinite
00121   inline double rap() const {
00122     _ensure_valid_rap_phi();
00123     return _rap;
00124   }
00125 
00126   /// the same as rap()
00127   inline double rapidity() const {return rap();} // like CLHEP
00128 
00129   /// returns the pseudo-rapidity or some large value when the
00130   /// rapidity is infinite
00131   double pseudorapidity() const;
00132   double eta() const {return pseudorapidity();}
00133 
00134   /// returns the squared transverse momentum
00135   inline double pt2() const {return _kt2;}
00136   /// returns the scalar transverse momentum
00137   inline double  pt() const {return sqrt(_kt2);} 
00138   /// returns the squared transverse momentum
00139   inline double perp2() const {return _kt2;}  // like CLHEP
00140   /// returns the scalar transverse momentum
00141   inline double  perp() const {return sqrt(_kt2);}    // like CLHEP
00142   /// returns the squared transverse momentum
00143   inline double kt2() const {return _kt2;} // for bkwds compatibility
00144 
00145   /// returns the squared invariant mass // like CLHEP
00146   inline double  m2() const {return (_E+_pz)*(_E-_pz)-_kt2;}    
00147   /// returns the invariant mass 
00148   /// (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP)
00149   inline double  m() const;    
00150 
00151   /// returns the squared transverse mass = kt^2+m^2
00152   inline double mperp2() const {return (_E+_pz)*(_E-_pz);}
00153   /// returns the transverse mass = sqrt(kt^2+m^2)
00154   inline double mperp() const {return sqrt(std::abs(mperp2()));}
00155   /// returns the squared transverse mass = kt^2+m^2
00156   inline double mt2() const {return (_E+_pz)*(_E-_pz);}
00157   /// returns the transverse mass = sqrt(kt^2+m^2)
00158   inline double mt() const {return sqrt(std::abs(mperp2()));}
00159 
00160   /// return the squared 3-vector modulus = px^2+py^2+pz^2
00161   inline double modp2() const {return _kt2+_pz*_pz;}
00162   /// return the 3-vector modulus = sqrt(px^2+py^2+pz^2)
00163   inline double modp() const {return sqrt(_kt2+_pz*_pz);}
00164 
00165   /// return the transverse energy
00166   inline double Et() const {return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);}
00167   /// return the transverse energy squared
00168   inline double Et2() const {return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);}
00169 
00170   /// returns component i, where X==0, Y==1, Z==2, E==3
00171   double operator () (int i) const ; 
00172   /// returns component i, where X==0, Y==1, Z==2, E==3
00173   inline double operator [] (int i) const { return (*this)(i); }; // this too
00174 
00175 
00176 
00177   /// returns kt distance (R=1) between this jet and another
00178   double kt_distance(const PseudoJet & other) const;
00179 
00180   /// returns squared cylinder (rap-phi) distance between this jet and another
00181   double plain_distance(const PseudoJet & other) const;
00182   /// returns squared cylinder (rap-phi) distance between this jet and
00183   /// another
00184   inline double squared_distance(const PseudoJet & other) const {
00185     return plain_distance(other);}
00186 
00187   /// return the cylinder (rap-phi) distance between this jet and another,
00188   /// \f$\Delta_R = \sqrt{\Delta y^2 + \Delta \phi^2}\f$.
00189   inline double delta_R(const PseudoJet & other) const {
00190     return sqrt(squared_distance(other));
00191   }
00192 
00193   /// returns other.phi() - this.phi(), constrained to be in 
00194   /// range -pi .. pi
00195   double delta_phi_to(const PseudoJet & other) const;
00196 
00197   //// this seemed to compile except if it was used
00198   //friend inline double 
00199   //  kt_distance(const PseudoJet & jet1, const PseudoJet & jet2) { 
00200   //                                      return jet1.kt_distance(jet2);}
00201 
00202   /// returns distance between this jet and the beam
00203   inline double beam_distance() const {return _kt2;}
00204 
00205   /// return a valarray containing the four-momentum (components 0-2
00206   /// are 3-mom, component 3 is energy).
00207   std::valarray<double> four_mom() const;
00208 
00209   //\}  ------- end of kinematic access functions
00210 
00211   // taken from CLHEP
00212   enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
00213 
00214 
00215   //----------------------------------------------------------------------
00216   /// @name Kinematic modification functions
00217   //\{
00218   //----------------------------------------------------------------------
00219   /// transform this jet (given in the rest frame of prest) into a jet
00220   /// in the lab frame [NOT FULLY TESTED]
00221   PseudoJet & boost(const PseudoJet & prest);
00222   /// transform this jet (given in lab) into a jet in the rest
00223   /// frame of prest  [NOT FULLY TESTED]
00224   PseudoJet & unboost(const PseudoJet & prest);
00225 
00226   void operator*=(double);
00227   void operator/=(double);
00228   void operator+=(const PseudoJet &);
00229   void operator-=(const PseudoJet &);
00230 
00231   /// reset the 4-momentum according to the supplied components and
00232   /// put the user and history indices back to their default values
00233   inline void reset(double px, double py, double pz, double E);
00234 
00235   /// reset the PseudoJet to be equal to psjet (including its
00236   /// indices); NB if the argument is derived from a PseudoJet then
00237   /// the "reset" used will be the templated version
00238   ///
00239   /// Note: this is included on top of the templated version because
00240   /// PseudoJet is not "derived" from PseudoJet, so the templated
00241   /// reset would not handle this case properly.
00242   inline void reset(const PseudoJet & psjet) {
00243     (*this) = psjet;
00244   }
00245 
00246   /// reset the 4-momentum according to the supplied generic 4-vector
00247   /// (accessible via indexing, [0]==px,...[3]==E) and put the user
00248   /// and history indices back to their default values.
00249   template <class L> inline void reset(const L & some_four_vector) {
00250     // check if some_four_vector can be cast to a PseudoJet
00251     //
00252     // Note that a regular dynamic_cast would not work here because
00253     // there is no guarantee that L is polymorphic. We use a more
00254     // complex construct here that works also in such a case. As for
00255     // dynamic_cast, NULL is returned if L is not derived from
00256     // PseudoJet
00257     //
00258     // Note the explicit request for fastjet::cast_if_derived; when
00259     // combining fastjet and fjcore, this avoids ambiguity in which of
00260     // the two cast_if_derived calls to use.
00261     const PseudoJet * pj = fastjet::cast_if_derived<const PseudoJet>(&some_four_vector);
00262 
00263     if (pj){
00264       (*this) = *pj;
00265     } else {
00266       reset(some_four_vector[0], some_four_vector[1],
00267             some_four_vector[2], some_four_vector[3]);
00268     }
00269   }
00270 
00271   /// reset the PseudoJet according to the specified pt, rapidity,
00272   /// azimuth and mass (also resetting indices, etc.)
00273   /// (phi should satisfy -2pi<phi<4pi)
00274   inline void reset_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in=0.0) {
00275     reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in);
00276     _reset_indices();
00277   }
00278 
00279   /// reset the 4-momentum according to the supplied components 
00280   /// but leave all other information (indices, user info, etc.)
00281   /// untouched
00282   inline void reset_momentum(double px, double py, double pz, double E);
00283 
00284   /// reset the 4-momentum according to the components of the supplied
00285   /// PseudoJet, including cached components; note that the template
00286   /// version (below) will be called for classes derived from PJ.
00287   inline void reset_momentum(const PseudoJet & pj);
00288 
00289   /// reset the 4-momentum according to the specified pt, rapidity,
00290   /// azimuth and mass (phi should satisfy -2pi<phi<4pi)
00291   void reset_momentum_PtYPhiM(double pt, double y, double phi, double m=0.0);
00292 
00293   /// reset the 4-momentum according to the supplied generic 4-vector
00294   /// (accessible via indexing, [0]==px,...[3]==E), but leave all
00295   /// other information (indices, user info, etc.)  untouched
00296   template <class L> inline void reset_momentum(const L & some_four_vector) {
00297     reset_momentum(some_four_vector[0], some_four_vector[1],
00298                    some_four_vector[2], some_four_vector[3]);
00299   }
00300 
00301   /// in some cases when setting a 4-momentum, the user/program knows
00302   /// what rapidity and azimuth are associated with that 4-momentum;
00303   /// by calling this routine the user can provide the information
00304   /// directly to the PseudoJet and avoid expensive rap-phi
00305   /// recalculations.
00306   ///
00307   /// - \param rap  rapidity
00308   /// - \param phi  (in range -twopi...4*pi)
00309   ///
00310   /// USE WITH CAUTION: there are no checks that the rapidity and
00311   /// azimuth supplied are sensible, nor does this reset the
00312   /// 4-momentum components if things don't match.
00313   void set_cached_rap_phi(double rap, double phi);
00314 
00315 
00316   //\} --- end of kin mod functions ------------------------------------
00317 
00318   //----------------------------------------------------------------------
00319   /// @name User index functions
00320   ///
00321   /// To allow the user to set and access an integer index which can
00322   /// be exploited by the user to associate extra information with a
00323   /// particle/jet (for example pdg id, or an indication of a
00324   /// particle's origin within the user's analysis)
00325   //
00326   //\{
00327 
00328   /// return the user_index, 
00329   inline int user_index() const {return _user_index;}
00330   /// set the user_index, intended to allow the user to add simple
00331   /// identifying information to a particle/jet
00332   inline void set_user_index(const int index) {_user_index = index;}
00333 
00334   //\} ----- end of use index functions ---------------------------------
00335 
00336   //----------------------------------------------------------------------
00337   /// @name User information types and functions
00338   ///
00339   /// Allows PseudoJet to carry extra user info (as an object derived from
00340   /// UserInfoBase).
00341   //\{
00342 
00343   /// @ingroup user_info
00344   /// \class UserInfoBase
00345   /// a base class to hold extra user information in a PseudoJet
00346   ///
00347   /// This is a base class to help associate extra user information
00348   /// with a jet. The user should store their information in a class
00349   /// derived from this. This allows information of arbitrary
00350   /// complexity to be easily associated with a PseudoJet (in contrast
00351   /// to the user index). For example, in a Monte Carlo simulation,
00352   /// the user information might include the PDG ID, and the position
00353   /// of the production vertex for the particle.
00354   ///
00355   /// The PseudoJet is able to store a shared pointer to any object
00356   /// derived from UserInfo. The use of a shared pointer frees the
00357   /// user of the need to handle the memory management associated with
00358   /// the information.
00359   ///
00360   /// Having the user information derive from a common base class also
00361   /// facilitates dynamic casting, etc.
00362   ///
00363   class UserInfoBase{
00364   public:
00365     // dummy ctor
00366     UserInfoBase(){};
00367 
00368     // dummy virtual dtor
00369     // makes it polymorphic to allow for dynamic_cast
00370     virtual ~UserInfoBase(){}; 
00371   };
00372 
00373   /// error class to be thrown if accessing user info when it doesn't
00374   /// exist
00375   class InexistentUserInfo : public Error {
00376   public:
00377     InexistentUserInfo();
00378   };
00379 
00380   /// sets the internal shared pointer to the user information.
00381   ///
00382   /// Note that the PseudoJet will now _own_ the pointer, and delete
00383   /// the corresponding object when it (the jet, and any copies of the jet)
00384   /// goes out of scope. 
00385   void set_user_info(UserInfoBase * user_info_in) {
00386     _user_info.reset(user_info_in);
00387   }
00388 
00389   /// returns a reference to the dynamic cast conversion of user_info
00390   /// to type L.
00391   ///
00392   /// Usage: suppose you have previously set the user info with a pointer
00393   /// to an object of type MyInfo, 
00394   ///
00395   ///   class MyInfo: public PseudoJet::UserInfoBase {
00396   ///      MyInfo(int id) : _pdg_id(id);
00397   ///      int pdg_id() const {return _pdg_id;}
00398   ///      int _pdg_id;
00399   ///   };
00400   ///
00401   ///   PseudoJet particle(...);
00402   ///   particle.set_user_info(new MyInfo(its_pdg_id));
00403   ///
00404   /// Then you would access that pdg_id() as
00405   ///
00406   ///   particle.user_info<MyInfo>().pdg_id();
00407   ///
00408   /// It's overkill for just a single integer, but scales easily to
00409   /// more extensive information.
00410   ///
00411   /// Note that user_info() throws an InexistentUserInfo() error if
00412   /// there is no user info; throws a std::bad_cast if the conversion
00413   /// doesn't work
00414   ///
00415   /// If this behaviour does not fit your needs, use instead the the
00416   /// user_info_ptr() or user_info_shared_ptr() member functions.
00417   template<class L>
00418   const L & user_info() const{
00419     if (_user_info.get() == 0) throw InexistentUserInfo();
00420     return dynamic_cast<const L &>(* _user_info.get());
00421   }
00422 
00423   /// returns true if the PseudoJet has user information
00424   bool has_user_info() const{
00425     return _user_info.get();
00426   }
00427 
00428   /// returns true if the PseudoJet has user information than can be
00429   /// cast to the template argument type.
00430   template<class L>
00431   bool has_user_info() const{
00432     return _user_info.get() && dynamic_cast<const L *>(_user_info.get());
00433   }
00434 
00435   /// retrieve a pointer to the (const) user information
00436   const UserInfoBase * user_info_ptr() const{
00437     if (!_user_info()) return NULL;
00438     return _user_info.get();
00439   }
00440 
00441 
00442   /// retrieve a (const) shared pointer to the user information
00443   const SharedPtr<UserInfoBase> & user_info_shared_ptr() const{
00444     return _user_info;
00445   }
00446 
00447   /// retrieve a (non-const) shared pointer to the user information;
00448   /// you can use this, for example, to set the shared pointer, eg
00449   ///
00450   /// \code
00451   ///   p2.user_info_shared_ptr() = p1.user_info_shared_ptr();
00452   /// \endcode
00453   ///
00454   /// or 
00455   ///
00456   /// \code
00457   ///   SharedPtr<PseudoJet::UserInfoBase> info_shared(new MyInfo(...));
00458   ///   p2.user_info_shared_ptr() = info_shared;
00459   /// \endcode
00460   SharedPtr<UserInfoBase> & user_info_shared_ptr(){
00461     return _user_info;
00462   }
00463 
00464   // \} --- end of extra info functions ---------------------------------
00465 
00466   //----------------------------------------------------------------------
00467   /// @name Description
00468   ///
00469   /// Since a PseudoJet can have a structure that contains a variety
00470   /// of information, we provide a description that allows one to check
00471   /// exactly what kind of PseudoJet we are dealing with
00472   //
00473   //\{
00474 
00475   /// return a string describing what kind of PseudoJet we are dealing with 
00476   std::string description() const;
00477 
00478   //\} ----- end of description functions ---------------------------------
00479 
00480   //-------------------------------------------------------------
00481   /// @name Access to the associated ClusterSequence object.
00482   ///
00483   /// In addition to having kinematic information, jets may contain a
00484   /// reference to an associated ClusterSequence (this is the case,
00485   /// for example, if the jet has been returned by a ClusterSequence
00486   /// member function).
00487   //\{
00488   //-------------------------------------------------------------
00489   /// returns true if this PseudoJet has an associated ClusterSequence.
00490   bool has_associated_cluster_sequence() const;
00491   /// shorthand for has_associated_cluster_sequence()
00492   bool has_associated_cs() const {return has_associated_cluster_sequence();}
00493 
00494   /// returns true if this PseudoJet has an associated and still
00495   /// valid(ated) ClusterSequence.
00496   bool has_valid_cluster_sequence() const;
00497   /// shorthand for has_valid_cluster_sequence()
00498   bool has_valid_cs() const {return has_valid_cluster_sequence();}
00499 
00500   /// get a (const) pointer to the parent ClusterSequence (NULL if
00501   /// inexistent)
00502   const ClusterSequence* associated_cluster_sequence() const;
00503   // shorthand for associated_cluster_sequence()
00504   const ClusterSequence* associated_cs() const {return associated_cluster_sequence();}
00505 
00506   /// if the jet has a valid associated cluster sequence then return a
00507   /// pointer to it; otherwise throw an error
00508   inline const ClusterSequence * validated_cluster_sequence() const {
00509     return validated_cs();
00510   }
00511   /// shorthand for validated_cluster_sequence()
00512   const ClusterSequence * validated_cs() const;
00513 
00514 #ifndef __FJCORE__
00515   /// if the jet has valid area information then return a pointer to
00516   /// the associated ClusterSequenceAreaBase object; otherwise throw an error
00517   inline const ClusterSequenceAreaBase * validated_cluster_sequence_area_base() const {
00518     return validated_csab();
00519   }
00520 
00521   /// shorthand for validated_cluster_sequence_area_base()
00522   const ClusterSequenceAreaBase * validated_csab() const;
00523 #endif  //  __FJCORE__
00524 
00525   //\}
00526 
00527   //-------------------------------------------------------------
00528   /// @name Access to the associated PseudoJetStructureBase object.
00529   ///
00530   /// In addition to having kinematic information, jets may contain a
00531   /// reference to an associated ClusterSequence (this is the case,
00532   /// for example, if the jet has been returned by a ClusterSequence
00533   /// member function).
00534   //\{
00535   //-------------------------------------------------------------
00536 
00537   /// set the associated structure
00538   void set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure);
00539 
00540   /// return true if there is some structure associated with this PseudoJet
00541   bool has_structure() const;
00542 
00543   /// return a pointer to the structure (of type
00544   /// PseudoJetStructureBase*) associated with this PseudoJet.
00545   ///
00546   /// return NULL if there is no associated structure
00547   const PseudoJetStructureBase* structure_ptr() const;
00548   
00549   /// return a non-const pointer to the structure (of type
00550   /// PseudoJetStructureBase*) associated with this PseudoJet.
00551   ///
00552   /// return NULL if there is no associated structure
00553   ///
00554   /// Only use this if you know what you are doing. In any case,
00555   /// prefer the 'structure_ptr()' (the const version) to this method,
00556   /// unless you really need a write access to the PseudoJet's
00557   /// underlying structure.
00558   PseudoJetStructureBase* structure_non_const_ptr();
00559   
00560   /// return a pointer to the structure (of type
00561   /// PseudoJetStructureBase*) associated with this PseudoJet.
00562   ///
00563   /// throw an error if there is no associated structure
00564   const PseudoJetStructureBase* validated_structure_ptr() const;
00565   
00566   /// return a reference to the shared pointer to the
00567   /// PseudoJetStructureBase associated with this PseudoJet
00568   const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const;
00569 
00570   /// returns a reference to the structure casted to the requested
00571   /// structure type
00572   ///
00573   /// If there is no sructure associated, an Error is thrown.
00574   /// If the type is not met, a std::bad_cast error is thrown.
00575   template<typename StructureType>
00576   const StructureType & structure() const;
00577 
00578   /// check if the PseudoJet has the structure resulting from a Transformer 
00579   /// (that is, its structure is compatible with a Transformer::StructureType).
00580   /// If there is no structure, false is returned.
00581   template<typename TransformerType>
00582   bool has_structure_of() const;
00583 
00584   /// this is a helper to access any structure created by a Transformer 
00585   /// (that is, of type Transformer::StructureType).
00586   ///
00587   /// If there is no structure, or if the structure is not compatible
00588   /// with TransformerType, an error is thrown.
00589   template<typename TransformerType>
00590   const typename TransformerType::StructureType & structure_of() const;
00591 
00592   //\}
00593 
00594   //-------------------------------------------------------------
00595   /// @name Methods for access to information about jet structure
00596   ///
00597   /// These allow access to jet constituents, and other jet
00598   /// subtructure information. They only work if the jet is associated
00599   /// with a ClusterSequence.
00600   //-------------------------------------------------------------
00601   //\{
00602 
00603   /// check if it has been recombined with another PseudoJet in which
00604   /// case, return its partner through the argument. Otherwise,
00605   /// 'partner' is set to 0.
00606   ///
00607   /// an Error is thrown if this PseudoJet has no currently valid
00608   /// associated ClusterSequence
00609   virtual bool has_partner(PseudoJet &partner) const;
00610 
00611   /// check if it has been recombined with another PseudoJet in which
00612   /// case, return its child through the argument. Otherwise, 'child'
00613   /// is set to 0.
00614   /// 
00615   /// an Error is thrown if this PseudoJet has no currently valid
00616   /// associated ClusterSequence
00617   virtual bool has_child(PseudoJet &child) const;
00618 
00619   /// check if it is the product of a recombination, in which case
00620   /// return the 2 parents through the 'parent1' and 'parent2'
00621   /// arguments. Otherwise, set these to 0.
00622   ///
00623   /// an Error is thrown if this PseudoJet has no currently valid
00624   /// associated ClusterSequence
00625   virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const;
00626 
00627   /// check if the current PseudoJet contains the one passed as
00628   /// argument.
00629   ///
00630   /// an Error is thrown if this PseudoJet has no currently valid
00631   /// associated ClusterSequence
00632   virtual bool contains(const PseudoJet &constituent) const;
00633 
00634   /// check if the current PseudoJet is contained the one passed as
00635   /// argument.
00636   ///
00637   /// an Error is thrown if this PseudoJet has no currently valid
00638   /// associated ClusterSequence
00639   virtual bool is_inside(const PseudoJet &jet) const;
00640 
00641 
00642   /// returns true if the PseudoJet has constituents
00643   virtual bool has_constituents() const;
00644 
00645   /// retrieve the constituents. 
00646   ///
00647   /// an Error is thrown if this PseudoJet has no currently valid
00648   /// associated ClusterSequence or other substructure information
00649   virtual std::vector<PseudoJet> constituents() const;
00650 
00651 
00652   /// returns true if the PseudoJet has support for exclusive subjets
00653   virtual bool has_exclusive_subjets() const;
00654 
00655   /// return a vector of all subjets of the current jet (in the sense
00656   /// of the exclusive algorithm) that would be obtained when running
00657   /// the algorithm with the given dcut. 
00658   ///
00659   /// Time taken is O(m ln m), where m is the number of subjets that
00660   /// are found. If m gets to be of order of the total number of
00661   /// constituents in the jet, this could be substantially slower than
00662   /// just getting that list of constituents.
00663   ///
00664   /// an Error is thrown if this PseudoJet has no currently valid
00665   /// associated ClusterSequence
00666   std::vector<PseudoJet> exclusive_subjets (const double & dcut) const;
00667 
00668   /// return the size of exclusive_subjets(...); still n ln n with same
00669   /// coefficient, but marginally more efficient than manually taking
00670   /// exclusive_subjets.size()
00671   ///
00672   /// an Error is thrown if this PseudoJet has no currently valid
00673   /// associated ClusterSequence
00674   int n_exclusive_subjets(const double & dcut) const;
00675 
00676   /// return the list of subjets obtained by unclustering the supplied
00677   /// jet down to nsub subjets. Throws an error if there are fewer than
00678   /// nsub particles in the jet.
00679   ///
00680   /// For ClusterSequence type jets, requires nsub ln nsub time
00681   ///
00682   /// An Error is thrown if this PseudoJet has no currently valid
00683   /// associated ClusterSequence
00684   std::vector<PseudoJet> exclusive_subjets (int nsub) const;
00685 
00686   /// return the list of subjets obtained by unclustering the supplied
00687   /// jet down to nsub subjets (or all constituents if there are fewer
00688   /// than nsub).
00689   ///
00690   /// For ClusterSequence type jets, requires nsub ln nsub time
00691   ///
00692   /// An Error is thrown if this PseudoJet has no currently valid
00693   /// associated ClusterSequence
00694   std::vector<PseudoJet> exclusive_subjets_up_to (int nsub) const;
00695 
00696   /// return the dij that was present in the merging nsub+1 -> nsub 
00697   /// subjets inside this jet.
00698   ///
00699   /// an Error is thrown if this PseudoJet has no currently valid
00700   /// associated ClusterSequence
00701   double exclusive_subdmerge(int nsub) const;
00702 
00703   /// return the maximum dij that occurred in the whole event at the
00704   /// stage that the nsub+1 -> nsub merge of subjets occurred inside 
00705   /// this jet.
00706   ///
00707   /// an Error is thrown if this PseudoJet has no currently valid
00708   /// associated ClusterSequence
00709   double exclusive_subdmerge_max(int nsub) const;
00710 
00711 
00712   /// returns true if a jet has pieces
00713   ///
00714   /// By default a single particle or a jet coming from a
00715   /// ClusterSequence have no pieces and this methos will return false.
00716   ///
00717   /// In practice, this is equivalent to have an structure of type
00718   /// CompositeJetStructure.
00719   virtual bool has_pieces() const;
00720 
00721 
00722   /// retrieve the pieces that make up the jet. 
00723   ///
00724   /// If the jet does not support pieces, an error is throw
00725   virtual std::vector<PseudoJet> pieces() const;
00726 
00727 
00728   // the following ones require a computation of the area in the
00729   // parent ClusterSequence (See ClusterSequenceAreaBase for details)
00730   //------------------------------------------------------------------
00731 #ifndef __FJCORE__
00732 
00733   /// check if it has a defined area
00734   virtual bool has_area() const;
00735 
00736   /// return the jet (scalar) area.
00737   /// throws an Error if there is no support for area in the parent CS
00738   virtual double area() const;
00739 
00740   /// return the error (uncertainty) associated with the determination
00741   /// of the area of this jet.
00742   /// throws an Error if there is no support for area in the parent CS
00743   virtual double area_error() const;
00744 
00745   /// return the jet 4-vector area.
00746   /// throws an Error if there is no support for area in the parent CS
00747   virtual PseudoJet area_4vector() const;
00748 
00749   /// true if this jet is made exclusively of ghosts.
00750   /// throws an Error if there is no support for area in the parent CS
00751   virtual bool is_pure_ghost() const;
00752 
00753 #endif  // __FJCORE__
00754   //\} --- end of jet structure -------------------------------------
00755 
00756 
00757 
00758   //----------------------------------------------------------------------
00759   /// @name Members mainly intended for internal use
00760   //----------------------------------------------------------------------
00761   //\{
00762   /// return the cluster_hist_index, intended to be used by clustering
00763   /// routines.
00764   inline int cluster_hist_index() const {return _cluster_hist_index;}
00765   /// set the cluster_hist_index, intended to be used by clustering routines.
00766   inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;}
00767 
00768   /// alternative name for cluster_hist_index() [perhaps more meaningful]
00769   inline int cluster_sequence_history_index() const {
00770     return cluster_hist_index();}
00771   /// alternative name for set_cluster_hist_index(...) [perhaps more
00772   /// meaningful]
00773   inline void set_cluster_sequence_history_index(const int index) {
00774     set_cluster_hist_index(index);}
00775 
00776   //\} ---- end of internal use functions ---------------------------
00777 
00778  protected:  
00779 
00780   SharedPtr<PseudoJetStructureBase> _structure;
00781   SharedPtr<UserInfoBase> _user_info;
00782 
00783 
00784  private: 
00785   // NB: following order must be kept for things to behave sensibly...
00786   double _px,_py,_pz,_E;
00787   mutable double _phi, _rap;
00788   double _kt2; 
00789   int    _cluster_hist_index, _user_index;
00790 
00791   /// calculate phi, rap, kt2 based on the 4-momentum components
00792   void _finish_init();
00793   /// set the indices to default values
00794   void _reset_indices();
00795 
00796   /// ensure that the internal values for rapidity and phi 
00797   /// correspond to 4-momentum structure
00798   inline void _ensure_valid_rap_phi() const {
00799     if (_phi == pseudojet_invalid_phi) _set_rap_phi();
00800   }
00801 
00802   /// set cached rapidity and phi values
00803   void _set_rap_phi() const;
00804 };
00805 
00806 
00807 //----------------------------------------------------------------------
00808 // routines for basic binary operations
00809 
00810 PseudoJet operator+(const PseudoJet &, const PseudoJet &);
00811 PseudoJet operator-(const PseudoJet &, const PseudoJet &);
00812 PseudoJet operator*(double, const PseudoJet &);
00813 PseudoJet operator*(const PseudoJet &, double);
00814 PseudoJet operator/(const PseudoJet &, double);
00815 
00816 /// returns true if the 4 momentum components of the two PseudoJets
00817 /// are identical and all the internal indices (user, cluster_history)
00818 /// + structure and user-info shared pointers are too
00819 bool operator==(const PseudoJet &, const PseudoJet &);
00820 
00821 /// inequality test which is exact opposite of operator==
00822 inline bool operator!=(const PseudoJet & a, const PseudoJet & b) {return !(a==b);}
00823 
00824 /// Can only be used with val=0 and tests whether all four
00825 /// momentum components are equal to val (=0.0)
00826 bool operator==(const PseudoJet & jet, const double val);
00827 
00828 /// Can only be used with val=0 and tests whether at least one of the
00829 /// four momentum components is different from val (=0.0)
00830 inline bool operator!=(const PseudoJet & a, const double & val) {return !(a==val);}
00831 
00832 inline double dot_product(const PseudoJet & a, const PseudoJet & b) {
00833   return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz();
00834 }
00835 
00836 /// returns true if the momenta of the two input jets are identical
00837 bool have_same_momentum(const PseudoJet &, const PseudoJet &);
00838 
00839 /// return a pseudojet with the given pt, y, phi and mass
00840 /// (phi should satisfy -2pi<phi<4pi)
00841 PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0);
00842 
00843 //----------------------------------------------------------------------
00844 // Routines to do with providing sorted arrays of vectors.
00845 
00846 /// return a vector of jets sorted into decreasing transverse momentum
00847 std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets);
00848 
00849 /// return a vector of jets sorted into increasing rapidity
00850 std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets);
00851 
00852 /// return a vector of jets sorted into decreasing energy
00853 std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets);
00854 
00855 /// return a vector of jets sorted into increasing pz
00856 std::vector<PseudoJet> sorted_by_pz(const std::vector<PseudoJet> & jets);
00857 
00858 //----------------------------------------------------------------------
00859 // some code to help sorting
00860 
00861 /// sort the indices so that values[indices[0->n-1]] is sorted
00862 /// into increasing order 
00863 void sort_indices(std::vector<int> & indices, 
00864                   const std::vector<double> & values);
00865 
00866 /// given a vector of values with a one-to-one correspondence with the
00867 /// vector of objects, sort objects into an order such that the
00868 /// associated values would be in increasing order (but don't actually
00869 /// touch the values vector in the process).
00870 template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects, 
00871                                               const std::vector<double> & values);
00872 
00873 /// \if internal_doc
00874 /// @ingroup internal
00875 /// \class IndexedSortHelper
00876 /// a class that helps us carry out indexed sorting.
00877 /// \endif
00878 class IndexedSortHelper {
00879 public:
00880   inline IndexedSortHelper (const std::vector<double> * reference_values) {
00881     _ref_values = reference_values;
00882   };
00883   inline int operator() (const int & i1, const int & i2) const {
00884     return  (*_ref_values)[i1] < (*_ref_values)[i2];
00885   };
00886 private:
00887   const std::vector<double> * _ref_values;
00888 };
00889 
00890 
00891 //----------------------------------------------------------------------
00892 /// constructor from any object that has px,py,pz,E = some_four_vector[0--3],
00893 // NB: do not know if it really needs to be inline, but when it wasn't
00894 //     linking failed with g++ (who knows what was wrong...)
00895 template <class L> inline  PseudoJet::PseudoJet(const L & some_four_vector) {
00896   reset(some_four_vector);
00897 }
00898 
00899 //----------------------------------------------------------------------
00900 inline void PseudoJet::_reset_indices() { 
00901   set_cluster_hist_index(-1);
00902   set_user_index(-1);
00903   _structure.reset();
00904   _user_info.reset();
00905 }
00906 
00907 
00908 // taken literally from CLHEP
00909 inline double PseudoJet::m() const {
00910   double mm = m2();
00911   return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
00912 }
00913 
00914 
00915 inline void PseudoJet::reset(double px_in, double py_in, double pz_in, double E_in) {
00916   _px = px_in;
00917   _py = py_in;
00918   _pz = pz_in;
00919   _E  = E_in;
00920   _finish_init();
00921   _reset_indices();
00922 }
00923 
00924 inline void PseudoJet::reset_momentum(double px_in, double py_in, double pz_in, double E_in) {
00925   _px = px_in;
00926   _py = py_in;
00927   _pz = pz_in;
00928   _E  = E_in;
00929   _finish_init();
00930 }
00931 
00932 inline void PseudoJet::reset_momentum(const PseudoJet & pj) {
00933   _px  = pj._px ;
00934   _py  = pj._py ;
00935   _pz  = pj._pz ;
00936   _E   = pj._E  ;
00937   _phi = pj._phi;
00938   _rap = pj._rap;
00939   _kt2 = pj._kt2;
00940 }
00941 
00942 //-------------------------------------------------------------------------------
00943 // implementation of the templated accesses to the underlying structyre
00944 //-------------------------------------------------------------------------------
00945 
00946 // returns a reference to the structure casted to the requested
00947 // structure type
00948 //
00949 // If there is no sructure associated, an Error is thrown.
00950 // If the type is not met, a std::bad_cast error is thrown.
00951 template<typename StructureType>
00952 const StructureType & PseudoJet::structure() const{
00953   return dynamic_cast<const StructureType &>(* validated_structure_ptr());
00954   
00955 }
00956 
00957 // check if the PseudoJet has the structure resulting from a Transformer 
00958 // (that is, its structure is compatible with a Transformer::StructureType)
00959 template<typename TransformerType>
00960 bool PseudoJet::has_structure_of() const{
00961   if (!_structure()) return false;
00962 
00963   return dynamic_cast<const typename TransformerType::StructureType *>(_structure.get()) != 0;
00964 }
00965 
00966 // this is a helper to access a structure created by a Transformer 
00967 // (that is, of type Transformer::StructureType)
00968 // NULL is returned if the corresponding type is not met
00969 template<typename TransformerType>
00970 const typename TransformerType::StructureType & PseudoJet::structure_of() const{
00971   if (!_structure()) 
00972     throw Error("Trying to access the structure of a PseudoJet without an associated structure");
00973 
00974   return dynamic_cast<const typename TransformerType::StructureType &>(*_structure);
00975 }
00976 
00977 
00978 
00979 //-------------------------------------------------------------------------------
00980 // helper functions to build a jet made of pieces
00981 //
00982 // Note that there are more complete versions of these functions, with
00983 // an additional argument for a recombination scheme, in
00984 // JetDefinition.hh
00985 // -------------------------------------------------------------------------------
00986 
00987 /// build a "CompositeJet" from the vector of its pieces
00988 ///
00989 /// In this case, E-scheme recombination is assumed to compute the
00990 /// total momentum
00991 PseudoJet join(const std::vector<PseudoJet> & pieces);
00992 
00993 /// build a MergedJet from a single PseudoJet
00994 PseudoJet join(const PseudoJet & j1);
00995 
00996 /// build a MergedJet from 2 PseudoJet
00997 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2);
00998 
00999 /// build a MergedJet from 3 PseudoJet
01000 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3);
01001 
01002 /// build a MergedJet from 4 PseudoJet
01003 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4);
01004 
01005 
01006 
01007 FASTJET_END_NAMESPACE
01008 
01009 #endif // __FASTJET_PSEUDOJET_HH__
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends