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