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