FastJet 3.4.0
PseudoJet.cc
1//FJSTARTHEADER
2// $Id$
3//
4// Copyright (c) 2005-2021, 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#include "fastjet/Error.hh"
33#include "fastjet/PseudoJet.hh"
34#include "fastjet/ClusterSequence.hh"
35#ifndef __FJCORE__
36#include "fastjet/ClusterSequenceAreaBase.hh"
37#endif // __FJCORE__
38#include "fastjet/CompositeJetStructure.hh"
39#include<valarray>
40#include<iostream>
41#include<sstream>
42#include<cmath>
43#include<algorithm>
44#include <cstdarg>
45
46FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
47
48using namespace std;
49
50
51//----------------------------------------------------------------------
52// another constructor...
53PseudoJet::PseudoJet(const double px_in, const double py_in, const double pz_in, const double E_in) {
54
55 _E = E_in ;
56 _px = px_in;
57 _py = py_in;
58 _pz = pz_in;
59
60 this->_finish_init();
61
62 // some default values for the history and user indices
63 // note: no reset of shared pointers needed
64 _reset_indices();
65
66}
67
68#ifdef FASTJET_HAVE_THREAD_SAFETY
69/// copy-assignmemt
70///
71/// this has to be explicitly specified since atomic does not support it.
72PseudoJet & PseudoJet::operator=(const PseudoJet & other_pj){
73 _structure = other_pj._structure;
74 _user_info = other_pj._user_info;
75
76 _kt2 = other_pj._kt2;
77 _cluster_hist_index = other_pj._cluster_hist_index;
78 _user_index = other_pj._user_index;
79
80 _px = other_pj._px;
81 _py = other_pj._py;
82 _pz = other_pj._pz;
83 _E = other_pj._E;
84
85 _phi = other_pj._phi;
86 _rap = other_pj._rap;
87
88 _init_status.store(other_pj._init_status);
89
90 return *this;
91}
92#endif // FASTJET_HAVE_THREAD_SAFETY
93
94//std::shared_ptr: //----------------------------------------------------------------------
95//std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
96//std::shared_ptr: PseudoJet::~PseudoJet(){
97//std::shared_ptr: _release_jet_from_cs();
98//std::shared_ptr: }
99//std::shared_ptr:
100//std::shared_ptr: // this has to be called everytime one tries to alter the jet
101//std::shared_ptr: // structural info
102//std::shared_ptr: void PseudoJet::_release_jet_from_cs(){
103//std::shared_ptr: // check if the jet has the structure of type CSstruct in which case
104//std::shared_ptr: // we have to check if there is a need for self-deletion of the CS
105//std::shared_ptr: ClusterSequenceStructure * assoc_css = dynamic_cast<ClusterSequenceStructure *>(_structure.get());
106//std::shared_ptr: if (assoc_css) {
107//std::shared_ptr: const ClusterSequence * assoc_cs = assoc_css->associated_cluster_sequence();
108//std::shared_ptr: if (assoc_cs) assoc_cs->release_pseudojet(*this);
109//std::shared_ptr: }
110//std::shared_ptr:
111//std::shared_ptr: // slightly less efficient version
112//std::shared_ptr: // if ((has_structure_of<ClusterSequence>()) && (has_valid_cluster_sequence())){
113//std::shared_ptr: // associated_cs()->release_pseudojet(*this);
114//std::shared_ptr: // }
115//std::shared_ptr: }
116//std::shared_ptr:
117//std::shared_ptr: #endif // FASTJET_HAVE_THREAD_SAFETY
118
119//----------------------------------------------------------------------
120/// do standard end of initialisation
121void PseudoJet::_finish_init () {
122 _kt2 = this->px()*this->px() + this->py()*this->py();
123 _phi = pseudojet_invalid_phi;
124 // strictly speaking, _rap does not need initialising, because
125 // it's never used as long as _phi == pseudojet_invalid_phi
126 // (and gets set when _phi is requested). However ATLAS
127 // 2013-03-28 complained that they sometimes have NaN's in
128 // _rap and this interferes with some of their internal validation.
129 // So we initialise it; penalty is about 0.3ns per PseudoJet out of about
130 // 10ns total initialisation time (on a intel Core i7 2.7GHz)
131 _rap = pseudojet_invalid_rap;
132
133#ifdef FASTJET_HAVE_THREAD_SAFETY
134 _init_status = Init_NotDone;
135#endif
136}
137
138//----------------------------------------------------------------------
139#ifdef FASTJET_HAVE_THREAD_SAFETY
140void PseudoJet::_ensure_valid_rap_phi() const{
141 //TODO: for _init_status we can use memory_order_release when
142 //writing and memory_order_acquire when reading. Check if that has
143 //an impact on timing
144 if (_init_status!=Init_Done){
145 int expected = Init_NotDone;
146 // if it's 0 (not done), set thing to "operation in progress"
147 // (-1) and do the init
148 //
149 // Note:
150 // - this has to be the strong version so we can make a single
151 // test
152 // - we need to make sure that the -1 is loaded properly
153 // (success), hence a std::memory_order_seq_cst model
154 // - the failure model can be anything since we're not going
155 // to use the value anyway
156 if (_init_status.compare_exchange_strong(expected, Init_InProgress,
157 std::memory_order_seq_cst,
158 std::memory_order_relaxed)){
159 _set_rap_phi();
160 _init_status = Init_Done; // can safely be done after all physics varlables are set
161 } else {
162 // wait until done
163 do{
164 // the operation below will reset expected to whatever is in
165 // init_state if the test fails, so we need to reset it to
166 // the 1 (aka init_done) we want!
167 expected = Init_Done;
168
169 // the next line
170 // - here we could potentially use the weak form
171 // - on success the value is unchanged so I think we can use relaxed ordering
172 // - expected will be reinitialised anyway so again, relaxed ordering should be fi
173
174 //} while (!_init_state.compare_exchange_strong(expected, 1));
175 } while (!_init_status.compare_exchange_weak(expected, Init_Done,
176 std::memory_order_relaxed,
177 std::memory_order_relaxed));
178 }
179
180 }
181}
182#endif
183
184//----------------------------------------------------------------------
185void PseudoJet::_set_rap_phi() const {
186
187 if (_kt2 == 0.0) {
188 _phi = 0.0; }
189 else {
190 _phi = atan2(this->py(),this->px());
191 }
192 if (_phi < 0.0) {_phi += twopi;}
193 if (_phi >= twopi) {_phi -= twopi;} // can happen if phi=-|eps<1e-15|?
194 if (this->E() == abs(this->pz()) && _kt2 == 0) {
195 // Point has infinite rapidity -- convert that into a very large
196 // number, but in such a way that different 0-pt momenta will have
197 // different rapidities (so as to lift the degeneracy between
198 // them) [this can be relevant at parton-level]
199 double MaxRapHere = MaxRap + abs(this->pz());
200 if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;}
201 } else {
202 // get the rapidity in a way that's modestly insensitive to roundoff
203 // error when things pz,E are large (actually the best we can do without
204 // explicit knowledge of mass)
205 double effective_m2 = max(0.0,m2()); // force non tachyonic mass
206 double E_plus_pz = _E + abs(_pz); // the safer of p+, p-
207 // p+/p- = (p+ p-) / (p-)^2 = (kt^2+m^2)/(p-)^2
208 _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
209 if (_pz > 0) {_rap = - _rap;}
210 }
211
212}
213
214
215//----------------------------------------------------------------------
216// return a valarray four-momentum
217valarray<double> PseudoJet::four_mom() const {
218 valarray<double> mom(4);
219 mom[0] = _px;
220 mom[1] = _py;
221 mom[2] = _pz;
222 mom[3] = _E ;
223 return mom;
224}
225
226//----------------------------------------------------------------------
227// Return the component corresponding to the specified index.
228// taken from CLHEP
229double PseudoJet::operator () (int i) const {
230 switch(i) {
231 case X:
232 return px();
233 case Y:
234 return py();
235 case Z:
236 return pz();
237 case T:
238 return e();
239 default:
240 ostringstream err;
241 err << "PseudoJet subscripting: bad index (" << i << ")";
242 throw Error(err.str());
243 }
244 return 0.;
245}
246
247//----------------------------------------------------------------------
248// return the pseudorapidity
249double PseudoJet::pseudorapidity() const {
250 if (px() == 0.0 && py() ==0.0) return MaxRap;
251 if (pz() == 0.0) return 0.0;
252
253 double theta = atan(perp()/pz());
254 if (theta < 0) theta += pi;
255 return -log(tan(theta/2));
256}
257
258//----------------------------------------------------------------------
259// return "sum" of two pseudojets
260PseudoJet operator+ (const PseudoJet & jet1, const PseudoJet & jet2) {
261 //return PseudoJet(jet1.four_mom()+jet2.four_mom());
262 return PseudoJet(jet1.px()+jet2.px(),
263 jet1.py()+jet2.py(),
264 jet1.pz()+jet2.pz(),
265 jet1.E() +jet2.E() );
266}
267
268//----------------------------------------------------------------------
269// return difference of two pseudojets
270PseudoJet operator- (const PseudoJet & jet1, const PseudoJet & jet2) {
271 //return PseudoJet(jet1.four_mom()-jet2.four_mom());
272 return PseudoJet(jet1.px()-jet2.px(),
273 jet1.py()-jet2.py(),
274 jet1.pz()-jet2.pz(),
275 jet1.E() -jet2.E() );
276}
277
278//----------------------------------------------------------------------
279// return the product, coeff * jet
280PseudoJet operator* (double coeff, const PseudoJet & jet) {
281 // see the comment in operator*= about ensuring valid rap phi
282 // before a multiplication to handle case of multiplication by
283 // zero, while maintaining rapidity and phi
284 jet._ensure_valid_rap_phi();
285 //return PseudoJet(coeff*jet.four_mom());
286 // the following code is hopefully more efficient
287 PseudoJet coeff_times_jet(jet);
288 coeff_times_jet *= coeff;
289 return coeff_times_jet;
290}
291
292//----------------------------------------------------------------------
293// return the product, coeff * jet
294PseudoJet operator* (const PseudoJet & jet, double coeff) {
295 return coeff*jet;
296}
297
298//----------------------------------------------------------------------
299// return the ratio, jet / coeff
300PseudoJet operator/ (const PseudoJet & jet, double coeff) {
301 return (1.0/coeff)*jet;
302}
303
304//----------------------------------------------------------------------
305/// multiply the jet's momentum by the coefficient
306PseudoJet & PseudoJet::operator*=(double coeff) {
307 // operator*= aims to maintain the rapidity and azimuth
308 // for the PseudoJet; if they have already been evaluated
309 // this is fine, but if they haven't and coeff is sufficiently
310 // small as to cause a zero or underflow result, then a subsequent
311 // invocation of rap or phi will lead to a non-sensical result.
312 // So, here, we preemptively ensure that rapidity and phi
313 // are correctly cached
314 _ensure_valid_rap_phi();
315 _px *= coeff;
316 _py *= coeff;
317 _pz *= coeff;
318 _E *= coeff;
319 _kt2*= coeff*coeff;
320 // phi and rap are unchanged
321 return *this;
322}
323
324//----------------------------------------------------------------------
325/// divide the jet's momentum by the coefficient
326PseudoJet & PseudoJet::operator/=(double coeff) {
327 (*this) *= 1.0/coeff;
328 return *this;
329}
330
331
332//----------------------------------------------------------------------
333/// add the other jet's momentum to this jet
334PseudoJet & PseudoJet::operator+=(const PseudoJet & other_jet) {
335 _px += other_jet._px;
336 _py += other_jet._py;
337 _pz += other_jet._pz;
338 _E += other_jet._E ;
339 _finish_init(); // we need to recalculate phi,rap,kt2
340 return *this;
341}
342
343
344//----------------------------------------------------------------------
345/// subtract the other jet's momentum from this jet
346PseudoJet & PseudoJet::operator-=(const PseudoJet & other_jet) {
347 _px -= other_jet._px;
348 _py -= other_jet._py;
349 _pz -= other_jet._pz;
350 _E -= other_jet._E ;
351 _finish_init(); // we need to recalculate phi,rap,kt2
352 return *this;
353}
354
355//----------------------------------------------------------------------
356bool operator==(const PseudoJet & a, const PseudoJet & b) {
357 if (a.px() != b.px()) return false;
358 if (a.py() != b.py()) return false;
359 if (a.pz() != b.pz()) return false;
360 if (a.E () != b.E ()) return false;
361
362 if (a.user_index() != b.user_index()) return false;
363 if (a.cluster_hist_index() != b.cluster_hist_index()) return false;
364 if (a.user_info_ptr() != b.user_info_ptr()) return false;
365 if (a.structure_ptr() != b.structure_ptr()) return false;
366
367 return true;
368}
369
370//----------------------------------------------------------------------
371// check if the jet has zero momentum
372bool operator==(const PseudoJet & jet, const double val) {
373 if (val != 0)
374 throw Error("comparing a PseudoJet with a non-zero constant (double) is not allowed.");
375 return (jet.px() == 0 && jet.py() == 0 &&
376 jet.pz() == 0 && jet.E() == 0);
377}
378
379
380
381//----------------------------------------------------------------------
382/// transform this jet (given in the rest frame of prest) into a jet
383/// in the lab frame
384//
385// NB: code adapted from that in herwig f77 (checked how it worked
386// long ago)
387PseudoJet & PseudoJet::boost(const PseudoJet & prest) {
388
389 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
390 return *this;
391
392 double m_local = prest.m();
393 assert(m_local != 0);
394
395 double pf4 = ( px()*prest.px() + py()*prest.py()
396 + pz()*prest.pz() + E()*prest.E() )/m_local;
397 double fn = (pf4 + E()) / (prest.E() + m_local);
398 _px += fn*prest.px();
399 _py += fn*prest.py();
400 _pz += fn*prest.pz();
401 _E = pf4;
402
403 _finish_init(); // we need to recalculate phi,rap,kt2
404 return *this;
405}
406
407
408//----------------------------------------------------------------------
409/// transform this jet (given in lab) into a jet in the rest
410/// frame of prest
411//
412// NB: code adapted from that in herwig f77 (checked how it worked
413// long ago)
414PseudoJet & PseudoJet::unboost(const PseudoJet & prest) {
415
416 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
417 return *this;
418
419 double m_local = prest.m();
420 assert(m_local != 0);
421
422 double pf4 = ( -px()*prest.px() - py()*prest.py()
423 - pz()*prest.pz() + E()*prest.E() )/m_local;
424 double fn = (pf4 + E()) / (prest.E() + m_local);
425 _px -= fn*prest.px();
426 _py -= fn*prest.py();
427 _pz -= fn*prest.pz();
428 _E = pf4;
429
430 _finish_init(); // we need to recalculate phi,rap,kt2
431 return *this;
432}
433
434
435//----------------------------------------------------------------------
436/// returns true if the momenta of the two input jets are identical
437bool have_same_momentum(const PseudoJet & jeta, const PseudoJet & jetb) {
438 return jeta.px() == jetb.px()
439 && jeta.py() == jetb.py()
440 && jeta.pz() == jetb.pz()
441 && jeta.E() == jetb.E();
442}
443
444//----------------------------------------------------------------------
445void PseudoJet::set_cached_rap_phi(double rap_in, double phi_in) {
446 _rap = rap_in; _phi = phi_in;
447 if (_phi >= twopi) _phi -= twopi;
448 if (_phi < 0) _phi += twopi;
449#ifdef FASTJET_HAVE_THREAD_SAFETY
450 _init_status = Init_Done;
451#endif
452}
453
454//----------------------------------------------------------------------
455void PseudoJet::reset_momentum_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in) {
456 assert(phi_in < 2*twopi && phi_in > -twopi);
457 double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in);
458 double exprap = exp(y_in);
459 double pminus = ptm/exprap;
460 double pplus = ptm*exprap;
461 double px_local = pt_in*cos(phi_in);
462 double py_local = pt_in*sin(phi_in);
463 reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus));
464 set_cached_rap_phi(y_in,phi_in);
465}
466
467//----------------------------------------------------------------------
468/// return a pseudojet with the given pt, y, phi and mass
469PseudoJet PtYPhiM(double pt, double y, double phi, double m) {
470 assert(phi < 2*twopi && phi > -twopi);
471 double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m);
472 double exprap = exp(y);
473 double pminus = ptm/exprap;
474 double pplus = ptm*exprap;
475 double px = pt*cos(phi);
476 double py = pt*sin(phi);
477 PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
478 mom.set_cached_rap_phi(y,phi);
479 return mom;
480 //return PseudoJet(pt*cos(phi), pt*sin(phi), ptm*sinh(y), ptm*cosh(y));
481}
482
483
484//----------------------------------------------------------------------
485// return kt-distance between this jet and another one
486double PseudoJet::kt_distance(const PseudoJet & other) const {
487 //double distance = min(this->kt2(), other.kt2());
488 double distance = min(_kt2, other._kt2);
489 double dphi = abs(phi() - other.phi());
490 if (dphi > pi) {dphi = twopi - dphi;}
491 double drap = rap() - other.rap();
492 distance = distance * (dphi*dphi + drap*drap);
493 return distance;
494}
495
496
497//----------------------------------------------------------------------
498// return squared cylinder (eta-phi) distance between this jet and another one
499double PseudoJet::plain_distance(const PseudoJet & other) const {
500 double dphi = abs(phi() - other.phi());
501 if (dphi > pi) {dphi = twopi - dphi;}
502 double drap = rap() - other.rap();
503 return (dphi*dphi + drap*drap);
504}
505
506//----------------------------------------------------------------------
507/// returns other.phi() - this.phi(), i.e. the phi distance to
508/// other, constrained to be in range -pi .. pi
509double PseudoJet::delta_phi_to(const PseudoJet & other) const {
510 double dphi = other.phi() - phi();
511 if (dphi > pi) dphi -= twopi;
512 if (dphi < -pi) dphi += twopi;
513 return dphi;
514}
515
516
517string PseudoJet::description() const{
518 // the "default" case of a PJ which does not belong to any cluster sequence
519 if (!_structure)
520 return "standard PseudoJet (with no associated clustering information)";
521
522 // for all the other cases, the description comes from the structure
523 return _structure->description();
524}
525
526
527
528//----------------------------------------------------------------------
529//
530// The following methods access the associated jet structure (if any)
531//
532//----------------------------------------------------------------------
533
534
535//----------------------------------------------------------------------
536// check whether this PseudoJet has an associated parent
537// ClusterSequence
538bool PseudoJet::has_associated_cluster_sequence() const{
539 return (_structure) && (_structure->has_associated_cluster_sequence());
540}
541
542//----------------------------------------------------------------------
543// get a (const) pointer to the associated ClusterSequence (NULL if
544// inexistent)
545const ClusterSequence* PseudoJet::associated_cluster_sequence() const{
546 if (! has_associated_cluster_sequence()) return NULL;
547
548 return _structure->associated_cluster_sequence();
549}
550
551
552//----------------------------------------------------------------------
553// check whether this PseudoJet has an associated parent
554// ClusterSequence that is still valid
555bool PseudoJet::has_valid_cluster_sequence() const{
556 return (_structure) && (_structure->has_valid_cluster_sequence());
557}
558
559//----------------------------------------------------------------------
560// If there is a valid cluster sequence associated with this jet,
561// returns a pointer to it; otherwise throws an Error.
562//
563// Open question: should these errors be upgraded to classes of their
564// own so that they can be caught? [Maybe, but later]
565const ClusterSequence * PseudoJet::validated_cs() const {
566 return validated_structure_ptr()->validated_cs();
567}
568
569
570//----------------------------------------------------------------------
571// set the associated structure
572void PseudoJet::set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure_in){
573//std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
574//std::shared_ptr: // if the jet currently belongs to a cs, we need to release it before any chenge
575//std::shared_ptr: _release_jet_from_cs();
576//std::shared_ptr: #endif //FASTJET_HAVE_THREAD_SAFETY
577 _structure = structure_in;
578}
579
580//----------------------------------------------------------------------
581// return true if there is some structure associated with this PseudoJet
582bool PseudoJet::has_structure() const{
583 return (bool) _structure; // cast to bool has to be made explicit
584}
585
586//----------------------------------------------------------------------
587// return a pointer to the structure (of type
588// PseudoJetStructureBase*) associated with this PseudoJet.
589//
590// return NULL if there is no associated structure
591const PseudoJetStructureBase* PseudoJet::structure_ptr() const {
592 //if (!_structure) return NULL;
593 return _structure.get();
594}
595
596//----------------------------------------------------------------------
597// return a non-const pointer to the structure (of type
598// PseudoJetStructureBase*) associated with this PseudoJet.
599//
600// return NULL if there is no associated structure
601//
602// Only use this if you know what you are doing. In any case,
603// prefer the 'structure_ptr()' (the const version) to this method,
604// unless you really need a write access to the PseudoJet's
605// underlying structure.
606PseudoJetStructureBase* PseudoJet::structure_non_const_ptr(){
607 //if (!_structure) return NULL;
608 return _structure.get();
609}
610
611//----------------------------------------------------------------------
612// return a pointer to the structure (of type
613// PseudoJetStructureBase*) associated with this PseudoJet.
614//
615// throw an error if there is no associated structure
616const PseudoJetStructureBase* PseudoJet::validated_structure_ptr() const {
617 if (!_structure)
618 throw Error("Trying to access the structure of a PseudoJet which has no associated structure");
619 return _structure.get();
620}
621
622//----------------------------------------------------------------------
623// return a reference to the shared pointer to the
624// PseudoJetStructureBase associated with this PseudoJet
625const SharedPtr<PseudoJetStructureBase> & PseudoJet::structure_shared_ptr() const {
626 return _structure;
627}
628
629
630//----------------------------------------------------------------------
631// check if it has been recombined with another PseudoJet in which
632// case, return its partner through the argument. Otherwise,
633// 'partner' is set to 0.
634//
635// false is also returned if this PseudoJet has no associated
636// ClusterSequence
637bool PseudoJet::has_partner(PseudoJet &partner) const{
638 return validated_structure_ptr()->has_partner(*this, partner);
639}
640
641//----------------------------------------------------------------------
642// check if it has been recombined with another PseudoJet in which
643// case, return its child through the argument. Otherwise, 'child'
644// is set to 0.
645//
646// false is also returned if this PseudoJet has no associated
647// ClusterSequence, with the child set to 0
648bool PseudoJet::has_child(PseudoJet &child) const{
649 return validated_structure_ptr()->has_child(*this, child);
650}
651
652//----------------------------------------------------------------------
653// check if it is the product of a recombination, in which case
654// return the 2 parents through the 'parent1' and 'parent2'
655// arguments. Otherwise, set these to 0.
656//
657// false is also returned if this PseudoJet has no parent
658// ClusterSequence
659bool PseudoJet::has_parents(PseudoJet &parent1, PseudoJet &parent2) const{
660 return validated_structure_ptr()->has_parents(*this, parent1, parent2);
661}
662
663//----------------------------------------------------------------------
664// check if the current PseudoJet contains the one passed as
665// argument
666//
667// false is also returned if this PseudoJet has no associated
668// ClusterSequence.
669bool PseudoJet::contains(const PseudoJet &constituent) const{
670 return validated_structure_ptr()->object_in_jet(constituent, *this);
671}
672
673//----------------------------------------------------------------------
674// check if the current PseudoJet is contained the one passed as
675// argument
676//
677// false is also returned if this PseudoJet has no associated
678// ClusterSequence
679bool PseudoJet::is_inside(const PseudoJet &jet) const{
680 return validated_structure_ptr()->object_in_jet(*this, jet);
681}
682
683
684//----------------------------------------------------------------------
685// returns true if the PseudoJet has constituents
686bool PseudoJet::has_constituents() const{
687 return (_structure) && (_structure->has_constituents());
688}
689
690//----------------------------------------------------------------------
691// retrieve the constituents.
692vector<PseudoJet> PseudoJet::constituents() const{
693 return validated_structure_ptr()->constituents(*this);
694}
695
696
697//----------------------------------------------------------------------
698// returns true if the PseudoJet has support for exclusive subjets
699bool PseudoJet::has_exclusive_subjets() const{
700 return (_structure) && (_structure->has_exclusive_subjets());
701}
702
703//----------------------------------------------------------------------
704// return a vector of all subjets of the current jet (in the sense
705// of the exclusive algorithm) that would be obtained when running
706// the algorithm with the given dcut.
707//
708// Time taken is O(m ln m), where m is the number of subjets that
709// are found. If m gets to be of order of the total number of
710// constituents in the jet, this could be substantially slower than
711// just getting that list of constituents.
712//
713// an Error is thrown if this PseudoJet has no currently valid
714// associated ClusterSequence
715std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double dcut) const {
716 return validated_structure_ptr()->exclusive_subjets(*this, dcut);
717}
718
719//----------------------------------------------------------------------
720// return the size of exclusive_subjets(...); still n ln n with same
721// coefficient, but marginally more efficient than manually taking
722// exclusive_subjets.size()
723//
724// an Error is thrown if this PseudoJet has no currently valid
725// associated ClusterSequence
726int PseudoJet::n_exclusive_subjets(const double dcut) const {
727 return validated_structure_ptr()->n_exclusive_subjets(*this, dcut);
728}
729
730//----------------------------------------------------------------------
731// return the list of subjets obtained by unclustering the supplied
732// jet down to n subjets (or all constituents if there are fewer
733// than n).
734//
735// requires n ln n time
736//
737// an Error is thrown if this PseudoJet has no currently valid
738// associated ClusterSequence
739std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (int nsub) const {
740 return validated_structure_ptr()->exclusive_subjets_up_to(*this, nsub);
741}
742
743//----------------------------------------------------------------------
744// Same as exclusive_subjets_up_to but throws an error if there are
745// fewer than nsub particles in the jet
746std::vector<PseudoJet> PseudoJet::exclusive_subjets (int nsub) const {
747 vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub);
748 if (int(subjets.size()) < nsub) {
749 ostringstream err;
750 err << "Requested " << nsub << " exclusive subjets, but there were only "
751 << subjets.size() << " particles in the jet";
752 throw Error(err.str());
753 }
754 return subjets;
755}
756
757//----------------------------------------------------------------------
758// return the dij that was present in the merging nsub+1 -> nsub
759// subjets inside this jet.
760//
761// an Error is thrown if this PseudoJet has no currently valid
762// associated ClusterSequence
763double PseudoJet::exclusive_subdmerge(int nsub) const {
764 return validated_structure_ptr()->exclusive_subdmerge(*this, nsub);
765}
766
767//----------------------------------------------------------------------
768// return the maximum dij that occurred in the whole event at the
769// stage that the nsub+1 -> nsub merge of subjets occurred inside
770// this jet.
771//
772// an Error is thrown if this PseudoJet has no currently valid
773// associated ClusterSequence
774double PseudoJet::exclusive_subdmerge_max(int nsub) const {
775 return validated_structure_ptr()->exclusive_subdmerge_max(*this, nsub);
776}
777
778
779// returns true if a jet has pieces
780//
781// By default a single particle or a jet coming from a
782// ClusterSequence have no pieces and this methos will return false.
783bool PseudoJet::has_pieces() const{
784 return ((_structure) && (_structure->has_pieces(*this)));
785}
786
787// retrieve the pieces that make up the jet.
788//
789// By default a jet does not have pieces.
790// If the underlying interface supports "pieces" retrieve the
791// pieces from there.
792std::vector<PseudoJet> PseudoJet::pieces() const{
793 return validated_structure_ptr()->pieces(*this);
794 // if (!has_pieces())
795 // throw Error("Trying to retrieve the pieces of a PseudoJet that has no support for pieces.");
796 //
797 // return _structure->pieces(*this);
798}
799
800
801//----------------------------------------------------------------------
802// the following ones require a computation of the area in the
803// associated ClusterSequence (See ClusterSequenceAreaBase for details)
804//----------------------------------------------------------------------
805
806#ifndef __FJCORE__
807
808//----------------------------------------------------------------------
809// if possible, return a valid ClusterSequenceAreaBase pointer; otherwise
810// throw an error
811const ClusterSequenceAreaBase * PseudoJet::validated_csab() const {
812 const ClusterSequenceAreaBase *csab = dynamic_cast<const ClusterSequenceAreaBase*>(validated_cs());
813 if (csab == NULL) throw Error("you requested jet-area related information, but the PseudoJet does not have associated area information.");
814 return csab;
815}
816
817
818//----------------------------------------------------------------------
819// check if it has a defined area
820bool PseudoJet::has_area() const{
821 //if (! has_associated_cluster_sequence()) return false;
822 if (! has_structure()) return false;
823 return (validated_structure_ptr()->has_area() != 0);
824}
825
826//----------------------------------------------------------------------
827// return the jet (scalar) area.
828// throw an Error if there is no support for area in the associated CS
829double PseudoJet::area() const{
830 return validated_structure_ptr()->area(*this);
831}
832
833//----------------------------------------------------------------------
834// return the error (uncertainty) associated with the determination
835// of the area of this jet.
836// throws an Error if there is no support for area in the associated CS
837double PseudoJet::area_error() const{
838 return validated_structure_ptr()->area_error(*this);
839}
840
841//----------------------------------------------------------------------
842// return the jet 4-vector area
843// throws an Error if there is no support for area in the associated CS
844PseudoJet PseudoJet::area_4vector() const{
845 return validated_structure_ptr()->area_4vector(*this);
846}
847
848//----------------------------------------------------------------------
849// true if this jet is made exclusively of ghosts
850// throws an Error if there is no support for area in the associated CS
851bool PseudoJet::is_pure_ghost() const{
852 return validated_structure_ptr()->is_pure_ghost(*this);
853}
854
855#endif // __FJCORE__
856
857//----------------------------------------------------------------------
858//
859// end of the methods accessing the information in the associated
860// Cluster Sequence
861//
862//----------------------------------------------------------------------
863
864//----------------------------------------------------------------------
865/// provide a meaningful error message for InexistentUserInfo
866PseudoJet::InexistentUserInfo::InexistentUserInfo() : Error("you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null")
867{}
868
869
870//----------------------------------------------------------------------
871// sort the indices so that values[indices[0..n-1]] is sorted
872// into increasing order
873void sort_indices(vector<int> & indices,
874 const vector<double> & values) {
875 IndexedSortHelper index_sort_helper(&values);
876 sort(indices.begin(), indices.end(), index_sort_helper);
877}
878
879
880//----------------------------------------------------------------------
881/// return a vector of jets sorted into decreasing kt2
882vector<PseudoJet> sorted_by_pt(const vector<PseudoJet> & jets) {
883 vector<double> minus_kt2(jets.size());
884 for (size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
885 return objects_sorted_by_values(jets, minus_kt2);
886}
887
888//----------------------------------------------------------------------
889/// return a vector of jets sorted into increasing rapidity
890vector<PseudoJet> sorted_by_rapidity(const vector<PseudoJet> & jets) {
891 vector<double> rapidities(jets.size());
892 for (size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
893 return objects_sorted_by_values(jets, rapidities);
894}
895
896//----------------------------------------------------------------------
897/// return a vector of jets sorted into decreasing energy
898vector<PseudoJet> sorted_by_E(const vector<PseudoJet> & jets) {
899 vector<double> energies(jets.size());
900 for (size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
901 return objects_sorted_by_values(jets, energies);
902}
903
904//----------------------------------------------------------------------
905/// return a vector of jets sorted into increasing pz
906vector<PseudoJet> sorted_by_pz(const vector<PseudoJet> & jets) {
907 vector<double> pz(jets.size());
908 for (size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();}
909 return objects_sorted_by_values(jets, pz);
910}
911
912
913
914//-------------------------------------------------------------------------------
915// helper functions to build a jet made of pieces
916//-------------------------------------------------------------------------------
917
918// build a "CompositeJet" from the vector of its pieces
919//
920// In this case, E-scheme recombination is assumed to compute the
921// total momentum
922PseudoJet join(const vector<PseudoJet> & pieces){
923 // compute the total momentum
924 //--------------------------------------------------
925 PseudoJet result; // automatically initialised to 0
926 for (unsigned int i=0; i<pieces.size(); i++)
927 result += pieces[i];
928
929 // attach a CompositeJetStructure to the result
930 //--------------------------------------------------
931 CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces);
932
933 result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
934
935 return result;
936}
937
938// build a "CompositeJet" from a single PseudoJet
939PseudoJet join(const PseudoJet & j1){
940 return join(vector<PseudoJet>(1,j1));
941}
942
943// build a "CompositeJet" from two PseudoJet
944PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
945 vector<PseudoJet> pieces;
946 pieces.reserve(2);
947 pieces.push_back(j1);
948 pieces.push_back(j2);
949 return join(pieces);
950}
951
952// build a "CompositeJet" from 3 PseudoJet
953PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3){
954 vector<PseudoJet> pieces;
955 pieces.reserve(3);
956 pieces.push_back(j1);
957 pieces.push_back(j2);
958 pieces.push_back(j3);
959 return join(pieces);
960}
961
962// build a "CompositeJet" from 4 PseudoJet
963PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4){
964 vector<PseudoJet> pieces;
965 pieces.reserve(4);
966 pieces.push_back(j1);
967 pieces.push_back(j2);
968 pieces.push_back(j3);
969 pieces.push_back(j4);
970 return join(pieces);
971}
972
973
974
975
976FASTJET_END_NAMESPACE
977
base class that sets interface for extensions of ClusterSequence that provide information about the a...
virtual double area(const PseudoJet &) const
return the area associated with the given jet; this base class returns 0.
deals with clustering
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:52
Contains any information related to the clustering that should be directly accessible to PseudoJet.
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition: PseudoJet.hh:68
PseudoJet()
default constructor, which as of FJ3.0 provides an object for which all operations are now valid and ...
Definition: PseudoJet.hh:82
void set_cached_rap_phi(double rap, double phi)
in some cases when setting a 4-momentum, the user/program knows what rapidity and azimuth are associa...
Definition: PseudoJet.cc:445
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:138
const UserInfoBase * user_info_ptr() const
retrieve a pointer to the (const) user information
Definition: PseudoJet.hh:496
virtual bool is_pure_ghost() const
true if this jet is made exclusively of ghosts.
Definition: PseudoJet.cc:851
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:123
const PseudoJetStructureBase * structure_ptr() const
return a pointer to the structure (of type PseudoJetStructureBase*) associated with this PseudoJet.
Definition: PseudoJet.cc:591
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
Definition: PseudoJet.hh:830
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP)
Definition: PseudoJet.hh:1060
int user_index() const
return the user_index,
Definition: PseudoJet.hh:389
virtual std::vector< PseudoJet > pieces() const
retrieve the pieces that make up the jet.
Definition: PseudoJet.cc:792
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
Definition: PseudoJet.cc:844
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
Definition: SharedPtr.hh:341
Selector operator*(const Selector &s1, const Selector &s2)
successive application of 2 selectors
Definition: Selector.cc:559
vector< PseudoJet > sorted_by_rapidity(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing rapidity
Definition: PseudoJet.cc:890
const double MaxRap
Used to protect against parton-level events where pt can be zero for some partons,...
Definition: PseudoJet.hh:53
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
Definition: PseudoJet.cc:898
bool have_same_momentum(const PseudoJet &jeta, const PseudoJet &jetb)
returns true if the momenta of the two input jets are identical
Definition: PseudoJet.cc:437
bool operator==(const PseudoJet &a, const PseudoJet &b)
returns true if the 4 momentum components of the two PseudoJets are identical and all the internal in...
Definition: PseudoJet.cc:356
std::vector< T > objects_sorted_by_values(const std::vector< T > &objects, const std::vector< double > &values)
given a vector of values with a one-to-one correspondence with the vector of objects,...
Definition: PseudoJet.hh:984
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:882
double theta(const PseudoJet &a, const PseudoJet &b)
returns the angle between a and b
Definition: PseudoJet.hh:946
PseudoJet PtYPhiM(double pt, double y, double phi, double m)
return a pseudojet with the given pt, y, phi and mass
Definition: PseudoJet.cc:469
vector< PseudoJet > sorted_by_pz(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing pz
Definition: PseudoJet.cc:906