Fork me on GitHub

source: svn/trunk/external/fastjet/PseudoJet.hh@ 1397

Last change on this file since 1397 was 1332, checked in by Pavel Demin, 11 years ago

upgrade fastjet to 3.0.6

File size: 37.1 KB
Line 
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
44FASTJET_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.
50const double MaxRap = 1e5;
51
52/// default value for phi, meaning it (and rapidity) have yet to be calculated)
53const double pseudojet_invalid_phi = -100.0;
54const double pseudojet_invalid_rap = -1e200;
55
56#ifndef __FJCORE__
57// forward definition
58class ClusterSequenceAreaBase;
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.
65class 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 ///
363 class UserInfoBase{
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:
377 InexistentUserInfo();
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
443 const SharedPtr<UserInfoBase> & user_info_shared_ptr() const{
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
460 SharedPtr<UserInfoBase> & user_info_shared_ptr(){
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
508 inline const ClusterSequence * validated_cluster_sequence() const {
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
517 inline const ClusterSequenceAreaBase * validated_cluster_sequence_area_base() const {
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
780 SharedPtr<PseudoJetStructureBase> _structure;
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
810PseudoJet operator+(const PseudoJet &, const PseudoJet &);
811PseudoJet operator-(const PseudoJet &, const PseudoJet &);
812PseudoJet operator*(double, const PseudoJet &);
813PseudoJet operator*(const PseudoJet &, double);
814PseudoJet 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
819bool operator==(const PseudoJet &, const PseudoJet &);
820
821/// inequality test which is exact opposite of operator==
822inline 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)
826bool 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)
830inline bool operator!=(const PseudoJet & a, const double & val) {return !(a==val);}
831
832inline 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
837bool 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)
841PseudoJet 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
847std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets);
848
849/// return a vector of jets sorted into increasing rapidity
850std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets);
851
852/// return a vector of jets sorted into decreasing energy
853std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets);
854
855/// return a vector of jets sorted into increasing pz
856std::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
863void 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).
870template<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
878class IndexedSortHelper {
879public:
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 };
886private:
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...)
895template <class L> inline PseudoJet::PseudoJet(const L & some_four_vector) {
896 reset(some_four_vector);
897}
898
899//----------------------------------------------------------------------
900inline 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
909inline double PseudoJet::m() const {
910 double mm = m2();
911 return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
912}
913
914
915inline 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
924inline 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
932inline 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.
951template<typename StructureType>
952const 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)
959template<typename TransformerType>
960bool 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
969template<typename TransformerType>
970const 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
991PseudoJet join(const std::vector<PseudoJet> & pieces);
992
993/// build a MergedJet from a single PseudoJet
994PseudoJet join(const PseudoJet & j1);
995
996/// build a MergedJet from 2 PseudoJet
997PseudoJet join(const PseudoJet & j1, const PseudoJet & j2);
998
999/// build a MergedJet from 3 PseudoJet
1000PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3);
1001
1002/// build a MergedJet from 4 PseudoJet
1003PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4);
1004
1005
1006
1007FASTJET_END_NAMESPACE
1008
1009#endif // __FASTJET_PSEUDOJET_HH__
Note: See TracBrowser for help on using the repository browser.