FastJet 3.0.0
|
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__