FastJet  3.3.2
PseudoJet.cc
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 
46 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
47 
48 using namespace std;
49 
50 
51 //----------------------------------------------------------------------
52 // another constructor...
53 PseudoJet::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
70 void 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 //----------------------------------------------------------------------
84 void 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
116 valarray<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
128 double 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
148 double 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
159 PseudoJet 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
169 PseudoJet 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
179 PseudoJet 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
193 PseudoJet operator* (const PseudoJet & jet, double coeff) {
194  return coeff*jet;
195 }
196 
197 //----------------------------------------------------------------------
198 // return the ratio, jet / coeff
199 PseudoJet operator/ (const PseudoJet & jet, double coeff) {
200  return (1.0/coeff)*jet;
201 }
202 
203 //----------------------------------------------------------------------
204 /// multiply the jet's momentum by the coefficient
205 PseudoJet & 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
225 PseudoJet & 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
233 PseudoJet & 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
245 PseudoJet & 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 //----------------------------------------------------------------------
255 bool 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
271 bool 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)
286 PseudoJet & 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)
313 PseudoJet & 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
336 bool 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 //----------------------------------------------------------------------
344 void 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 //----------------------------------------------------------------------
351 void 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
365 PseudoJet 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
382 double 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
395 double 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
405 double 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 
413 string 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
434 bool 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)
441 const 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
451 bool 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]
461 const ClusterSequence * PseudoJet::validated_cs() const {
462  return validated_structure_ptr()->validated_cs();
463 }
464 
465 
466 //----------------------------------------------------------------------
467 // set the associated structure
468 void 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
474 bool 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
483 const 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.
498 PseudoJetStructureBase* 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
508 const 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
517 const 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
529 bool 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
540 bool 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
551 bool 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.
561 bool 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
571 bool 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
578 bool PseudoJet::has_constituents() const{
579  return (_structure) && (_structure->has_constituents());
580 }
581 
582 //----------------------------------------------------------------------
583 // retrieve the constituents.
584 vector<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
591 bool 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
607 std::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
618 int 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
631 std::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
638 std::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
655 double 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
666 double 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.
675 bool 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.
684 std::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
703 const 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
712 bool 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
721 double 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
729 double 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
736 PseudoJet 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
743 bool 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
758 PseudoJet::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
765 void 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
774 vector<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
782 vector<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
790 vector<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
798 vector<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
814 PseudoJet 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 
826 
827  return result;
828 }
829 
830 // build a "CompositeJet" from a single PseudoJet
831 PseudoJet join(const PseudoJet & j1){
832  return join(vector<PseudoJet>(1,j1));
833 }
834 
835 // build a "CompositeJet" from two PseudoJet
836 PseudoJet 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
845 PseudoJet 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
855 PseudoJet 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 
868 FASTJET_END_NAMESPACE
869 
PseudoJet PtYPhiM(double pt, double y, double phi, double m)
return a pseudojet with the given pt, y, phi and mass
Definition: PseudoJet.cc:365
double theta(const PseudoJet &a, const PseudoJet &b)
returns the angle between a and b
Definition: PseudoJet.hh:869
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:344
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:774
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, sort objects into an order such that the associated values would be in increasing order (but don&#39;t actually touch the values vector in the process).
Definition: PseudoJet.hh:907
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
Definition: PseudoJet.cc:736
deals with clustering
const double MaxRap
Used to protect against parton-level events where pt can be zero for some partons, giving rapidity=infinity.
Definition: PseudoJet.hh:52
const UserInfoBase * user_info_ptr() const
retrieve a pointer to the (const) user information
Definition: PseudoJet.hh:450
Selector operator*(const Selector &s1, const Selector &s2)
successive application of 2 selectors
Definition: Selector.cc:559
vector< PseudoJet > sorted_by_pz(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing pz
Definition: PseudoJet.cc:798
vector< PseudoJet > sorted_by_rapidity(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing rapidity
Definition: PseudoJet.cc:782
The structure for a jet made of pieces.
virtual double area(const PseudoJet &) const
return the area associated with the given jet; this base class returns 0.
Contains any information related to the clustering that should be directly accessible to PseudoJet...
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:255
virtual bool is_pure_ghost() const
true if this jet is made exclusively of ghosts.
Definition: PseudoJet.cc:743
const PseudoJetStructureBase * structure_ptr() const
return a pointer to the structure (of type PseudoJetStructureBase*) associated with this PseudoJet...
Definition: PseudoJet.cc:483
base class that sets interface for extensions of ClusterSequence that provide information about the a...
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
Definition: PseudoJet.cc:790
virtual std::vector< PseudoJet > pieces() const
retrieve the pieces that make up the jet.
Definition: PseudoJet.cc:684
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) ...
Definition: PseudoJet.hh:970
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
Definition: PseudoJet.hh:784
an implementation of C++0x shared pointers (or boost&#39;s)
Definition: SharedPtr.hh:121
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:110
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:125
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:336
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure)
set the associated structure
Definition: PseudoJet.cc:468
const double pseudojet_invalid_phi
default value for phi, meaning it (and rapidity) have yet to be calculated)
Definition: PseudoJet.hh:55
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
int user_index() const
return the user_index,
Definition: PseudoJet.hh:343