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