Fork me on GitHub

source: git/external/fastjet/PseudoJet.hh@ 0e6b14e

ImprovedOutputFile Timing
Last change on this file since 0e6b14e was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 7 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

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