[1] | 1 | //STARTHEADER |
---|
| 2 | // $Id: PseudoJet.hh 2728 2011-11-20 14:18:59Z salam $ |
---|
| 3 | // |
---|
| 4 | // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez |
---|
| 5 | // |
---|
| 6 | //---------------------------------------------------------------------- |
---|
| 7 | // This file is part of FastJet. |
---|
| 8 | // |
---|
| 9 | // FastJet is free software; you can redistribute it and/or modify |
---|
| 10 | // it under the terms of the GNU General Public License as published by |
---|
| 11 | // the Free Software Foundation; either version 2 of the License, or |
---|
| 12 | // (at your option) any later version. |
---|
| 13 | // |
---|
| 14 | // The algorithms that underlie FastJet have required considerable |
---|
| 15 | // development and are described in hep-ph/0512210. If you use |
---|
| 16 | // FastJet as part of work towards a scientific publication, please |
---|
| 17 | // include a citation to the FastJet paper. |
---|
| 18 | // |
---|
| 19 | // FastJet is distributed in the hope that it will be useful, |
---|
| 20 | // but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 21 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
| 22 | // GNU General Public License for more details. |
---|
| 23 | // |
---|
| 24 | // You should have received a copy of the GNU General Public License |
---|
| 25 | // along with FastJet. If not, see <http://www.gnu.org/licenses/>. |
---|
| 26 | //---------------------------------------------------------------------- |
---|
| 27 | //ENDHEADER |
---|
| 28 | |
---|
| 29 | |
---|
| 30 | #ifndef __FASTJET_PSEUDOJET_HH__ |
---|
| 31 | #define __FASTJET_PSEUDOJET_HH__ |
---|
| 32 | |
---|
| 33 | #include<valarray> |
---|
| 34 | #include<vector> |
---|
| 35 | #include<cassert> |
---|
| 36 | #include<cmath> |
---|
| 37 | #include<iostream> |
---|
| 38 | #include "fastjet/internal/numconsts.hh" |
---|
| 39 | #include "fastjet/internal/IsBase.hh" |
---|
| 40 | #include "fastjet/SharedPtr.hh" |
---|
| 41 | #include "fastjet/Error.hh" |
---|
| 42 | #include "fastjet/PseudoJetStructureBase.hh" |
---|
| 43 | |
---|
| 44 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
| 45 | |
---|
| 46 | //using namespace std; |
---|
| 47 | |
---|
| 48 | /// Used to protect against parton-level events where pt can be zero |
---|
| 49 | /// for some partons, giving rapidity=infinity. KtJet fails in those cases. |
---|
| 50 | const double MaxRap = 1e5; |
---|
| 51 | |
---|
| 52 | /// default value for phi, meaning it (and rapidity) have yet to be calculated) |
---|
| 53 | const double pseudojet_invalid_phi = -100.0; |
---|
| 54 | |
---|
| 55 | // forward definition |
---|
| 56 | class ClusterSequenceAreaBase; |
---|
| 57 | |
---|
| 58 | /// @ingroup basic_classes |
---|
| 59 | /// \class PseudoJet |
---|
| 60 | /// Class to contain pseudojets, including minimal information of use to |
---|
| 61 | /// jet-clustering routines. |
---|
| 62 | class PseudoJet { |
---|
| 63 | |
---|
| 64 | public: |
---|
| 65 | //---------------------------------------------------------------------- |
---|
| 66 | /// @name Constructors and destructor |
---|
| 67 | //\{ |
---|
| 68 | /// default constructor, which as of FJ3.0 provides an object for |
---|
| 69 | /// which all operations are now valid and which has zero momentum |
---|
| 70 | /// |
---|
| 71 | // (cf. this is actually OK from a timing point of view and in some |
---|
| 72 | // cases better than just having the default constructor for the |
---|
| 73 | // internal shared pointer: see PJtiming.cc and the notes therein) |
---|
| 74 | PseudoJet() : _px(0), _py(0), _pz(0), _E(0) {_finish_init(); _reset_indices();} |
---|
| 75 | /// construct a pseudojet from explicit components |
---|
| 76 | PseudoJet(const double px, const double py, const double pz, const double E); |
---|
| 77 | |
---|
| 78 | /// constructor from any object that has px,py,pz,E = some_four_vector[0--3], |
---|
| 79 | template <class L> PseudoJet(const L & some_four_vector); |
---|
| 80 | |
---|
| 81 | // Constructor that performs minimal initialisation (only that of |
---|
| 82 | // the shared pointers), of use in certain speed-critical contexts |
---|
| 83 | // |
---|
| 84 | // NB: "dummy" is commented to avoid unused-variable compiler warnings |
---|
| 85 | PseudoJet(bool /* dummy */) {} |
---|
| 86 | |
---|
| 87 | /// default (virtual) destructor |
---|
| 88 | virtual ~PseudoJet(){}; |
---|
| 89 | //\} ---- end of constructors and destructors -------------------------- |
---|
| 90 | |
---|
| 91 | //---------------------------------------------------------------------- |
---|
| 92 | /// @name Kinematic access functions |
---|
| 93 | //\{ |
---|
| 94 | //---------------------------------------------------------------------- |
---|
| 95 | inline double E() const {return _E;} |
---|
| 96 | inline double e() const {return _E;} // like CLHEP |
---|
| 97 | inline double px() const {return _px;} |
---|
| 98 | inline double py() const {return _py;} |
---|
| 99 | inline double pz() const {return _pz;} |
---|
| 100 | |
---|
| 101 | /// returns phi (in the range 0..2pi) |
---|
| 102 | inline double phi() const {return phi_02pi();} |
---|
| 103 | |
---|
| 104 | /// returns phi in the range -pi..pi |
---|
| 105 | inline double phi_std() const { |
---|
| 106 | _ensure_valid_rap_phi(); |
---|
| 107 | return _phi > pi ? _phi-twopi : _phi;} |
---|
| 108 | |
---|
| 109 | /// returns phi in the range 0..2pi |
---|
| 110 | inline double phi_02pi() const { |
---|
| 111 | _ensure_valid_rap_phi(); |
---|
| 112 | return _phi; |
---|
| 113 | } |
---|
| 114 | |
---|
| 115 | /// returns the rapidity or some large value when the rapidity |
---|
| 116 | /// is infinite |
---|
| 117 | inline double rap() const { |
---|
| 118 | _ensure_valid_rap_phi(); |
---|
| 119 | return _rap; |
---|
| 120 | } |
---|
| 121 | |
---|
| 122 | /// the same as rap() |
---|
| 123 | inline double rapidity() const {return rap();} // like CLHEP |
---|
| 124 | |
---|
| 125 | /// returns the pseudo-rapidity or some large value when the |
---|
| 126 | /// rapidity is infinite |
---|
| 127 | double pseudorapidity() const; |
---|
| 128 | double eta() const {return pseudorapidity();} |
---|
| 129 | |
---|
| 130 | /// returns the squared transverse momentum |
---|
| 131 | inline double pt2() const {return _kt2;} |
---|
| 132 | /// returns the scalar transverse momentum |
---|
| 133 | inline double pt() const {return sqrt(_kt2);} |
---|
| 134 | /// returns the squared transverse momentum |
---|
| 135 | inline double perp2() const {return _kt2;} // like CLHEP |
---|
| 136 | /// returns the scalar transverse momentum |
---|
| 137 | inline double perp() const {return sqrt(_kt2);} // like CLHEP |
---|
| 138 | /// returns the squared transverse momentum |
---|
| 139 | inline double kt2() const {return _kt2;} // for bkwds compatibility |
---|
| 140 | |
---|
| 141 | /// returns the squared invariant mass // like CLHEP |
---|
| 142 | inline double m2() const {return (_E+_pz)*(_E-_pz)-_kt2;} |
---|
| 143 | /// returns the invariant mass |
---|
| 144 | /// (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) |
---|
| 145 | inline double m() const; |
---|
| 146 | |
---|
| 147 | /// returns the squared transverse mass = kt^2+m^2 |
---|
| 148 | inline double mperp2() const {return (_E+_pz)*(_E-_pz);} |
---|
| 149 | /// returns the transverse mass = sqrt(kt^2+m^2) |
---|
| 150 | inline double mperp() const {return sqrt(std::abs(mperp2()));} |
---|
| 151 | /// returns the squared transverse mass = kt^2+m^2 |
---|
| 152 | inline double mt2() const {return (_E+_pz)*(_E-_pz);} |
---|
| 153 | /// returns the transverse mass = sqrt(kt^2+m^2) |
---|
| 154 | inline double mt() const {return sqrt(std::abs(mperp2()));} |
---|
| 155 | |
---|
| 156 | /// return the squared 3-vector modulus = px^2+py^2+pz^2 |
---|
| 157 | inline double modp2() const {return _kt2+_pz*_pz;} |
---|
| 158 | /// return the 3-vector modulus = sqrt(px^2+py^2+pz^2) |
---|
| 159 | inline double modp() const {return sqrt(_kt2+_pz*_pz);} |
---|
| 160 | |
---|
| 161 | /// return the transverse energy |
---|
| 162 | inline double Et() const {return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);} |
---|
| 163 | /// return the transverse energy squared |
---|
| 164 | inline double Et2() const {return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);} |
---|
| 165 | |
---|
| 166 | /// returns component i, where X==0, Y==1, Z==2, E==3 |
---|
| 167 | double operator () (int i) const ; |
---|
| 168 | /// returns component i, where X==0, Y==1, Z==2, E==3 |
---|
| 169 | inline double operator [] (int i) const { return (*this)(i); }; // this too |
---|
| 170 | |
---|
| 171 | |
---|
| 172 | |
---|
| 173 | /// returns kt distance (R=1) between this jet and another |
---|
| 174 | double kt_distance(const PseudoJet & other) const; |
---|
| 175 | |
---|
| 176 | /// returns squared cylinder (rap-phi) distance between this jet and another |
---|
| 177 | double plain_distance(const PseudoJet & other) const; |
---|
| 178 | /// returns squared cylinder (rap-phi) distance between this jet and |
---|
| 179 | /// another |
---|
| 180 | inline double squared_distance(const PseudoJet & other) const { |
---|
| 181 | return plain_distance(other);} |
---|
| 182 | |
---|
| 183 | /// return the cylinder (rap-phi) distance between this jet and another, |
---|
| 184 | /// \f$\Delta_R = \sqrt{\Delta y^2 + \Delta \phi^2}\f$. |
---|
| 185 | inline double delta_R(const PseudoJet & other) const { |
---|
| 186 | return sqrt(squared_distance(other)); |
---|
| 187 | } |
---|
| 188 | |
---|
| 189 | /// returns other.phi() - this.phi(), constrained to be in |
---|
| 190 | /// range -pi .. pi |
---|
| 191 | double delta_phi_to(const PseudoJet & other) const; |
---|
| 192 | |
---|
| 193 | //// this seemed to compile except if it was used |
---|
| 194 | //friend inline double |
---|
| 195 | // kt_distance(const PseudoJet & jet1, const PseudoJet & jet2) { |
---|
| 196 | // return jet1.kt_distance(jet2);} |
---|
| 197 | |
---|
| 198 | /// returns distance between this jet and the beam |
---|
| 199 | inline double beam_distance() const {return _kt2;} |
---|
| 200 | |
---|
| 201 | /// return a valarray containing the four-momentum (components 0-2 |
---|
| 202 | /// are 3-mom, component 3 is energy). |
---|
| 203 | std::valarray<double> four_mom() const; |
---|
| 204 | |
---|
| 205 | //\} ------- end of kinematic access functions |
---|
| 206 | |
---|
| 207 | // taken from CLHEP |
---|
| 208 | enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES }; |
---|
| 209 | |
---|
| 210 | |
---|
| 211 | //---------------------------------------------------------------------- |
---|
| 212 | /// @name Kinematic modification functions |
---|
| 213 | //\{ |
---|
| 214 | //---------------------------------------------------------------------- |
---|
| 215 | /// transform this jet (given in the rest frame of prest) into a jet |
---|
| 216 | /// in the lab frame [NOT FULLY TESTED] |
---|
| 217 | PseudoJet & boost(const PseudoJet & prest); |
---|
| 218 | /// transform this jet (given in lab) into a jet in the rest |
---|
| 219 | /// frame of prest [NOT FULLY TESTED] |
---|
| 220 | PseudoJet & unboost(const PseudoJet & prest); |
---|
| 221 | |
---|
| 222 | void operator*=(double); |
---|
| 223 | void operator/=(double); |
---|
| 224 | void operator+=(const PseudoJet &); |
---|
| 225 | void operator-=(const PseudoJet &); |
---|
| 226 | |
---|
| 227 | /// reset the 4-momentum according to the supplied components and |
---|
| 228 | /// put the user and history indices back to their default values |
---|
| 229 | inline void reset(double px, double py, double pz, double E); |
---|
| 230 | |
---|
| 231 | /// reset the PseudoJet to be equal to psjet (including its |
---|
| 232 | /// indices); NB if the argument is derived from a PseudoJet then |
---|
| 233 | /// the "reset" used will be the templated version |
---|
| 234 | /// |
---|
| 235 | /// Note: this is included on top of the templated version because |
---|
| 236 | /// PseudoJet is not "derived" from PseudoJet, so the templated |
---|
| 237 | /// reset would not handle this case properly. |
---|
| 238 | inline void reset(const PseudoJet & psjet) { |
---|
| 239 | (*this) = psjet; |
---|
| 240 | } |
---|
| 241 | |
---|
| 242 | /// reset the 4-momentum according to the supplied generic 4-vector |
---|
| 243 | /// (accessible via indexing, [0]==px,...[3]==E) and put the user |
---|
| 244 | /// and history indices back to their default values. |
---|
| 245 | template <class L> inline void reset(const L & some_four_vector) { |
---|
| 246 | // check if some_four_vector can be cast to a PseudoJet |
---|
| 247 | // |
---|
| 248 | // Note that a regular dynamic_cast would not work here because |
---|
| 249 | // there is no guarantee that L is polymorphic. We use a more |
---|
| 250 | // complex construct here that works also in such a case. As for |
---|
| 251 | // dynamic_cast, NULL is returned if L is not derived from |
---|
| 252 | // PseudoJet |
---|
| 253 | const PseudoJet * pj = cast_if_derived<const PseudoJet>(&some_four_vector); |
---|
| 254 | |
---|
| 255 | if (pj){ |
---|
| 256 | (*this) = *pj; |
---|
| 257 | } else { |
---|
| 258 | reset(some_four_vector[0], some_four_vector[1], |
---|
| 259 | some_four_vector[2], some_four_vector[3]); |
---|
| 260 | } |
---|
| 261 | } |
---|
| 262 | |
---|
| 263 | /// reset the PseudoJet according to the specified pt, rapidity, |
---|
| 264 | /// azimuth and mass (also resetting indices, etc.) |
---|
| 265 | /// (phi should satisfy -2pi<phi<4pi) |
---|
| 266 | inline void reset_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in=0.0) { |
---|
| 267 | reset_momentum_PtYPhiM(pt_in, y_in, phi_in, m_in); |
---|
| 268 | _reset_indices(); |
---|
| 269 | } |
---|
| 270 | |
---|
| 271 | /// reset the 4-momentum according to the supplied components |
---|
| 272 | /// but leave all other information (indices, user info, etc.) |
---|
| 273 | /// untouched |
---|
| 274 | inline void reset_momentum(double px, double py, double pz, double E); |
---|
| 275 | |
---|
| 276 | /// reset the 4-momentum according to the components of the supplied |
---|
| 277 | /// PseudoJet, including cached components; note that the template |
---|
| 278 | /// version (below) will be called for classes derived from PJ. |
---|
| 279 | inline void reset_momentum(const PseudoJet & pj); |
---|
| 280 | |
---|
| 281 | /// reset the 4-momentum according to the specified pt, rapidity, |
---|
| 282 | /// azimuth and mass (phi should satisfy -2pi<phi<4pi) |
---|
| 283 | void reset_momentum_PtYPhiM(double pt, double y, double phi, double m=0.0); |
---|
| 284 | |
---|
| 285 | /// reset the 4-momentum according to the supplied generic 4-vector |
---|
| 286 | /// (accessible via indexing, [0]==px,...[3]==E), but leave all |
---|
| 287 | /// other information (indices, user info, etc.) untouched |
---|
| 288 | template <class L> inline void reset_momentum(const L & some_four_vector) { |
---|
| 289 | reset_momentum(some_four_vector[0], some_four_vector[1], |
---|
| 290 | some_four_vector[2], some_four_vector[3]); |
---|
| 291 | } |
---|
| 292 | |
---|
| 293 | /// in some cases when setting a 4-momentum, the user/program knows |
---|
| 294 | /// what rapidity and azimuth are associated with that 4-momentum; |
---|
| 295 | /// by calling this routine the user can provide the information |
---|
| 296 | /// directly to the PseudoJet and avoid expensive rap-phi |
---|
| 297 | /// recalculations. |
---|
| 298 | /// |
---|
| 299 | /// - \param rap rapidity |
---|
| 300 | /// - \param phi (in range -twopi...4*pi) |
---|
| 301 | /// |
---|
| 302 | /// USE WITH CAUTION: there are no checks that the rapidity and |
---|
| 303 | /// azimuth supplied are sensible, nor does this reset the |
---|
| 304 | /// 4-momentum components if things don't match. |
---|
| 305 | void set_cached_rap_phi(double rap, double phi); |
---|
| 306 | |
---|
| 307 | |
---|
| 308 | //\} --- end of kin mod functions ------------------------------------ |
---|
| 309 | |
---|
| 310 | //---------------------------------------------------------------------- |
---|
| 311 | /// @name User index functions |
---|
| 312 | /// |
---|
| 313 | /// To allow the user to set and access an integer index which can |
---|
| 314 | /// be exploited by the user to associate extra information with a |
---|
| 315 | /// particle/jet (for example pdg id, or an indication of a |
---|
| 316 | /// particle's origin within the user's analysis) |
---|
| 317 | // |
---|
| 318 | //\{ |
---|
| 319 | |
---|
| 320 | /// return the user_index, |
---|
| 321 | inline int user_index() const {return _user_index;} |
---|
| 322 | /// set the user_index, intended to allow the user to add simple |
---|
| 323 | /// identifying information to a particle/jet |
---|
| 324 | inline void set_user_index(const int index) {_user_index = index;} |
---|
| 325 | |
---|
| 326 | //\} ----- end of use index functions --------------------------------- |
---|
| 327 | |
---|
| 328 | //---------------------------------------------------------------------- |
---|
| 329 | /// @name User information types and functions |
---|
| 330 | /// |
---|
| 331 | /// Allows PseudoJet to carry extra user info (as an object derived from |
---|
| 332 | /// UserInfoBase). |
---|
| 333 | //\{ |
---|
| 334 | |
---|
| 335 | /// @ingroup user_info |
---|
| 336 | /// \class UserInfoBase |
---|
| 337 | /// a base class to hold extra user information in a PseudoJet |
---|
| 338 | /// |
---|
| 339 | /// This is a base class to help associate extra user information |
---|
| 340 | /// with a jet. The user should store their information in a class |
---|
| 341 | /// derived from this. This allows information of arbitrary |
---|
| 342 | /// complexity to be easily associated with a PseudoJet (in contrast |
---|
| 343 | /// to the user index). For example, in a Monte Carlo simulation, |
---|
| 344 | /// the user information might include the PDG ID, and the position |
---|
| 345 | /// of the production vertex for the particle. |
---|
| 346 | /// |
---|
| 347 | /// The PseudoJet is able to store a shared pointer to any object |
---|
| 348 | /// derived from UserInfo. The use of a shared pointer frees the |
---|
| 349 | /// user of the need to handle the memory management associated with |
---|
| 350 | /// the information. |
---|
| 351 | /// |
---|
| 352 | /// Having the user information derive from a common base class also |
---|
| 353 | /// facilitates dynamic casting, etc. |
---|
| 354 | /// |
---|
| 355 | class UserInfoBase{ |
---|
| 356 | public: |
---|
| 357 | // dummy ctor |
---|
| 358 | UserInfoBase(){}; |
---|
| 359 | |
---|
| 360 | // dummy virtual dtor |
---|
| 361 | // makes it polymorphic to allow for dynamic_cast |
---|
| 362 | virtual ~UserInfoBase(){}; |
---|
| 363 | }; |
---|
| 364 | |
---|
| 365 | /// error class to be thrown if accessing user info when it doesn't |
---|
| 366 | /// exist |
---|
| 367 | class InexistentUserInfo : public Error { |
---|
| 368 | public: |
---|
| 369 | InexistentUserInfo(); |
---|
| 370 | }; |
---|
| 371 | |
---|
| 372 | /// sets the internal shared pointer to the user information. |
---|
| 373 | /// |
---|
| 374 | /// Note that the PseudoJet will now _own_ the pointer, and delete |
---|
| 375 | /// the corresponding object when it (the jet, and any copies of the jet) |
---|
| 376 | /// goes out of scope. |
---|
| 377 | void set_user_info(UserInfoBase * user_info_in) { |
---|
| 378 | _user_info.reset(user_info_in); |
---|
| 379 | } |
---|
| 380 | |
---|
| 381 | /// returns a reference to the dynamic cast conversion of user_info |
---|
| 382 | /// to type L. |
---|
| 383 | /// |
---|
| 384 | /// Usage: suppose you have previously set the user info with a pointer |
---|
| 385 | /// to an object of type MyInfo, |
---|
| 386 | /// |
---|
| 387 | /// class MyInfo: public PseudoJet::UserInfoBase { |
---|
| 388 | /// MyInfo(int id) : _pdg_id(id); |
---|
| 389 | /// int pdg_id() const {return _pdg_id;} |
---|
| 390 | /// int _pdg_id; |
---|
| 391 | /// }; |
---|
| 392 | /// |
---|
| 393 | /// PseudoJet particle(...); |
---|
| 394 | /// particle.set_user_info(new MyInfo(its_pdg_id)); |
---|
| 395 | /// |
---|
| 396 | /// Then you would access that pdg_id() as |
---|
| 397 | /// |
---|
| 398 | /// particle.user_info<MyInfo>().pdg_id(); |
---|
| 399 | /// |
---|
| 400 | /// It's overkill for just a single integer, but scales easily to |
---|
| 401 | /// more extensive information. |
---|
| 402 | /// |
---|
| 403 | /// Note that user_info() throws an InexistentUserInfo() error if |
---|
| 404 | /// there is no user info; throws a std::bad_cast if the conversion |
---|
| 405 | /// doesn't work |
---|
| 406 | /// |
---|
| 407 | /// If this behaviour does not fit your needs, use instead the the |
---|
| 408 | /// user_info_ptr() or user_info_shared_ptr() member functions. |
---|
| 409 | template<class L> |
---|
| 410 | const L & user_info() const{ |
---|
| 411 | if (_user_info.get() == 0) throw InexistentUserInfo(); |
---|
| 412 | return dynamic_cast<const L &>(* _user_info.get()); |
---|
| 413 | } |
---|
| 414 | |
---|
| 415 | /// returns true if the PseudoJet has user information |
---|
| 416 | bool has_user_info() const{ |
---|
| 417 | return _user_info.get(); |
---|
| 418 | } |
---|
| 419 | |
---|
| 420 | /// returns true if the PseudoJet has user information than can be |
---|
| 421 | /// cast to the template argument type. |
---|
| 422 | template<class L> |
---|
| 423 | bool has_user_info() const{ |
---|
| 424 | return _user_info.get() && dynamic_cast<const L *>(_user_info.get()); |
---|
| 425 | } |
---|
| 426 | |
---|
| 427 | /// retrieve a pointer to the (const) user information |
---|
| 428 | const UserInfoBase * user_info_ptr() const{ |
---|
| 429 | if (!_user_info()) return NULL; |
---|
| 430 | return _user_info.get(); |
---|
| 431 | } |
---|
| 432 | |
---|
| 433 | |
---|
| 434 | /// retrieve a (const) shared pointer to the user information |
---|
| 435 | const SharedPtr<UserInfoBase> & user_info_shared_ptr() const{ |
---|
| 436 | return _user_info; |
---|
| 437 | } |
---|
| 438 | |
---|
| 439 | /// retrieve a (non-const) shared pointer to the user information; |
---|
| 440 | /// you can use this, for example, to set the shared pointer, eg |
---|
| 441 | /// |
---|
| 442 | /// \code |
---|
| 443 | /// p2.user_info_shared_ptr() = p1.user_info_shared_ptr(); |
---|
| 444 | /// \endcode |
---|
| 445 | /// |
---|
| 446 | /// or |
---|
| 447 | /// |
---|
| 448 | /// \code |
---|
| 449 | /// SharedPtr<PseudoJet::UserInfoBase> info_shared(new MyInfo(...)); |
---|
| 450 | /// p2.user_info_shared_ptr() = info_shared; |
---|
| 451 | /// \endcode |
---|
| 452 | SharedPtr<UserInfoBase> & user_info_shared_ptr(){ |
---|
| 453 | return _user_info; |
---|
| 454 | } |
---|
| 455 | |
---|
| 456 | // \} --- end of extra info functions --------------------------------- |
---|
| 457 | |
---|
| 458 | //---------------------------------------------------------------------- |
---|
| 459 | /// @name Description |
---|
| 460 | /// |
---|
| 461 | /// Since a PseudoJet can have a structure that contains a variety |
---|
| 462 | /// of information, we provide a description that allows one to check |
---|
| 463 | /// exactly what kind of PseudoJet we are dealing with |
---|
| 464 | // |
---|
| 465 | //\{ |
---|
| 466 | |
---|
| 467 | /// return a string describing what kind of PseudoJet we are dealing with |
---|
| 468 | std::string description() const; |
---|
| 469 | |
---|
| 470 | //\} ----- end of description functions --------------------------------- |
---|
| 471 | |
---|
| 472 | //------------------------------------------------------------- |
---|
| 473 | /// @name Access to the associated ClusterSequence object. |
---|
| 474 | /// |
---|
| 475 | /// In addition to having kinematic information, jets may contain a |
---|
| 476 | /// reference to an associated ClusterSequence (this is the case, |
---|
| 477 | /// for example, if the jet has been returned by a ClusterSequence |
---|
| 478 | /// member function). |
---|
| 479 | //\{ |
---|
| 480 | //------------------------------------------------------------- |
---|
| 481 | /// returns true if this PseudoJet has an associated ClusterSequence. |
---|
| 482 | bool has_associated_cluster_sequence() const; |
---|
| 483 | /// shorthand for has_associated_cluster_sequence() |
---|
| 484 | bool has_associated_cs() const {return has_associated_cluster_sequence();} |
---|
| 485 | |
---|
| 486 | /// returns true if this PseudoJet has an associated and still |
---|
| 487 | /// valid(ated) ClusterSequence. |
---|
| 488 | bool has_valid_cluster_sequence() const; |
---|
| 489 | /// shorthand for has_valid_cluster_sequence() |
---|
| 490 | bool has_valid_cs() const {return has_valid_cluster_sequence();} |
---|
| 491 | |
---|
| 492 | /// get a (const) pointer to the parent ClusterSequence (NULL if |
---|
| 493 | /// inexistent) |
---|
| 494 | const ClusterSequence* associated_cluster_sequence() const; |
---|
| 495 | // shorthand for associated_cluster_sequence() |
---|
| 496 | const ClusterSequence* associated_cs() const {return associated_cluster_sequence();} |
---|
| 497 | |
---|
| 498 | /// if the jet has a valid associated cluster sequence then return a |
---|
| 499 | /// pointer to it; otherwise throw an error |
---|
| 500 | inline const ClusterSequence * validated_cluster_sequence() const { |
---|
| 501 | return validated_cs(); |
---|
| 502 | } |
---|
| 503 | /// shorthand for validated_cluster_sequence() |
---|
| 504 | const ClusterSequence * validated_cs() const; |
---|
| 505 | |
---|
| 506 | /// if the jet has valid area information then return a pointer to |
---|
| 507 | /// the associated ClusterSequenceAreaBase object; otherwise throw an error |
---|
| 508 | inline const ClusterSequenceAreaBase * validated_cluster_sequence_area_base() const { |
---|
| 509 | return validated_csab(); |
---|
| 510 | } |
---|
| 511 | |
---|
| 512 | /// shorthand for validated_cluster_sequence_area_base() |
---|
| 513 | const ClusterSequenceAreaBase * validated_csab() const; |
---|
| 514 | //\} |
---|
| 515 | |
---|
| 516 | //------------------------------------------------------------- |
---|
| 517 | /// @name Access to the associated PseudoJetStructureBase object. |
---|
| 518 | /// |
---|
| 519 | /// In addition to having kinematic information, jets may contain a |
---|
| 520 | /// reference to an associated ClusterSequence (this is the case, |
---|
| 521 | /// for example, if the jet has been returned by a ClusterSequence |
---|
| 522 | /// member function). |
---|
| 523 | //\{ |
---|
| 524 | //------------------------------------------------------------- |
---|
| 525 | |
---|
| 526 | /// set the associated structure |
---|
| 527 | void set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure); |
---|
| 528 | |
---|
| 529 | /// return true if there is some structure associated with this PseudoJet |
---|
| 530 | bool has_structure() const; |
---|
| 531 | |
---|
| 532 | /// return a pointer to the structure (of type |
---|
| 533 | /// PseudoJetStructureBase*) associated with this PseudoJet. |
---|
| 534 | /// |
---|
| 535 | /// return NULL if there is no associated structure |
---|
| 536 | const PseudoJetStructureBase* structure_ptr() const; |
---|
| 537 | |
---|
| 538 | /// return a non-const pointer to the structure (of type |
---|
| 539 | /// PseudoJetStructureBase*) associated with this PseudoJet. |
---|
| 540 | /// |
---|
| 541 | /// return NULL if there is no associated structure |
---|
| 542 | /// |
---|
| 543 | /// Only use this if you know what you are doing. In any case, |
---|
| 544 | /// prefer the 'structure_ptr()' (the const version) to this method, |
---|
| 545 | /// unless you really need a write access to the PseudoJet's |
---|
| 546 | /// underlying structure. |
---|
| 547 | PseudoJetStructureBase* structure_non_const_ptr(); |
---|
| 548 | |
---|
| 549 | /// return a pointer to the structure (of type |
---|
| 550 | /// PseudoJetStructureBase*) associated with this PseudoJet. |
---|
| 551 | /// |
---|
| 552 | /// throw an error if there is no associated structure |
---|
| 553 | const PseudoJetStructureBase* validated_structure_ptr() const; |
---|
| 554 | |
---|
| 555 | /// return a reference to the shared pointer to the |
---|
| 556 | /// PseudoJetStructureBase associated with this PseudoJet |
---|
| 557 | const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const; |
---|
| 558 | |
---|
| 559 | /// returns a reference to the structure casted to the requested |
---|
| 560 | /// structure type |
---|
| 561 | /// |
---|
| 562 | /// If there is no sructure associated, an Error is thrown. |
---|
| 563 | /// If the type is not met, a std::bad_cast error is thrown. |
---|
| 564 | template<typename StructureType> |
---|
| 565 | const StructureType & structure() const; |
---|
| 566 | |
---|
| 567 | /// check if the PseudoJet has the structure resulting from a Transformer |
---|
| 568 | /// (that is, its structure is compatible with a Transformer::StructureType). |
---|
| 569 | /// If there is no structure, false is returned. |
---|
| 570 | template<typename TransformerType> |
---|
| 571 | bool has_structure_of() const; |
---|
| 572 | |
---|
| 573 | /// this is a helper to access any structure created by a Transformer |
---|
| 574 | /// (that is, of type Transformer::StructureType). |
---|
| 575 | /// |
---|
| 576 | /// If there is no structure, or if the structure is not compatible |
---|
| 577 | /// with TransformerType, an error is thrown. |
---|
| 578 | template<typename TransformerType> |
---|
| 579 | const typename TransformerType::StructureType & structure_of() const; |
---|
| 580 | |
---|
| 581 | //\} |
---|
| 582 | |
---|
| 583 | //------------------------------------------------------------- |
---|
| 584 | /// @name Methods for access to information about jet structure |
---|
| 585 | /// |
---|
| 586 | /// These allow access to jet constituents, and other jet |
---|
| 587 | /// subtructure information. They only work if the jet is associated |
---|
| 588 | /// with a ClusterSequence. |
---|
| 589 | //------------------------------------------------------------- |
---|
| 590 | //\{ |
---|
| 591 | |
---|
| 592 | /// check if it has been recombined with another PseudoJet in which |
---|
| 593 | /// case, return its partner through the argument. Otherwise, |
---|
| 594 | /// 'partner' is set to 0. |
---|
| 595 | /// |
---|
| 596 | /// an Error is thrown if this PseudoJet has no currently valid |
---|
| 597 | /// associated ClusterSequence |
---|
| 598 | virtual bool has_partner(PseudoJet &partner) const; |
---|
| 599 | |
---|
| 600 | /// check if it has been recombined with another PseudoJet in which |
---|
| 601 | /// case, return its child through the argument. Otherwise, 'child' |
---|
| 602 | /// is set to 0. |
---|
| 603 | /// |
---|
| 604 | /// an Error is thrown if this PseudoJet has no currently valid |
---|
| 605 | /// associated ClusterSequence |
---|
| 606 | virtual bool has_child(PseudoJet &child) const; |
---|
| 607 | |
---|
| 608 | /// check if it is the product of a recombination, in which case |
---|
| 609 | /// return the 2 parents through the 'parent1' and 'parent2' |
---|
| 610 | /// arguments. Otherwise, set these to 0. |
---|
| 611 | /// |
---|
| 612 | /// an Error is thrown if this PseudoJet has no currently valid |
---|
| 613 | /// associated ClusterSequence |
---|
| 614 | virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const; |
---|
| 615 | |
---|
| 616 | /// check if the current PseudoJet contains the one passed as |
---|
| 617 | /// argument. |
---|
| 618 | /// |
---|
| 619 | /// an Error is thrown if this PseudoJet has no currently valid |
---|
| 620 | /// associated ClusterSequence |
---|
| 621 | virtual bool contains(const PseudoJet &constituent) const; |
---|
| 622 | |
---|
| 623 | /// check if the current PseudoJet is contained the one passed as |
---|
| 624 | /// argument. |
---|
| 625 | /// |
---|
| 626 | /// an Error is thrown if this PseudoJet has no currently valid |
---|
| 627 | /// associated ClusterSequence |
---|
| 628 | virtual bool is_inside(const PseudoJet &jet) const; |
---|
| 629 | |
---|
| 630 | |
---|
| 631 | /// returns true if the PseudoJet has constituents |
---|
| 632 | virtual bool has_constituents() const; |
---|
| 633 | |
---|
| 634 | /// retrieve the constituents. |
---|
| 635 | /// |
---|
| 636 | /// an Error is thrown if this PseudoJet has no currently valid |
---|
| 637 | /// associated ClusterSequence or other substructure information |
---|
| 638 | virtual std::vector<PseudoJet> constituents() const; |
---|
| 639 | |
---|
| 640 | |
---|
| 641 | /// returns true if the PseudoJet has support for exclusive subjets |
---|
| 642 | virtual bool has_exclusive_subjets() const; |
---|
| 643 | |
---|
| 644 | /// return a vector of all subjets of the current jet (in the sense |
---|
| 645 | /// of the exclusive algorithm) that would be obtained when running |
---|
| 646 | /// the algorithm with the given dcut. |
---|
| 647 | /// |
---|
| 648 | /// Time taken is O(m ln m), where m is the number of subjets that |
---|
| 649 | /// are found. If m gets to be of order of the total number of |
---|
| 650 | /// constituents in the jet, this could be substantially slower than |
---|
| 651 | /// just getting that list of constituents. |
---|
| 652 | /// |
---|
| 653 | /// an Error is thrown if this PseudoJet has no currently valid |
---|
| 654 | /// associated ClusterSequence |
---|
| 655 | std::vector<PseudoJet> exclusive_subjets (const double & dcut) const; |
---|
| 656 | |
---|
| 657 | /// return the size of exclusive_subjets(...); still n ln n with same |
---|
| 658 | /// coefficient, but marginally more efficient than manually taking |
---|
| 659 | /// exclusive_subjets.size() |
---|
| 660 | /// |
---|
| 661 | /// an Error is thrown if this PseudoJet has no currently valid |
---|
| 662 | /// associated ClusterSequence |
---|
| 663 | int n_exclusive_subjets(const double & dcut) const; |
---|
| 664 | |
---|
| 665 | /// return the list of subjets obtained by unclustering the supplied |
---|
| 666 | /// jet down to nsub subjets. Throws an error if there are fewer than |
---|
| 667 | /// nsub particles in the jet. |
---|
| 668 | /// |
---|
| 669 | /// For ClusterSequence type jets, requires nsub ln nsub time |
---|
| 670 | /// |
---|
| 671 | /// An Error is thrown if this PseudoJet has no currently valid |
---|
| 672 | /// associated ClusterSequence |
---|
| 673 | std::vector<PseudoJet> exclusive_subjets (int nsub) const; |
---|
| 674 | |
---|
| 675 | /// return the list of subjets obtained by unclustering the supplied |
---|
| 676 | /// jet down to nsub subjets (or all constituents if there are fewer |
---|
| 677 | /// than nsub). |
---|
| 678 | /// |
---|
| 679 | /// For ClusterSequence type jets, requires nsub ln nsub time |
---|
| 680 | /// |
---|
| 681 | /// An Error is thrown if this PseudoJet has no currently valid |
---|
| 682 | /// associated ClusterSequence |
---|
| 683 | std::vector<PseudoJet> exclusive_subjets_up_to (int nsub) const; |
---|
| 684 | |
---|
| 685 | /// return the dij that was present in the merging nsub+1 -> nsub |
---|
| 686 | /// subjets inside this jet. |
---|
| 687 | /// |
---|
| 688 | /// an Error is thrown if this PseudoJet has no currently valid |
---|
| 689 | /// associated ClusterSequence |
---|
| 690 | double exclusive_subdmerge(int nsub) const; |
---|
| 691 | |
---|
| 692 | /// return the maximum dij that occurred in the whole event at the |
---|
| 693 | /// stage that the nsub+1 -> nsub merge of subjets occurred inside |
---|
| 694 | /// this jet. |
---|
| 695 | /// |
---|
| 696 | /// an Error is thrown if this PseudoJet has no currently valid |
---|
| 697 | /// associated ClusterSequence |
---|
| 698 | double exclusive_subdmerge_max(int nsub) const; |
---|
| 699 | |
---|
| 700 | |
---|
| 701 | /// returns true if a jet has pieces |
---|
| 702 | /// |
---|
| 703 | /// By default a single particle or a jet coming from a |
---|
| 704 | /// ClusterSequence have no pieces and this methos will return false. |
---|
| 705 | /// |
---|
| 706 | /// In practice, this is equivalent to have an structure of type |
---|
| 707 | /// CompositeJetStructure. |
---|
| 708 | virtual bool has_pieces() const; |
---|
| 709 | |
---|
| 710 | |
---|
| 711 | /// retrieve the pieces that make up the jet. |
---|
| 712 | /// |
---|
| 713 | /// If the jet does not support pieces, an error is throw |
---|
| 714 | virtual std::vector<PseudoJet> pieces() const; |
---|
| 715 | |
---|
| 716 | |
---|
| 717 | // the following ones require a computation of the area in the |
---|
| 718 | // parent ClusterSequence (See ClusterSequenceAreaBase for details) |
---|
| 719 | //------------------------------------------------------------------ |
---|
| 720 | |
---|
| 721 | /// check if it has a defined area |
---|
| 722 | virtual bool has_area() const; |
---|
| 723 | |
---|
| 724 | /// return the jet (scalar) area. |
---|
| 725 | /// throws an Error if there is no support for area in the parent CS |
---|
| 726 | virtual double area() const; |
---|
| 727 | |
---|
| 728 | /// return the error (uncertainty) associated with the determination |
---|
| 729 | /// of the area of this jet. |
---|
| 730 | /// throws an Error if there is no support for area in the parent CS |
---|
| 731 | virtual double area_error() const; |
---|
| 732 | |
---|
| 733 | /// return the jet 4-vector area. |
---|
| 734 | /// throws an Error if there is no support for area in the parent CS |
---|
| 735 | virtual PseudoJet area_4vector() const; |
---|
| 736 | |
---|
| 737 | /// true if this jet is made exclusively of ghosts. |
---|
| 738 | /// throws an Error if there is no support for area in the parent CS |
---|
| 739 | virtual bool is_pure_ghost() const; |
---|
| 740 | |
---|
| 741 | //\} --- end of jet structure ------------------------------------- |
---|
| 742 | |
---|
| 743 | |
---|
| 744 | |
---|
| 745 | //---------------------------------------------------------------------- |
---|
| 746 | /// @name Members mainly intended for internal use |
---|
| 747 | //---------------------------------------------------------------------- |
---|
| 748 | //\{ |
---|
| 749 | /// return the cluster_hist_index, intended to be used by clustering |
---|
| 750 | /// routines. |
---|
| 751 | inline int cluster_hist_index() const {return _cluster_hist_index;} |
---|
| 752 | /// set the cluster_hist_index, intended to be used by clustering routines. |
---|
| 753 | inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;} |
---|
| 754 | |
---|
| 755 | /// alternative name for cluster_hist_index() [perhaps more meaningful] |
---|
| 756 | inline int cluster_sequence_history_index() const { |
---|
| 757 | return cluster_hist_index();} |
---|
| 758 | /// alternative name for set_cluster_hist_index(...) [perhaps more |
---|
| 759 | /// meaningful] |
---|
| 760 | inline void set_cluster_sequence_history_index(const int index) { |
---|
| 761 | set_cluster_hist_index(index);} |
---|
| 762 | |
---|
| 763 | //\} ---- end of internal use functions --------------------------- |
---|
| 764 | |
---|
| 765 | protected: |
---|
| 766 | |
---|
| 767 | SharedPtr<PseudoJetStructureBase> _structure; |
---|
| 768 | SharedPtr<UserInfoBase> _user_info; |
---|
| 769 | |
---|
| 770 | |
---|
| 771 | private: |
---|
| 772 | // NB: following order must be kept for things to behave sensibly... |
---|
| 773 | double _px,_py,_pz,_E; |
---|
| 774 | mutable double _phi, _rap; |
---|
| 775 | double _kt2; |
---|
| 776 | int _cluster_hist_index, _user_index; |
---|
| 777 | |
---|
| 778 | /// calculate phi, rap, kt2 based on the 4-momentum components |
---|
| 779 | void _finish_init(); |
---|
| 780 | /// set the indices to default values |
---|
| 781 | void _reset_indices(); |
---|
| 782 | |
---|
| 783 | /// ensure that the internal values for rapidity and phi |
---|
| 784 | /// correspond to 4-momentum structure |
---|
| 785 | inline void _ensure_valid_rap_phi() const { |
---|
| 786 | if (_phi == pseudojet_invalid_phi) _set_rap_phi(); |
---|
| 787 | } |
---|
| 788 | |
---|
| 789 | /// set cached rapidity and phi values |
---|
| 790 | void _set_rap_phi() const; |
---|
| 791 | }; |
---|
| 792 | |
---|
| 793 | |
---|
| 794 | //---------------------------------------------------------------------- |
---|
| 795 | // routines for basic binary operations |
---|
| 796 | |
---|
| 797 | PseudoJet operator+(const PseudoJet &, const PseudoJet &); |
---|
| 798 | PseudoJet operator-(const PseudoJet &, const PseudoJet &); |
---|
| 799 | PseudoJet operator*(double, const PseudoJet &); |
---|
| 800 | PseudoJet operator*(const PseudoJet &, double); |
---|
| 801 | PseudoJet operator/(const PseudoJet &, double); |
---|
| 802 | |
---|
| 803 | /// returns true if the 4 momentum components of the two PseudoJets |
---|
| 804 | /// are identical and all the internal indices (user, cluster_history) |
---|
| 805 | /// + structure and user-info shared pointers are too |
---|
| 806 | bool operator==(const PseudoJet &, const PseudoJet &); |
---|
| 807 | |
---|
| 808 | /// inequality test which is exact opposite of operator== |
---|
| 809 | inline bool operator!=(const PseudoJet & a, const PseudoJet & b) {return !(a==b);} |
---|
| 810 | |
---|
| 811 | /// Can only be used with val=0 and tests whether all four |
---|
| 812 | /// momentum components are equal to val (=0.0) |
---|
| 813 | bool operator==(const PseudoJet & jet, const double val); |
---|
| 814 | |
---|
| 815 | /// Can only be used with val=0 and tests whether at least one of the |
---|
| 816 | /// four momentum components is different from val (=0.0) |
---|
| 817 | inline bool operator!=(const PseudoJet & a, const double & val) {return !(a==val);} |
---|
| 818 | |
---|
| 819 | inline double dot_product(const PseudoJet & a, const PseudoJet & b) { |
---|
| 820 | return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz(); |
---|
| 821 | } |
---|
| 822 | |
---|
| 823 | /// returns true if the momenta of the two input jets are identical |
---|
| 824 | bool have_same_momentum(const PseudoJet &, const PseudoJet &); |
---|
| 825 | |
---|
| 826 | /// return a pseudojet with the given pt, y, phi and mass |
---|
| 827 | /// (phi should satisfy -2pi<phi<4pi) |
---|
| 828 | PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0); |
---|
| 829 | |
---|
| 830 | //---------------------------------------------------------------------- |
---|
| 831 | // Routines to do with providing sorted arrays of vectors. |
---|
| 832 | |
---|
| 833 | /// return a vector of jets sorted into decreasing transverse momentum |
---|
| 834 | std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets); |
---|
| 835 | |
---|
| 836 | /// return a vector of jets sorted into increasing rapidity |
---|
| 837 | std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets); |
---|
| 838 | |
---|
| 839 | /// return a vector of jets sorted into decreasing energy |
---|
| 840 | std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets); |
---|
| 841 | |
---|
| 842 | /// return a vector of jets sorted into increasing pz |
---|
| 843 | std::vector<PseudoJet> sorted_by_pz(const std::vector<PseudoJet> & jets); |
---|
| 844 | |
---|
| 845 | //---------------------------------------------------------------------- |
---|
| 846 | // some code to help sorting |
---|
| 847 | |
---|
| 848 | /// sort the indices so that values[indices[0->n-1]] is sorted |
---|
| 849 | /// into increasing order |
---|
| 850 | void sort_indices(std::vector<int> & indices, |
---|
| 851 | const std::vector<double> & values); |
---|
| 852 | |
---|
| 853 | /// given a vector of values with a one-to-one correspondence with the |
---|
| 854 | /// vector of objects, sort objects into an order such that the |
---|
| 855 | /// associated values would be in increasing order (but don't actually |
---|
| 856 | /// touch the values vector in the process). |
---|
| 857 | template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects, |
---|
| 858 | const std::vector<double> & values); |
---|
| 859 | |
---|
| 860 | /// \if internal_doc |
---|
| 861 | /// @ingroup internal |
---|
| 862 | /// \class IndexedSortHelper |
---|
| 863 | /// a class that helps us carry out indexed sorting. |
---|
| 864 | /// \endif |
---|
| 865 | class IndexedSortHelper { |
---|
| 866 | public: |
---|
| 867 | inline IndexedSortHelper (const std::vector<double> * reference_values) { |
---|
| 868 | _ref_values = reference_values; |
---|
| 869 | }; |
---|
| 870 | inline int operator() (const int & i1, const int & i2) const { |
---|
| 871 | return (*_ref_values)[i1] < (*_ref_values)[i2]; |
---|
| 872 | }; |
---|
| 873 | private: |
---|
| 874 | const std::vector<double> * _ref_values; |
---|
| 875 | }; |
---|
| 876 | |
---|
| 877 | |
---|
| 878 | //---------------------------------------------------------------------- |
---|
| 879 | /// constructor from any object that has px,py,pz,E = some_four_vector[0--3], |
---|
| 880 | // NB: do not know if it really needs to be inline, but when it wasn't |
---|
| 881 | // linking failed with g++ (who knows what was wrong...) |
---|
| 882 | template <class L> inline PseudoJet::PseudoJet(const L & some_four_vector) { |
---|
| 883 | reset(some_four_vector); |
---|
| 884 | } |
---|
| 885 | |
---|
| 886 | //---------------------------------------------------------------------- |
---|
| 887 | inline void PseudoJet::_reset_indices() { |
---|
| 888 | set_cluster_hist_index(-1); |
---|
| 889 | set_user_index(-1); |
---|
| 890 | _structure.reset(); |
---|
| 891 | _user_info.reset(); |
---|
| 892 | } |
---|
| 893 | |
---|
| 894 | |
---|
| 895 | // taken literally from CLHEP |
---|
| 896 | inline double PseudoJet::m() const { |
---|
| 897 | double mm = m2(); |
---|
| 898 | return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm); |
---|
| 899 | } |
---|
| 900 | |
---|
| 901 | |
---|
| 902 | inline void PseudoJet::reset(double px_in, double py_in, double pz_in, double E_in) { |
---|
| 903 | _px = px_in; |
---|
| 904 | _py = py_in; |
---|
| 905 | _pz = pz_in; |
---|
| 906 | _E = E_in; |
---|
| 907 | _finish_init(); |
---|
| 908 | _reset_indices(); |
---|
| 909 | } |
---|
| 910 | |
---|
| 911 | inline void PseudoJet::reset_momentum(double px_in, double py_in, double pz_in, double E_in) { |
---|
| 912 | _px = px_in; |
---|
| 913 | _py = py_in; |
---|
| 914 | _pz = pz_in; |
---|
| 915 | _E = E_in; |
---|
| 916 | _finish_init(); |
---|
| 917 | } |
---|
| 918 | |
---|
| 919 | inline void PseudoJet::reset_momentum(const PseudoJet & pj) { |
---|
| 920 | _px = pj._px ; |
---|
| 921 | _py = pj._py ; |
---|
| 922 | _pz = pj._pz ; |
---|
| 923 | _E = pj._E ; |
---|
| 924 | _phi = pj._phi; |
---|
| 925 | _rap = pj._rap; |
---|
| 926 | _kt2 = pj._kt2; |
---|
| 927 | } |
---|
| 928 | |
---|
| 929 | //------------------------------------------------------------------------------- |
---|
| 930 | // implementation of the templated accesses to the underlying structyre |
---|
| 931 | //------------------------------------------------------------------------------- |
---|
| 932 | |
---|
| 933 | // returns a reference to the structure casted to the requested |
---|
| 934 | // structure type |
---|
| 935 | // |
---|
| 936 | // If there is no sructure associated, an Error is thrown. |
---|
| 937 | // If the type is not met, a std::bad_cast error is thrown. |
---|
| 938 | template<typename StructureType> |
---|
| 939 | const StructureType & PseudoJet::structure() const{ |
---|
| 940 | return dynamic_cast<const StructureType &>(* validated_structure_ptr()); |
---|
| 941 | |
---|
| 942 | } |
---|
| 943 | |
---|
| 944 | // check if the PseudoJet has the structure resulting from a Transformer |
---|
| 945 | // (that is, its structure is compatible with a Transformer::StructureType) |
---|
| 946 | template<typename TransformerType> |
---|
| 947 | bool PseudoJet::has_structure_of() const{ |
---|
| 948 | if (!_structure()) return false; |
---|
| 949 | |
---|
| 950 | return dynamic_cast<const typename TransformerType::StructureType *>(_structure.get()) != 0; |
---|
| 951 | } |
---|
| 952 | |
---|
| 953 | // this is a helper to access a structure created by a Transformer |
---|
| 954 | // (that is, of type Transformer::StructureType) |
---|
| 955 | // NULL is returned if the corresponding type is not met |
---|
| 956 | template<typename TransformerType> |
---|
| 957 | const typename TransformerType::StructureType & PseudoJet::structure_of() const{ |
---|
| 958 | if (!_structure()) |
---|
| 959 | throw Error("Trying to access the structure of a PseudoJet without an associated structure"); |
---|
| 960 | |
---|
| 961 | return dynamic_cast<const typename TransformerType::StructureType &>(*_structure); |
---|
| 962 | } |
---|
| 963 | |
---|
| 964 | |
---|
| 965 | |
---|
| 966 | //------------------------------------------------------------------------------- |
---|
| 967 | // helper functions to build a jet made of pieces |
---|
| 968 | // |
---|
| 969 | // Note that there are more complete versions of these functions, with |
---|
| 970 | // an additional argument for a recombination scheme, in |
---|
| 971 | // JetDefinition.hh |
---|
| 972 | // ------------------------------------------------------------------------------- |
---|
| 973 | |
---|
| 974 | /// build a "CompositeJet" from the vector of its pieces |
---|
| 975 | /// |
---|
| 976 | /// In this case, E-scheme recombination is assumed to compute the |
---|
| 977 | /// total momentum |
---|
| 978 | PseudoJet join(const std::vector<PseudoJet> & pieces); |
---|
| 979 | |
---|
| 980 | /// build a MergedJet from a single PseudoJet |
---|
| 981 | PseudoJet join(const PseudoJet & j1); |
---|
| 982 | |
---|
| 983 | /// build a MergedJet from 2 PseudoJet |
---|
| 984 | PseudoJet join(const PseudoJet & j1, const PseudoJet & j2); |
---|
| 985 | |
---|
| 986 | /// build a MergedJet from 3 PseudoJet |
---|
| 987 | PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3); |
---|
| 988 | |
---|
| 989 | /// build a MergedJet from 4 PseudoJet |
---|
| 990 | PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4); |
---|
| 991 | |
---|
| 992 | |
---|
| 993 | |
---|
| 994 | FASTJET_END_NAMESPACE |
---|
| 995 | |
---|
| 996 | #endif // __FASTJET_PSEUDOJET_HH__ |
---|