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