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