FastJet 3.0.6
PseudoJet.cc
00001 //STARTHEADER
00002 // $Id: PseudoJet.cc 3083 2013-04-06 12:15:17Z cacciari $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 
00030 #include "fastjet/Error.hh"
00031 #include "fastjet/PseudoJet.hh"
00032 #include "fastjet/ClusterSequence.hh"
00033 #ifndef __FJCORE__
00034 #include "fastjet/ClusterSequenceAreaBase.hh"
00035 #endif  // __FJCORE__
00036 #include "fastjet/CompositeJetStructure.hh"
00037 #include<valarray>
00038 #include<iostream>
00039 #include<sstream>
00040 #include<cmath>
00041 #include<algorithm>
00042 #include <cstdarg>
00043 
00044 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00045 
00046 using namespace std;
00047 
00048 
00049 //----------------------------------------------------------------------
00050 // another constructor...
00051 PseudoJet::PseudoJet(const double px_in, const double py_in, const double pz_in, const double E_in) {
00052   
00053   _E  = E_in ;
00054   _px = px_in;
00055   _py = py_in;
00056   _pz = pz_in;
00057 
00058   this->_finish_init();
00059 
00060   // some default values for the history and user indices
00061   _reset_indices();
00062 
00063 }
00064 
00065 
00066 //----------------------------------------------------------------------
00067 /// do standard end of initialisation
00068 void PseudoJet::_finish_init () {
00069   _kt2 = this->px()*this->px() + this->py()*this->py();
00070   _phi = pseudojet_invalid_phi;
00071   // strictly speaking, _rap does not need initialising, because
00072   // it's never used as long as _phi == pseudojet_invalid_phi
00073   // (and gets set when _phi is requested). However ATLAS
00074   // 2013-03-28 complained that they sometimes have NaN's in
00075   // _rap and this interferes with some of their internal validation. 
00076   // So we initialise it; penalty is about 0.3ns per PseudoJet out of about
00077   // 10ns total initialisation time (on a intel Core i7 2.7GHz)
00078   _rap = pseudojet_invalid_rap;
00079 }
00080 
00081 //----------------------------------------------------------------------
00082 void PseudoJet::_set_rap_phi() const {
00083 
00084   if (_kt2 == 0.0) {
00085     _phi = 0.0; } 
00086   else {
00087     _phi = atan2(this->py(),this->px());
00088   }
00089   if (_phi < 0.0) {_phi += twopi;}
00090   if (_phi >= twopi) {_phi -= twopi;} // can happen if phi=-|eps<1e-15|?
00091   if (this->E() == abs(this->pz()) && _kt2 == 0) {
00092     // Point has infinite rapidity -- convert that into a very large
00093     // number, but in such a way that different 0-pt momenta will have
00094     // different rapidities (so as to lift the degeneracy between
00095     // them) [this can be relevant at parton-level]
00096     double MaxRapHere = MaxRap + abs(this->pz());
00097     if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;}
00098   } else {
00099     // get the rapidity in a way that's modestly insensitive to roundoff
00100     // error when things pz,E are large (actually the best we can do without
00101     // explicit knowledge of mass)
00102     double effective_m2 = max(0.0,m2()); // force non tachyonic mass
00103     double E_plus_pz    = _E + abs(_pz); // the safer of p+, p-
00104     // p+/p- = (p+ p-) / (p-)^2 = (kt^2+m^2)/(p-)^2
00105     _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
00106     if (_pz > 0) {_rap = - _rap;}
00107   }
00108 
00109 }
00110 
00111 
00112 //----------------------------------------------------------------------
00113 // return a valarray four-momentum
00114 valarray<double> PseudoJet::four_mom() const {
00115   valarray<double> mom(4);
00116   mom[0] = _px;
00117   mom[1] = _py;
00118   mom[2] = _pz;
00119   mom[3] = _E ;
00120   return mom;
00121 }
00122 
00123 //----------------------------------------------------------------------
00124 // Return the component corresponding to the specified index.
00125 // taken from CLHEP
00126 double PseudoJet::operator () (int i) const {
00127   switch(i) {
00128   case X:
00129     return px();
00130   case Y:
00131     return py();
00132   case Z:
00133     return pz();
00134   case T:
00135     return e();
00136   default:
00137     ostringstream err;
00138     err << "PseudoJet subscripting: bad index (" << i << ")";
00139     throw Error(err.str());
00140   }
00141   return 0.;
00142 }  
00143 
00144 //----------------------------------------------------------------------
00145 // return the pseudorapidity
00146 double PseudoJet::pseudorapidity() const {
00147   if (px() == 0.0 && py() ==0.0) return MaxRap;
00148   if (pz() == 0.0) return 0.0;
00149 
00150   double theta = atan(perp()/pz());
00151   if (theta < 0) theta += pi;
00152   return -log(tan(theta/2));
00153 }
00154 
00155 //----------------------------------------------------------------------
00156 // return "sum" of two pseudojets
00157 PseudoJet operator+ (const PseudoJet & jet1, const PseudoJet & jet2) {
00158   //return PseudoJet(jet1.four_mom()+jet2.four_mom());
00159   return PseudoJet(jet1.px()+jet2.px(),
00160                    jet1.py()+jet2.py(),
00161                    jet1.pz()+jet2.pz(),
00162                    jet1.E() +jet2.E()  );
00163 } 
00164 
00165 //----------------------------------------------------------------------
00166 // return difference of two pseudojets
00167 PseudoJet operator- (const PseudoJet & jet1, const PseudoJet & jet2) {
00168   //return PseudoJet(jet1.four_mom()-jet2.four_mom());
00169   return PseudoJet(jet1.px()-jet2.px(),
00170                    jet1.py()-jet2.py(),
00171                    jet1.pz()-jet2.pz(),
00172                    jet1.E() -jet2.E()  );
00173 } 
00174 
00175 //----------------------------------------------------------------------
00176 // return the product, coeff * jet
00177 PseudoJet operator* (double coeff, const PseudoJet & jet) {
00178   //return PseudoJet(coeff*jet.four_mom());
00179   // the following code is hopefully more efficient
00180   PseudoJet coeff_times_jet(jet);
00181   coeff_times_jet *= coeff;
00182   return coeff_times_jet;
00183 } 
00184 
00185 //----------------------------------------------------------------------
00186 // return the product, coeff * jet
00187 PseudoJet operator* (const PseudoJet & jet, double coeff) {
00188   return coeff*jet;
00189 } 
00190 
00191 //----------------------------------------------------------------------
00192 // return the ratio, jet / coeff
00193 PseudoJet operator/ (const PseudoJet & jet, double coeff) {
00194   return (1.0/coeff)*jet;
00195 } 
00196 
00197 //----------------------------------------------------------------------
00198 /// multiply the jet's momentum by the coefficient
00199 void PseudoJet::operator*=(double coeff) {
00200   _px *= coeff;
00201   _py *= coeff;
00202   _pz *= coeff;
00203   _E  *= coeff;
00204   _kt2*= coeff*coeff;
00205   // phi and rap are unchanged
00206 }
00207 
00208 //----------------------------------------------------------------------
00209 /// divide the jet's momentum by the coefficient
00210 void PseudoJet::operator/=(double coeff) {
00211   (*this) *= 1.0/coeff;
00212 }
00213 
00214 
00215 //----------------------------------------------------------------------
00216 /// add the other jet's momentum to this jet
00217 void PseudoJet::operator+=(const PseudoJet & other_jet) {
00218   _px += other_jet._px;
00219   _py += other_jet._py;
00220   _pz += other_jet._pz;
00221   _E  += other_jet._E ;
00222   _finish_init(); // we need to recalculate phi,rap,kt2
00223 }
00224 
00225 
00226 //----------------------------------------------------------------------
00227 /// subtract the other jet's momentum from this jet
00228 void PseudoJet::operator-=(const PseudoJet & other_jet) {
00229   _px -= other_jet._px;
00230   _py -= other_jet._py;
00231   _pz -= other_jet._pz;
00232   _E  -= other_jet._E ;
00233   _finish_init(); // we need to recalculate phi,rap,kt2
00234 }
00235 
00236 //----------------------------------------------------------------------
00237 bool operator==(const PseudoJet & a, const PseudoJet & b) {
00238   if (a.px() != b.px()) return false;
00239   if (a.py() != b.py()) return false;
00240   if (a.pz() != b.pz()) return false;
00241   if (a.E () != b.E ()) return false;
00242   
00243   if (a.user_index()    != b.user_index()) return false;
00244   if (a.cluster_hist_index() != b.cluster_hist_index()) return false;
00245   if (a.user_info_ptr() != b.user_info_ptr()) return false;
00246   if (a.structure_ptr() != b.structure_ptr()) return false;
00247 
00248   return true;
00249 }
00250 
00251 //----------------------------------------------------------------------
00252 // check if the jet has zero momentum
00253 bool operator==(const PseudoJet & jet, const double val) {
00254   if (val != 0) 
00255     throw Error("comparing a PseudoJet with a non-zero constant (double) is not allowed.");
00256   return (jet.px() == 0 && jet.py() == 0 && 
00257           jet.pz() == 0 && jet.E() == 0);
00258 }
00259 
00260 
00261 
00262 //----------------------------------------------------------------------
00263 /// transform this jet (given in lab) into a jet in the rest
00264 /// frame of prest
00265 //
00266 // NB: code adapted from that in herwig f77 (checked how it worked
00267 // long ago)
00268 PseudoJet & PseudoJet::boost(const PseudoJet & prest) {
00269   
00270   if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0) 
00271     return *this;
00272 
00273   double m_local = prest.m();
00274   assert(m_local != 0);
00275 
00276   double pf4  = (  px()*prest.px() + py()*prest.py()
00277                  + pz()*prest.pz() + E()*prest.E() )/m_local;
00278   double fn   = (pf4 + E()) / (prest.E() + m_local);
00279   _px +=  fn*prest.px();
00280   _py +=  fn*prest.py();
00281   _pz +=  fn*prest.pz();
00282   _E = pf4;
00283 
00284   _finish_init(); // we need to recalculate phi,rap,kt2
00285   return *this;
00286 }
00287 
00288 
00289 //----------------------------------------------------------------------
00290 /// transform this jet (given in the rest frame of prest) into a jet
00291 /// in the lab frame;
00292 //
00293 // NB: code adapted from that in herwig f77 (checked how it worked
00294 // long ago)
00295 PseudoJet & PseudoJet::unboost(const PseudoJet & prest) {
00296   
00297   if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0) 
00298     return *this;
00299 
00300   double m_local = prest.m();
00301   assert(m_local != 0);
00302 
00303   double pf4  = ( -px()*prest.px() - py()*prest.py()
00304                  - pz()*prest.pz() + E()*prest.E() )/m_local;
00305   double fn   = (pf4 + E()) / (prest.E() + m_local);
00306   _px -=  fn*prest.px();
00307   _py -=  fn*prest.py();
00308   _pz -=  fn*prest.pz();
00309   _E = pf4;
00310 
00311   _finish_init(); // we need to recalculate phi,rap,kt2
00312   return *this;
00313 }
00314 
00315 
00316 //----------------------------------------------------------------------
00317 /// returns true if the momenta of the two input jets are identical
00318 bool have_same_momentum(const PseudoJet & jeta, const PseudoJet & jetb) {
00319   return jeta.px() == jetb.px()
00320     &&   jeta.py() == jetb.py()
00321     &&   jeta.pz() == jetb.pz()
00322     &&   jeta.E()  == jetb.E();
00323 }
00324 
00325 //----------------------------------------------------------------------
00326 void PseudoJet::set_cached_rap_phi(double rap_in, double phi_in) {
00327   _rap = rap_in; _phi = phi_in;
00328   if (_phi >= twopi) _phi -= twopi;
00329   if (_phi < 0)      _phi += twopi;
00330 }
00331 
00332 //----------------------------------------------------------------------
00333 void PseudoJet::reset_momentum_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in) {
00334   assert(phi_in < 2*twopi && phi_in > -twopi);
00335   double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in);
00336   double exprap = exp(y_in);
00337   double pminus = ptm/exprap;
00338   double pplus  = ptm*exprap;
00339   double px_local = pt_in*cos(phi_in);
00340   double py_local = pt_in*sin(phi_in);
00341   reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus));
00342   set_cached_rap_phi(y_in,phi_in);
00343 }
00344 
00345 //----------------------------------------------------------------------
00346 /// return a pseudojet with the given pt, y, phi and mass
00347 PseudoJet PtYPhiM(double pt, double y, double phi, double m) {
00348   assert(phi < 2*twopi && phi > -twopi);
00349   double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m);
00350   double exprap = exp(y);
00351   double pminus = ptm/exprap;
00352   double pplus  = ptm*exprap;
00353   double px = pt*cos(phi);
00354   double py = pt*sin(phi);
00355   PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
00356   mom.set_cached_rap_phi(y,phi);
00357   return mom;
00358   //return PseudoJet(pt*cos(phi), pt*sin(phi), ptm*sinh(y), ptm*cosh(y));
00359 }
00360 
00361 
00362 //----------------------------------------------------------------------
00363 // return kt-distance between this jet and another one
00364 double PseudoJet::kt_distance(const PseudoJet & other) const {
00365   //double distance = min(this->kt2(), other.kt2());
00366   double distance = min(_kt2, other._kt2);
00367   double dphi = abs(phi() - other.phi());
00368   if (dphi > pi) {dphi = twopi - dphi;}
00369   double drap = rap() - other.rap();
00370   distance = distance * (dphi*dphi + drap*drap);
00371   return distance;
00372 }
00373 
00374 
00375 //----------------------------------------------------------------------
00376 // return squared cylinder (eta-phi) distance between this jet and another one
00377 double PseudoJet::plain_distance(const PseudoJet & other) const {
00378   double dphi = abs(phi() - other.phi());
00379   if (dphi > pi) {dphi = twopi - dphi;}
00380   double drap = rap() - other.rap();
00381   return (dphi*dphi + drap*drap);
00382 }
00383 
00384 //----------------------------------------------------------------------
00385 /// returns other.phi() - this.phi(), i.e. the phi distance to
00386 /// other, constrained to be in range -pi .. pi
00387 double PseudoJet::delta_phi_to(const PseudoJet & other) const {
00388   double dphi = other.phi() - phi();
00389   if (dphi >  pi) dphi -= twopi;
00390   if (dphi < -pi) dphi += twopi;
00391   return dphi;
00392 }
00393 
00394 
00395 string PseudoJet::description() const{
00396   // the "default" case of a PJ which does not belong to any cluster sequence
00397   if (!_structure())
00398     return "standard PseudoJet (with no associated clustering information)";
00399   
00400   // for all the other cases, the description comes from the structure
00401   return _structure()->description();
00402 }
00403 
00404 
00405 
00406 //----------------------------------------------------------------------
00407 //
00408 // The following methods access the associated jet structure (if any)
00409 //
00410 //----------------------------------------------------------------------
00411 
00412 
00413 //----------------------------------------------------------------------
00414 // check whether this PseudoJet has an associated parent
00415 // ClusterSequence
00416 bool PseudoJet::has_associated_cluster_sequence() const{
00417   return (_structure()) && (_structure->has_associated_cluster_sequence());
00418 }
00419 
00420 //----------------------------------------------------------------------
00421 // get a (const) pointer to the associated ClusterSequence (NULL if
00422 // inexistent)
00423 const ClusterSequence* PseudoJet::associated_cluster_sequence() const{
00424   if (! has_associated_cluster_sequence()) return NULL;
00425 
00426   return _structure->associated_cluster_sequence();
00427 }
00428 
00429 
00430 //----------------------------------------------------------------------
00431 // check whether this PseudoJet has an associated parent
00432 // ClusterSequence that is still valid
00433 bool PseudoJet::has_valid_cluster_sequence() const{
00434   return (_structure()) && (_structure->has_valid_cluster_sequence());
00435 }
00436 
00437 //----------------------------------------------------------------------
00438 // If there is a valid cluster sequence associated with this jet,
00439 // returns a pointer to it; otherwise throws an Error.
00440 //
00441 // Open question: should these errors be upgraded to classes of their
00442 // own so that they can be caught? [Maybe, but later]
00443 const ClusterSequence * PseudoJet::validated_cs() const {
00444   return validated_structure_ptr()->validated_cs();
00445 }
00446 
00447 
00448 //----------------------------------------------------------------------
00449 // set the associated structure
00450 void PseudoJet::set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure_in){
00451   _structure = structure_in;
00452 }
00453 
00454 //----------------------------------------------------------------------
00455 // return true if there is some strusture associated with this PseudoJet
00456 bool PseudoJet::has_structure() const{
00457   return _structure();
00458 }
00459 
00460 //----------------------------------------------------------------------
00461 // return a pointer to the structure (of type
00462 // PseudoJetStructureBase*) associated with this PseudoJet.
00463 //
00464 // return NULL if there is no associated structure
00465 const PseudoJetStructureBase* PseudoJet::structure_ptr() const {
00466   if (!_structure()) return NULL;
00467   return _structure();
00468 }
00469   
00470 //----------------------------------------------------------------------
00471 // return a non-const pointer to the structure (of type
00472 // PseudoJetStructureBase*) associated with this PseudoJet.
00473 //
00474 // return NULL if there is no associated structure
00475 //
00476 // Only use this if you know what you are doing. In any case,
00477 // prefer the 'structure_ptr()' (the const version) to this method,
00478 // unless you really need a write access to the PseudoJet's
00479 // underlying structure.
00480 PseudoJetStructureBase* PseudoJet::structure_non_const_ptr(){
00481   if (!_structure()) return NULL;
00482   return _structure();
00483 }
00484   
00485 //----------------------------------------------------------------------
00486 // return a pointer to the structure (of type
00487 // PseudoJetStructureBase*) associated with this PseudoJet.
00488 //
00489 // throw an error if there is no associated structure
00490 const PseudoJetStructureBase* PseudoJet::validated_structure_ptr() const {
00491   if (!_structure()) 
00492     throw Error("Trying to access the structure of a PseudoJet which has no associated structure");
00493   return _structure();
00494 }
00495   
00496 //----------------------------------------------------------------------
00497 // return a reference to the shared pointer to the
00498 // PseudoJetStructureBase associated with this PseudoJet
00499 const SharedPtr<PseudoJetStructureBase> & PseudoJet::structure_shared_ptr() const {
00500   return _structure;
00501 }
00502 
00503 
00504 //----------------------------------------------------------------------
00505 // check if it has been recombined with another PseudoJet in which
00506 // case, return its partner through the argument. Otherwise,
00507 // 'partner' is set to 0.
00508 //
00509 // false is also returned if this PseudoJet has no associated
00510 // ClusterSequence
00511 bool PseudoJet::has_partner(PseudoJet &partner) const{
00512   return validated_structure_ptr()->has_partner(*this, partner);
00513 }
00514 
00515 //----------------------------------------------------------------------
00516 // check if it has been recombined with another PseudoJet in which
00517 // case, return its child through the argument. Otherwise, 'child'
00518 // is set to 0.
00519 // 
00520 // false is also returned if this PseudoJet has no associated
00521 // ClusterSequence, with the child set to 0
00522 bool PseudoJet::has_child(PseudoJet &child) const{
00523   return validated_structure_ptr()->has_child(*this, child);
00524 }
00525 
00526 //----------------------------------------------------------------------
00527 // check if it is the product of a recombination, in which case
00528 // return the 2 parents through the 'parent1' and 'parent2'
00529 // arguments. Otherwise, set these to 0.
00530 //
00531 // false is also returned if this PseudoJet has no parent
00532 // ClusterSequence
00533 bool PseudoJet::has_parents(PseudoJet &parent1, PseudoJet &parent2) const{
00534   return validated_structure_ptr()->has_parents(*this, parent1, parent2);
00535 }
00536 
00537 //----------------------------------------------------------------------
00538 // check if the current PseudoJet contains the one passed as
00539 // argument
00540 //
00541 // false is also returned if this PseudoJet has no associated
00542 // ClusterSequence.
00543 bool PseudoJet::contains(const PseudoJet &constituent) const{
00544   return validated_structure_ptr()->object_in_jet(constituent, *this);
00545 }
00546 
00547 //----------------------------------------------------------------------
00548 // check if the current PseudoJet is contained the one passed as
00549 // argument
00550 //
00551 // false is also returned if this PseudoJet has no associated
00552 // ClusterSequence
00553 bool PseudoJet::is_inside(const PseudoJet &jet) const{
00554   return validated_structure_ptr()->object_in_jet(*this, jet);
00555 }
00556 
00557 
00558 //----------------------------------------------------------------------
00559 // returns true if the PseudoJet has constituents
00560 bool PseudoJet::has_constituents() const{
00561   return (_structure()) && (_structure->has_constituents());
00562 }
00563 
00564 //----------------------------------------------------------------------
00565 // retrieve the constituents.
00566 vector<PseudoJet> PseudoJet::constituents() const{
00567   return validated_structure_ptr()->constituents(*this);
00568 }
00569 
00570 
00571 //----------------------------------------------------------------------
00572 // returns true if the PseudoJet has support for exclusive subjets
00573 bool PseudoJet::has_exclusive_subjets() const{
00574   return (_structure()) && (_structure->has_exclusive_subjets());
00575 }
00576 
00577 //----------------------------------------------------------------------
00578 // return a vector of all subjets of the current jet (in the sense
00579 // of the exclusive algorithm) that would be obtained when running
00580 // the algorithm with the given dcut. 
00581 //
00582 // Time taken is O(m ln m), where m is the number of subjets that
00583 // are found. If m gets to be of order of the total number of
00584 // constituents in the jet, this could be substantially slower than
00585 // just getting that list of constituents.
00586 //
00587 // an Error is thrown if this PseudoJet has no currently valid
00588 // associated ClusterSequence
00589 std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double & dcut) const {
00590   return validated_structure_ptr()->exclusive_subjets(*this, dcut);
00591 }
00592 
00593 //----------------------------------------------------------------------
00594 // return the size of exclusive_subjets(...); still n ln n with same
00595 // coefficient, but marginally more efficient than manually taking
00596 // exclusive_subjets.size()
00597 //
00598 // an Error is thrown if this PseudoJet has no currently valid
00599 // associated ClusterSequence
00600 int PseudoJet::n_exclusive_subjets(const double & dcut) const {
00601   return validated_structure_ptr()->n_exclusive_subjets(*this, dcut);
00602 }
00603 
00604 //----------------------------------------------------------------------
00605 // return the list of subjets obtained by unclustering the supplied
00606 // jet down to n subjets (or all constituents if there are fewer
00607 // than n).
00608 //
00609 // requires n ln n time
00610 //
00611 // an Error is thrown if this PseudoJet has no currently valid
00612 // associated ClusterSequence
00613 std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (int nsub) const {
00614   return validated_structure_ptr()->exclusive_subjets_up_to(*this, nsub);
00615 }
00616 
00617 //----------------------------------------------------------------------
00618 // Same as exclusive_subjets_up_to but throws an error if there are
00619 // fewer than nsub particles in the jet
00620 std::vector<PseudoJet> PseudoJet::exclusive_subjets (int nsub) const {
00621   vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub);
00622   if (int(subjets.size()) < nsub) {
00623     ostringstream err;
00624     err << "Requested " << nsub << " exclusive subjets, but there were only " 
00625         << subjets.size() << " particles in the jet";
00626     throw Error(err.str());
00627   }
00628   return subjets;
00629 }
00630 
00631 //----------------------------------------------------------------------
00632 // return the dij that was present in the merging nsub+1 -> nsub 
00633 // subjets inside this jet.
00634 //
00635 // an Error is thrown if this PseudoJet has no currently valid
00636 // associated ClusterSequence
00637 double PseudoJet::exclusive_subdmerge(int nsub) const {
00638   return validated_structure_ptr()->exclusive_subdmerge(*this, nsub);
00639 }
00640 
00641 //----------------------------------------------------------------------
00642 // return the maximum dij that occurred in the whole event at the
00643 // stage that the nsub+1 -> nsub merge of subjets occurred inside 
00644 // this jet.
00645 //
00646 // an Error is thrown if this PseudoJet has no currently valid
00647 // associated ClusterSequence
00648 double PseudoJet::exclusive_subdmerge_max(int nsub) const {
00649   return validated_structure_ptr()->exclusive_subdmerge_max(*this, nsub);
00650 }
00651 
00652 
00653 // returns true if a jet has pieces
00654 //
00655 // By default a single particle or a jet coming from a
00656 // ClusterSequence have no pieces and this methos will return false.
00657 bool PseudoJet::has_pieces() const{
00658   return ((_structure()) && (_structure->has_pieces(*this)));
00659 }
00660 
00661 // retrieve the pieces that make up the jet. 
00662 //
00663 // By default a jet does not have pieces.
00664 // If the underlying interface supports "pieces" retrieve the
00665 // pieces from there.
00666 std::vector<PseudoJet> PseudoJet::pieces() const{
00667   return validated_structure_ptr()->pieces(*this);
00668   // if (!has_pieces())
00669   //   throw Error("Trying to retrieve the pieces of a PseudoJet that has no support for pieces.");
00670   //
00671   // return _structure->pieces(*this);
00672 }
00673 
00674 
00675 //----------------------------------------------------------------------
00676 // the following ones require a computation of the area in the
00677 // associated ClusterSequence (See ClusterSequenceAreaBase for details)
00678 //----------------------------------------------------------------------
00679 
00680 #ifndef __FJCORE__
00681 
00682 //----------------------------------------------------------------------
00683 // if possible, return a valid ClusterSequenceAreaBase pointer; otherwise
00684 // throw an error
00685 const ClusterSequenceAreaBase * PseudoJet::validated_csab() const {
00686   const ClusterSequenceAreaBase *csab = dynamic_cast<const ClusterSequenceAreaBase*>(validated_cs());
00687   if (csab == NULL) throw Error("you requested jet-area related information, but the PseudoJet does not have associated area information.");
00688   return csab;
00689 }
00690 
00691 
00692 //----------------------------------------------------------------------
00693 // check if it has a defined area
00694 bool PseudoJet::has_area() const{
00695   //if (! has_associated_cluster_sequence()) return false;
00696   if (! has_structure()) return false;
00697   return (validated_structure_ptr()->has_area() != 0);
00698 }
00699 
00700 //----------------------------------------------------------------------
00701 // return the jet (scalar) area.
00702 // throw an Error if there is no support for area in the associated CS
00703 double PseudoJet::area() const{
00704   return validated_structure_ptr()->area(*this);
00705 }
00706 
00707 //----------------------------------------------------------------------
00708 // return the error (uncertainty) associated with the determination
00709 // of the area of this jet.
00710 // throws an Error if there is no support for area in the associated CS
00711 double PseudoJet::area_error() const{
00712   return validated_structure_ptr()->area_error(*this);
00713 }
00714 
00715 //----------------------------------------------------------------------
00716 // return the jet 4-vector area
00717 // throws an Error if there is no support for area in the associated CS
00718 PseudoJet PseudoJet::area_4vector() const{
00719   return validated_structure_ptr()->area_4vector(*this);
00720 }
00721 
00722 //----------------------------------------------------------------------
00723 // true if this jet is made exclusively of ghosts
00724 // throws an Error if there is no support for area in the associated CS
00725 bool PseudoJet::is_pure_ghost() const{
00726   return validated_structure_ptr()->is_pure_ghost(*this);
00727 }
00728 
00729 #endif  // __FJCORE__
00730 
00731 //----------------------------------------------------------------------
00732 //
00733 // end of the methods accessing the information in the associated
00734 // Cluster Sequence
00735 //
00736 //----------------------------------------------------------------------
00737 
00738 //----------------------------------------------------------------------
00739 /// provide a meaningful error message for InexistentUserInfo
00740 PseudoJet::InexistentUserInfo::InexistentUserInfo() : Error("you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null")
00741 {}
00742 
00743 
00744 //----------------------------------------------------------------------
00745 // sort the indices so that values[indices[0..n-1]] is sorted
00746 // into increasing order 
00747 void sort_indices(vector<int> & indices, 
00748                          const vector<double> & values) {
00749   IndexedSortHelper index_sort_helper(&values);
00750   sort(indices.begin(), indices.end(), index_sort_helper);
00751 }
00752 
00753 
00754 
00755 //----------------------------------------------------------------------
00756 /// given a vector of values with a one-to-one correspondence with the
00757 /// vector of objects, sort objects into an order such that the
00758 /// associated values would be in increasing order
00759 template<class T> vector<T>  objects_sorted_by_values(
00760                        const vector<T> & objects, 
00761                        const vector<double> & values) {
00762 
00763   assert(objects.size() == values.size());
00764 
00765   // get a vector of indices
00766   vector<int> indices(values.size());
00767   for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
00768   
00769   // sort the indices
00770   sort_indices(indices, values);
00771   
00772   // copy the objects 
00773   vector<T> objects_sorted(objects.size());
00774   
00775   // place the objects in the correct order
00776   for (size_t i = 0; i < indices.size(); i++) {
00777     objects_sorted[i] = objects[indices[i]];
00778   }
00779 
00780   return objects_sorted;
00781 }
00782 
00783 //----------------------------------------------------------------------
00784 /// return a vector of jets sorted into decreasing kt2
00785 vector<PseudoJet> sorted_by_pt(const vector<PseudoJet> & jets) {
00786   vector<double> minus_kt2(jets.size());
00787   for (size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
00788   return objects_sorted_by_values(jets, minus_kt2);
00789 }
00790 
00791 //----------------------------------------------------------------------
00792 /// return a vector of jets sorted into increasing rapidity
00793 vector<PseudoJet> sorted_by_rapidity(const vector<PseudoJet> & jets) {
00794   vector<double> rapidities(jets.size());
00795   for (size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
00796   return objects_sorted_by_values(jets, rapidities);
00797 }
00798 
00799 //----------------------------------------------------------------------
00800 /// return a vector of jets sorted into decreasing energy
00801 vector<PseudoJet> sorted_by_E(const vector<PseudoJet> & jets) {
00802   vector<double> energies(jets.size());
00803   for (size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
00804   return objects_sorted_by_values(jets, energies);
00805 }
00806 
00807 //----------------------------------------------------------------------
00808 /// return a vector of jets sorted into increasing pz
00809 vector<PseudoJet> sorted_by_pz(const vector<PseudoJet> & jets) {
00810   vector<double> pz(jets.size());
00811   for (size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();}
00812   return objects_sorted_by_values(jets, pz);
00813 }
00814 
00815 
00816 
00817 //-------------------------------------------------------------------------------
00818 // helper functions to build a jet made of pieces
00819 //-------------------------------------------------------------------------------
00820 
00821 // build a "CompositeJet" from the vector of its pieces
00822 //
00823 // In this case, E-scheme recombination is assumed to compute the
00824 // total momentum
00825 PseudoJet join(const vector<PseudoJet> & pieces){
00826   // compute the total momentum
00827   //--------------------------------------------------
00828   PseudoJet result;  // automatically initialised to 0
00829   for (unsigned int i=0; i<pieces.size(); i++)
00830     result += pieces[i];
00831 
00832   // attach a CompositeJetStructure to the result
00833   //--------------------------------------------------
00834   CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces);
00835 
00836   result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
00837 
00838   return result;
00839 }
00840 
00841 // build a "CompositeJet" from a single PseudoJet
00842 PseudoJet join(const PseudoJet & j1){
00843   return join(vector<PseudoJet>(1,j1));
00844 }
00845 
00846 // build a "CompositeJet" from two PseudoJet
00847 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
00848   vector<PseudoJet> pieces;
00849   pieces.push_back(j1);
00850   pieces.push_back(j2);
00851   return join(pieces);
00852 }
00853 
00854 // build a "CompositeJet" from 3 PseudoJet
00855 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3){
00856   vector<PseudoJet> pieces;
00857   pieces.push_back(j1);
00858   pieces.push_back(j2);
00859   pieces.push_back(j3);
00860   return join(pieces);
00861 }
00862 
00863 // build a "CompositeJet" from 4 PseudoJet
00864 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4){
00865   vector<PseudoJet> pieces;
00866   pieces.push_back(j1);
00867   pieces.push_back(j2);
00868   pieces.push_back(j3);
00869   pieces.push_back(j4);
00870   return join(pieces);
00871 }
00872 
00873 
00874 
00875 
00876 FASTJET_END_NAMESPACE
00877 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends