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