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