Fork me on GitHub

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