Fork me on GitHub

source: git/external/fastjet/PseudoJet.hh@ a37e592

ImprovedOutputFile Timing dual_readout llp
Last change on this file since a37e592 was 1d208a2, checked in by Pavel Demin <pavel.demin@…>, 8 years ago

update FastJet library to 3.2.1 and Nsubjettiness library to 2.2.4

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