FastJet  3.1.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
PseudoJet.cc
1 //FJSTARTHEADER
2 // $Id: PseudoJet.cc 3652 2014-09-03 13:31:13Z salam $
3 //
4 // Copyright (c) 2005-2014, 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 void 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 }
221 
222 //----------------------------------------------------------------------
223 /// divide the jet's momentum by the coefficient
224 void PseudoJet::operator/=(double coeff) {
225  (*this) *= 1.0/coeff;
226 }
227 
228 
229 //----------------------------------------------------------------------
230 /// add the other jet's momentum to this jet
231 void PseudoJet::operator+=(const PseudoJet & other_jet) {
232  _px += other_jet._px;
233  _py += other_jet._py;
234  _pz += other_jet._pz;
235  _E += other_jet._E ;
236  _finish_init(); // we need to recalculate phi,rap,kt2
237 }
238 
239 
240 //----------------------------------------------------------------------
241 /// subtract the other jet's momentum from this jet
242 void PseudoJet::operator-=(const PseudoJet & other_jet) {
243  _px -= other_jet._px;
244  _py -= other_jet._py;
245  _pz -= other_jet._pz;
246  _E -= other_jet._E ;
247  _finish_init(); // we need to recalculate phi,rap,kt2
248 }
249 
250 //----------------------------------------------------------------------
251 bool operator==(const PseudoJet & a, const PseudoJet & b) {
252  if (a.px() != b.px()) return false;
253  if (a.py() != b.py()) return false;
254  if (a.pz() != b.pz()) return false;
255  if (a.E () != b.E ()) return false;
256 
257  if (a.user_index() != b.user_index()) return false;
258  if (a.cluster_hist_index() != b.cluster_hist_index()) return false;
259  if (a.user_info_ptr() != b.user_info_ptr()) return false;
260  if (a.structure_ptr() != b.structure_ptr()) return false;
261 
262  return true;
263 }
264 
265 //----------------------------------------------------------------------
266 // check if the jet has zero momentum
267 bool operator==(const PseudoJet & jet, const double val) {
268  if (val != 0)
269  throw Error("comparing a PseudoJet with a non-zero constant (double) is not allowed.");
270  return (jet.px() == 0 && jet.py() == 0 &&
271  jet.pz() == 0 && jet.E() == 0);
272 }
273 
274 
275 
276 //----------------------------------------------------------------------
277 /// transform this jet (given in the rest frame of prest) into a jet
278 /// in the lab frame
279 //
280 // NB: code adapted from that in herwig f77 (checked how it worked
281 // long ago)
282 PseudoJet & PseudoJet::boost(const PseudoJet & prest) {
283 
284  if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
285  return *this;
286 
287  double m_local = prest.m();
288  assert(m_local != 0);
289 
290  double pf4 = ( px()*prest.px() + py()*prest.py()
291  + pz()*prest.pz() + E()*prest.E() )/m_local;
292  double fn = (pf4 + E()) / (prest.E() + m_local);
293  _px += fn*prest.px();
294  _py += fn*prest.py();
295  _pz += fn*prest.pz();
296  _E = pf4;
297 
298  _finish_init(); // we need to recalculate phi,rap,kt2
299  return *this;
300 }
301 
302 
303 //----------------------------------------------------------------------
304 /// transform this jet (given in lab) into a jet in the rest
305 /// frame of prest
306 //
307 // NB: code adapted from that in herwig f77 (checked how it worked
308 // long ago)
309 PseudoJet & PseudoJet::unboost(const PseudoJet & prest) {
310 
311  if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
312  return *this;
313 
314  double m_local = prest.m();
315  assert(m_local != 0);
316 
317  double pf4 = ( -px()*prest.px() - py()*prest.py()
318  - pz()*prest.pz() + E()*prest.E() )/m_local;
319  double fn = (pf4 + E()) / (prest.E() + m_local);
320  _px -= fn*prest.px();
321  _py -= fn*prest.py();
322  _pz -= fn*prest.pz();
323  _E = pf4;
324 
325  _finish_init(); // we need to recalculate phi,rap,kt2
326  return *this;
327 }
328 
329 
330 //----------------------------------------------------------------------
331 /// returns true if the momenta of the two input jets are identical
332 bool have_same_momentum(const PseudoJet & jeta, const PseudoJet & jetb) {
333  return jeta.px() == jetb.px()
334  && jeta.py() == jetb.py()
335  && jeta.pz() == jetb.pz()
336  && jeta.E() == jetb.E();
337 }
338 
339 //----------------------------------------------------------------------
340 void PseudoJet::set_cached_rap_phi(double rap_in, double phi_in) {
341  _rap = rap_in; _phi = phi_in;
342  if (_phi >= twopi) _phi -= twopi;
343  if (_phi < 0) _phi += twopi;
344 }
345 
346 //----------------------------------------------------------------------
347 void PseudoJet::reset_momentum_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in) {
348  assert(phi_in < 2*twopi && phi_in > -twopi);
349  double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in);
350  double exprap = exp(y_in);
351  double pminus = ptm/exprap;
352  double pplus = ptm*exprap;
353  double px_local = pt_in*cos(phi_in);
354  double py_local = pt_in*sin(phi_in);
355  reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus));
356  set_cached_rap_phi(y_in,phi_in);
357 }
358 
359 //----------------------------------------------------------------------
360 /// return a pseudojet with the given pt, y, phi and mass
361 PseudoJet PtYPhiM(double pt, double y, double phi, double m) {
362  assert(phi < 2*twopi && phi > -twopi);
363  double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m);
364  double exprap = exp(y);
365  double pminus = ptm/exprap;
366  double pplus = ptm*exprap;
367  double px = pt*cos(phi);
368  double py = pt*sin(phi);
369  PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
370  mom.set_cached_rap_phi(y,phi);
371  return mom;
372  //return PseudoJet(pt*cos(phi), pt*sin(phi), ptm*sinh(y), ptm*cosh(y));
373 }
374 
375 
376 //----------------------------------------------------------------------
377 // return kt-distance between this jet and another one
378 double PseudoJet::kt_distance(const PseudoJet & other) const {
379  //double distance = min(this->kt2(), other.kt2());
380  double distance = min(_kt2, other._kt2);
381  double dphi = abs(phi() - other.phi());
382  if (dphi > pi) {dphi = twopi - dphi;}
383  double drap = rap() - other.rap();
384  distance = distance * (dphi*dphi + drap*drap);
385  return distance;
386 }
387 
388 
389 //----------------------------------------------------------------------
390 // return squared cylinder (eta-phi) distance between this jet and another one
391 double PseudoJet::plain_distance(const PseudoJet & other) const {
392  double dphi = abs(phi() - other.phi());
393  if (dphi > pi) {dphi = twopi - dphi;}
394  double drap = rap() - other.rap();
395  return (dphi*dphi + drap*drap);
396 }
397 
398 //----------------------------------------------------------------------
399 /// returns other.phi() - this.phi(), i.e. the phi distance to
400 /// other, constrained to be in range -pi .. pi
401 double PseudoJet::delta_phi_to(const PseudoJet & other) const {
402  double dphi = other.phi() - phi();
403  if (dphi > pi) dphi -= twopi;
404  if (dphi < -pi) dphi += twopi;
405  return dphi;
406 }
407 
408 
409 string PseudoJet::description() const{
410  // the "default" case of a PJ which does not belong to any cluster sequence
411  if (!_structure())
412  return "standard PseudoJet (with no associated clustering information)";
413 
414  // for all the other cases, the description comes from the structure
415  return _structure()->description();
416 }
417 
418 
419 
420 //----------------------------------------------------------------------
421 //
422 // The following methods access the associated jet structure (if any)
423 //
424 //----------------------------------------------------------------------
425 
426 
427 //----------------------------------------------------------------------
428 // check whether this PseudoJet has an associated parent
429 // ClusterSequence
430 bool PseudoJet::has_associated_cluster_sequence() const{
431  return (_structure()) && (_structure->has_associated_cluster_sequence());
432 }
433 
434 //----------------------------------------------------------------------
435 // get a (const) pointer to the associated ClusterSequence (NULL if
436 // inexistent)
437 const ClusterSequence* PseudoJet::associated_cluster_sequence() const{
438  if (! has_associated_cluster_sequence()) return NULL;
439 
440  return _structure->associated_cluster_sequence();
441 }
442 
443 
444 //----------------------------------------------------------------------
445 // check whether this PseudoJet has an associated parent
446 // ClusterSequence that is still valid
447 bool PseudoJet::has_valid_cluster_sequence() const{
448  return (_structure()) && (_structure->has_valid_cluster_sequence());
449 }
450 
451 //----------------------------------------------------------------------
452 // If there is a valid cluster sequence associated with this jet,
453 // returns a pointer to it; otherwise throws an Error.
454 //
455 // Open question: should these errors be upgraded to classes of their
456 // own so that they can be caught? [Maybe, but later]
457 const ClusterSequence * PseudoJet::validated_cs() const {
458  return validated_structure_ptr()->validated_cs();
459 }
460 
461 
462 //----------------------------------------------------------------------
463 // set the associated structure
464 void PseudoJet::set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure_in){
465  _structure = structure_in;
466 }
467 
468 //----------------------------------------------------------------------
469 // return true if there is some strusture associated with this PseudoJet
470 bool PseudoJet::has_structure() const{
471  return _structure();
472 }
473 
474 //----------------------------------------------------------------------
475 // return a pointer to the structure (of type
476 // PseudoJetStructureBase*) associated with this PseudoJet.
477 //
478 // return NULL if there is no associated structure
479 const PseudoJetStructureBase* PseudoJet::structure_ptr() const {
480  if (!_structure()) return NULL;
481  return _structure();
482 }
483 
484 //----------------------------------------------------------------------
485 // return a non-const pointer to the structure (of type
486 // PseudoJetStructureBase*) associated with this PseudoJet.
487 //
488 // return NULL if there is no associated structure
489 //
490 // Only use this if you know what you are doing. In any case,
491 // prefer the 'structure_ptr()' (the const version) to this method,
492 // unless you really need a write access to the PseudoJet's
493 // underlying structure.
494 PseudoJetStructureBase* PseudoJet::structure_non_const_ptr(){
495  if (!_structure()) return NULL;
496  return _structure();
497 }
498 
499 //----------------------------------------------------------------------
500 // return a pointer to the structure (of type
501 // PseudoJetStructureBase*) associated with this PseudoJet.
502 //
503 // throw an error if there is no associated structure
504 const PseudoJetStructureBase* PseudoJet::validated_structure_ptr() const {
505  if (!_structure())
506  throw Error("Trying to access the structure of a PseudoJet which has no associated structure");
507  return _structure();
508 }
509 
510 //----------------------------------------------------------------------
511 // return a reference to the shared pointer to the
512 // PseudoJetStructureBase associated with this PseudoJet
513 const SharedPtr<PseudoJetStructureBase> & PseudoJet::structure_shared_ptr() const {
514  return _structure;
515 }
516 
517 
518 //----------------------------------------------------------------------
519 // check if it has been recombined with another PseudoJet in which
520 // case, return its partner through the argument. Otherwise,
521 // 'partner' is set to 0.
522 //
523 // false is also returned if this PseudoJet has no associated
524 // ClusterSequence
525 bool PseudoJet::has_partner(PseudoJet &partner) const{
526  return validated_structure_ptr()->has_partner(*this, partner);
527 }
528 
529 //----------------------------------------------------------------------
530 // check if it has been recombined with another PseudoJet in which
531 // case, return its child through the argument. Otherwise, 'child'
532 // is set to 0.
533 //
534 // false is also returned if this PseudoJet has no associated
535 // ClusterSequence, with the child set to 0
536 bool PseudoJet::has_child(PseudoJet &child) const{
537  return validated_structure_ptr()->has_child(*this, child);
538 }
539 
540 //----------------------------------------------------------------------
541 // check if it is the product of a recombination, in which case
542 // return the 2 parents through the 'parent1' and 'parent2'
543 // arguments. Otherwise, set these to 0.
544 //
545 // false is also returned if this PseudoJet has no parent
546 // ClusterSequence
547 bool PseudoJet::has_parents(PseudoJet &parent1, PseudoJet &parent2) const{
548  return validated_structure_ptr()->has_parents(*this, parent1, parent2);
549 }
550 
551 //----------------------------------------------------------------------
552 // check if the current PseudoJet contains the one passed as
553 // argument
554 //
555 // false is also returned if this PseudoJet has no associated
556 // ClusterSequence.
557 bool PseudoJet::contains(const PseudoJet &constituent) const{
558  return validated_structure_ptr()->object_in_jet(constituent, *this);
559 }
560 
561 //----------------------------------------------------------------------
562 // check if the current PseudoJet is contained the one passed as
563 // argument
564 //
565 // false is also returned if this PseudoJet has no associated
566 // ClusterSequence
567 bool PseudoJet::is_inside(const PseudoJet &jet) const{
568  return validated_structure_ptr()->object_in_jet(*this, jet);
569 }
570 
571 
572 //----------------------------------------------------------------------
573 // returns true if the PseudoJet has constituents
574 bool PseudoJet::has_constituents() const{
575  return (_structure()) && (_structure->has_constituents());
576 }
577 
578 //----------------------------------------------------------------------
579 // retrieve the constituents.
580 vector<PseudoJet> PseudoJet::constituents() const{
581  return validated_structure_ptr()->constituents(*this);
582 }
583 
584 
585 //----------------------------------------------------------------------
586 // returns true if the PseudoJet has support for exclusive subjets
587 bool PseudoJet::has_exclusive_subjets() const{
588  return (_structure()) && (_structure->has_exclusive_subjets());
589 }
590 
591 //----------------------------------------------------------------------
592 // return a vector of all subjets of the current jet (in the sense
593 // of the exclusive algorithm) that would be obtained when running
594 // the algorithm with the given dcut.
595 //
596 // Time taken is O(m ln m), where m is the number of subjets that
597 // are found. If m gets to be of order of the total number of
598 // constituents in the jet, this could be substantially slower than
599 // just getting that list of constituents.
600 //
601 // an Error is thrown if this PseudoJet has no currently valid
602 // associated ClusterSequence
603 std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double dcut) const {
604  return validated_structure_ptr()->exclusive_subjets(*this, dcut);
605 }
606 
607 //----------------------------------------------------------------------
608 // return the size of exclusive_subjets(...); still n ln n with same
609 // coefficient, but marginally more efficient than manually taking
610 // exclusive_subjets.size()
611 //
612 // an Error is thrown if this PseudoJet has no currently valid
613 // associated ClusterSequence
614 int PseudoJet::n_exclusive_subjets(const double dcut) const {
615  return validated_structure_ptr()->n_exclusive_subjets(*this, dcut);
616 }
617 
618 //----------------------------------------------------------------------
619 // return the list of subjets obtained by unclustering the supplied
620 // jet down to n subjets (or all constituents if there are fewer
621 // than n).
622 //
623 // requires n ln n time
624 //
625 // an Error is thrown if this PseudoJet has no currently valid
626 // associated ClusterSequence
627 std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (int nsub) const {
628  return validated_structure_ptr()->exclusive_subjets_up_to(*this, nsub);
629 }
630 
631 //----------------------------------------------------------------------
632 // Same as exclusive_subjets_up_to but throws an error if there are
633 // fewer than nsub particles in the jet
634 std::vector<PseudoJet> PseudoJet::exclusive_subjets (int nsub) const {
635  vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub);
636  if (int(subjets.size()) < nsub) {
637  ostringstream err;
638  err << "Requested " << nsub << " exclusive subjets, but there were only "
639  << subjets.size() << " particles in the jet";
640  throw Error(err.str());
641  }
642  return subjets;
643 }
644 
645 //----------------------------------------------------------------------
646 // return the dij that was present in the merging nsub+1 -> nsub
647 // subjets inside this jet.
648 //
649 // an Error is thrown if this PseudoJet has no currently valid
650 // associated ClusterSequence
651 double PseudoJet::exclusive_subdmerge(int nsub) const {
652  return validated_structure_ptr()->exclusive_subdmerge(*this, nsub);
653 }
654 
655 //----------------------------------------------------------------------
656 // return the maximum dij that occurred in the whole event at the
657 // stage that the nsub+1 -> nsub merge of subjets occurred inside
658 // this jet.
659 //
660 // an Error is thrown if this PseudoJet has no currently valid
661 // associated ClusterSequence
662 double PseudoJet::exclusive_subdmerge_max(int nsub) const {
663  return validated_structure_ptr()->exclusive_subdmerge_max(*this, nsub);
664 }
665 
666 
667 // returns true if a jet has pieces
668 //
669 // By default a single particle or a jet coming from a
670 // ClusterSequence have no pieces and this methos will return false.
671 bool PseudoJet::has_pieces() const{
672  return ((_structure()) && (_structure->has_pieces(*this)));
673 }
674 
675 // retrieve the pieces that make up the jet.
676 //
677 // By default a jet does not have pieces.
678 // If the underlying interface supports "pieces" retrieve the
679 // pieces from there.
680 std::vector<PseudoJet> PseudoJet::pieces() const{
681  return validated_structure_ptr()->pieces(*this);
682  // if (!has_pieces())
683  // throw Error("Trying to retrieve the pieces of a PseudoJet that has no support for pieces.");
684  //
685  // return _structure->pieces(*this);
686 }
687 
688 
689 //----------------------------------------------------------------------
690 // the following ones require a computation of the area in the
691 // associated ClusterSequence (See ClusterSequenceAreaBase for details)
692 //----------------------------------------------------------------------
693 
694 #ifndef __FJCORE__
695 
696 //----------------------------------------------------------------------
697 // if possible, return a valid ClusterSequenceAreaBase pointer; otherwise
698 // throw an error
699 const ClusterSequenceAreaBase * PseudoJet::validated_csab() const {
700  const ClusterSequenceAreaBase *csab = dynamic_cast<const ClusterSequenceAreaBase*>(validated_cs());
701  if (csab == NULL) throw Error("you requested jet-area related information, but the PseudoJet does not have associated area information.");
702  return csab;
703 }
704 
705 
706 //----------------------------------------------------------------------
707 // check if it has a defined area
708 bool PseudoJet::has_area() const{
709  //if (! has_associated_cluster_sequence()) return false;
710  if (! has_structure()) return false;
711  return (validated_structure_ptr()->has_area() != 0);
712 }
713 
714 //----------------------------------------------------------------------
715 // return the jet (scalar) area.
716 // throw an Error if there is no support for area in the associated CS
717 double PseudoJet::area() const{
718  return validated_structure_ptr()->area(*this);
719 }
720 
721 //----------------------------------------------------------------------
722 // return the error (uncertainty) associated with the determination
723 // of the area of this jet.
724 // throws an Error if there is no support for area in the associated CS
725 double PseudoJet::area_error() const{
726  return validated_structure_ptr()->area_error(*this);
727 }
728 
729 //----------------------------------------------------------------------
730 // return the jet 4-vector area
731 // throws an Error if there is no support for area in the associated CS
732 PseudoJet PseudoJet::area_4vector() const{
733  return validated_structure_ptr()->area_4vector(*this);
734 }
735 
736 //----------------------------------------------------------------------
737 // true if this jet is made exclusively of ghosts
738 // throws an Error if there is no support for area in the associated CS
739 bool PseudoJet::is_pure_ghost() const{
740  return validated_structure_ptr()->is_pure_ghost(*this);
741 }
742 
743 #endif // __FJCORE__
744 
745 //----------------------------------------------------------------------
746 //
747 // end of the methods accessing the information in the associated
748 // Cluster Sequence
749 //
750 //----------------------------------------------------------------------
751 
752 //----------------------------------------------------------------------
753 /// provide a meaningful error message for InexistentUserInfo
754 PseudoJet::InexistentUserInfo::InexistentUserInfo() : Error("you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null")
755 {}
756 
757 
758 //----------------------------------------------------------------------
759 // sort the indices so that values[indices[0..n-1]] is sorted
760 // into increasing order
761 void sort_indices(vector<int> & indices,
762  const vector<double> & values) {
763  IndexedSortHelper index_sort_helper(&values);
764  sort(indices.begin(), indices.end(), index_sort_helper);
765 }
766 
767 
768 
769 //----------------------------------------------------------------------
770 /// given a vector of values with a one-to-one correspondence with the
771 /// vector of objects, sort objects into an order such that the
772 /// associated values would be in increasing order
773 template<class T> vector<T> objects_sorted_by_values(
774  const vector<T> & objects,
775  const vector<double> & values) {
776 
777  assert(objects.size() == values.size());
778 
779  // get a vector of indices
780  vector<int> indices(values.size());
781  for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
782 
783  // sort the indices
784  sort_indices(indices, values);
785 
786  // copy the objects
787  vector<T> objects_sorted(objects.size());
788 
789  // place the objects in the correct order
790  for (size_t i = 0; i < indices.size(); i++) {
791  objects_sorted[i] = objects[indices[i]];
792  }
793 
794  return objects_sorted;
795 }
796 
797 //----------------------------------------------------------------------
798 /// return a vector of jets sorted into decreasing kt2
799 vector<PseudoJet> sorted_by_pt(const vector<PseudoJet> & jets) {
800  vector<double> minus_kt2(jets.size());
801  for (size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
802  return objects_sorted_by_values(jets, minus_kt2);
803 }
804 
805 //----------------------------------------------------------------------
806 /// return a vector of jets sorted into increasing rapidity
807 vector<PseudoJet> sorted_by_rapidity(const vector<PseudoJet> & jets) {
808  vector<double> rapidities(jets.size());
809  for (size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
810  return objects_sorted_by_values(jets, rapidities);
811 }
812 
813 //----------------------------------------------------------------------
814 /// return a vector of jets sorted into decreasing energy
815 vector<PseudoJet> sorted_by_E(const vector<PseudoJet> & jets) {
816  vector<double> energies(jets.size());
817  for (size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
818  return objects_sorted_by_values(jets, energies);
819 }
820 
821 //----------------------------------------------------------------------
822 /// return a vector of jets sorted into increasing pz
823 vector<PseudoJet> sorted_by_pz(const vector<PseudoJet> & jets) {
824  vector<double> pz(jets.size());
825  for (size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();}
826  return objects_sorted_by_values(jets, pz);
827 }
828 
829 
830 
831 //-------------------------------------------------------------------------------
832 // helper functions to build a jet made of pieces
833 //-------------------------------------------------------------------------------
834 
835 // build a "CompositeJet" from the vector of its pieces
836 //
837 // In this case, E-scheme recombination is assumed to compute the
838 // total momentum
839 PseudoJet join(const vector<PseudoJet> & pieces){
840  // compute the total momentum
841  //--------------------------------------------------
842  PseudoJet result; // automatically initialised to 0
843  for (unsigned int i=0; i<pieces.size(); i++)
844  result += pieces[i];
845 
846  // attach a CompositeJetStructure to the result
847  //--------------------------------------------------
848  CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces);
849 
850  result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
851 
852  return result;
853 }
854 
855 // build a "CompositeJet" from a single PseudoJet
856 PseudoJet join(const PseudoJet & j1){
857  return join(vector<PseudoJet>(1,j1));
858 }
859 
860 // build a "CompositeJet" from two PseudoJet
861 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
862  vector<PseudoJet> pieces;
863  pieces.reserve(2);
864  pieces.push_back(j1);
865  pieces.push_back(j2);
866  return join(pieces);
867 }
868 
869 // build a "CompositeJet" from 3 PseudoJet
870 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3){
871  vector<PseudoJet> pieces;
872  pieces.reserve(3);
873  pieces.push_back(j1);
874  pieces.push_back(j2);
875  pieces.push_back(j3);
876  return join(pieces);
877 }
878 
879 // build a "CompositeJet" from 4 PseudoJet
880 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4){
881  vector<PseudoJet> pieces;
882  pieces.reserve(4);
883  pieces.push_back(j1);
884  pieces.push_back(j2);
885  pieces.push_back(j3);
886  pieces.push_back(j4);
887  return join(pieces);
888 }
889 
890 
891 
892 
893 FASTJET_END_NAMESPACE
894 
PseudoJet PtYPhiM(double pt, double y, double phi, double m)
return a pseudojet with the given pt, y, phi and mass
Definition: PseudoJet.cc:361
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:123
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:340
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:799
const UserInfoBase * user_info_ptr() const
retrieve a pointer to the (const) user information
Definition: PseudoJet.hh:438
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
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
Definition: PseudoJet.hh:770
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:823
vector< PseudoJet > sorted_by_rapidity(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing rapidity
Definition: PseudoJet.cc:807
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:680
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:773
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:251
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:815
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
an implementation of C++0x shared pointers (or boost's)
Definition: SharedPtr.hh:116
virtual bool is_pure_ghost() const
true if this jet is made exclusively of ghosts.
Definition: PseudoJet.cc:739
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:332
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:108
int user_index() const
return the user_index,
Definition: PseudoJet.hh:331
PseudoJet()
default constructor, which as of FJ3.0 provides an object for which all operations are now valid and ...
Definition: PseudoJet.hh:79
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) ...
Definition: PseudoJet.hh:920
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
const PseudoJetStructureBase * structure_ptr() const
return a pointer to the structure (of type PseudoJetStructureBase*) associated with this PseudoJet...
Definition: PseudoJet.cc:479
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
Definition: PseudoJet.cc:732