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