PseudoJet.cc

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: PseudoJet.cc 958 2007-11-28 17:43:48Z cacciari $
00003 //
00004 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
00026 //  Foundation, Inc.:
00027 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00028 //----------------------------------------------------------------------
00029 //ENDHEADER
00030 
00031 
00032 #include "fastjet/Error.hh"
00033 #include "fastjet/PseudoJet.hh"
00034 #include<valarray>
00035 #include<iostream>
00036 #include<sstream>
00037 #include<cmath>
00038 
00039 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00040 
00041 using namespace std;
00042 
00043 
00044 //----------------------------------------------------------------------
00045 // another constructor...
00046 PseudoJet::PseudoJet(const double px, const double py, const double pz, const double E) {
00047   
00048   _E  = E ;
00049   _px = px;
00050   _py = py;
00051   _pz = pz;
00052 
00053   this->_finish_init();
00054 
00055   // some default values for the history and user indices
00056   _reset_indices();
00057 
00058 }
00059 
00060 
00061 //----------------------------------------------------------------------
00063 void PseudoJet::_finish_init () {
00064   _kt2 = this->px()*this->px() + this->py()*this->py();
00065   if (_kt2 == 0.0) {
00066     _phi = 0.0; } 
00067   else {
00068     _phi = atan2(this->py(),this->px());
00069   }
00070   if (_phi < 0.0) {_phi += twopi;}
00071   if (_phi >= twopi) {_phi -= twopi;} // can happen if phi=-|eps<1e-15|?
00072   if (this->E() == abs(this->pz()) && _kt2 == 0) {
00073     // Point has infinite rapidity -- convert that into a very large
00074     // number, but in such a way that different 0-pt momenta will have
00075     // different rapidities (so as to lift the degeneracy between
00076     // them) [this can be relevant at parton-level]
00077     double MaxRapHere = MaxRap + abs(this->pz());
00078     if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;}
00079   } else {
00080     // get the rapidity in a way that's modestly insensitive to roundoff
00081     // error when things pz,E are large (actually the best we can do without
00082     // explicit knowledge of mass)
00083     double effective_m2 = max(0.0,m2()); // force non tachyonic mass
00084     double E_plus_pz    = _E + abs(_pz); // the safer of p+, p-
00085     // p+/p- = (p+ p-) / (p-)^2 = (kt^2+m^2)/(p-)^2
00086     _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
00087     if (_pz > 0) {_rap = - _rap;}
00088   }
00089 
00091   //if (this->E() != abs(this->pz())) {
00092   //  _rap = 0.5*log((this->E() + this->pz())/(this->E() - this->pz()));
00093   //    } else {
00094   //  // Overlapping points can give problems. Let's lift the degeneracy
00095   //  // in case of multiple 0-pT points (can be found at parton-level)
00096   //  double MaxRapHere = MaxRap + abs(this->pz());
00097   //  if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;}
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 //----------------------------------------------------------------------
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 //----------------------------------------------------------------------
00200 void PseudoJet::operator/=(double coeff) {
00201   (*this) *= 1.0/coeff;
00202 }
00203 
00204 
00205 //----------------------------------------------------------------------
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 //----------------------------------------------------------------------
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 //----------------------------------------------------------------------
00230 //
00231 // NB: code adapted from that in herwig f77 (checked how it worked
00232 // long ago)
00233 PseudoJet & PseudoJet::boost(const PseudoJet & prest) {
00234   
00235   if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0) 
00236     return *this;
00237 
00238   double m = prest.m();
00239   assert(m != 0);
00240 
00241   double pf4  = (  px()*prest.px() + py()*prest.py()
00242                  + pz()*prest.pz() + E()*prest.E() )/m;
00243   double fn   = (pf4 + E()) / (prest.E() + m);
00244   _px +=  fn*prest.px();
00245   _py +=  fn*prest.py();
00246   _pz +=  fn*prest.pz();
00247   _E = pf4;
00248 
00249   _finish_init(); // we need to recalculate phi,rap,kt2
00250   return *this;
00251 }
00252 
00253 
00254 //----------------------------------------------------------------------
00257 //
00258 // NB: code adapted from that in herwig f77 (checked how it worked
00259 // long ago)
00260 PseudoJet & PseudoJet::unboost(const PseudoJet & prest) {
00261   
00262   if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0) 
00263     return *this;
00264 
00265   double m = prest.m();
00266   assert(m != 0);
00267 
00268   double pf4  = ( -px()*prest.px() - py()*prest.py()
00269                  - pz()*prest.pz() + E()*prest.E() )/m;
00270   double fn   = (pf4 + E()) / (prest.E() + m);
00271   _px -=  fn*prest.px();
00272   _py -=  fn*prest.py();
00273   _pz -=  fn*prest.pz();
00274   _E = pf4;
00275 
00276   _finish_init(); // we need to recalculate phi,rap,kt2
00277   return *this;
00278 }
00279 
00280 
00281 //----------------------------------------------------------------------
00283 bool have_same_momentum(const PseudoJet & jeta, const PseudoJet & jetb) {
00284   return jeta.px() == jetb.px()
00285     &&   jeta.py() == jetb.py()
00286     &&   jeta.pz() == jetb.pz()
00287     &&   jeta.E()  == jetb.E();
00288 }
00289 
00290 
00291 //----------------------------------------------------------------------
00293 PseudoJet PtYPhiM(double pt, double y, double phi, double m) {
00294   double ptm = sqrt(pt*pt+m*m);
00295   return PseudoJet(pt*cos(phi), pt*sin(phi), ptm*sinh(y), ptm*cosh(y));
00296 }
00297 
00298 
00299 //----------------------------------------------------------------------
00300 // return kt-distance between this jet and another one
00301 double PseudoJet::kt_distance(const PseudoJet & other) const {
00302   //double distance = min(this->kt2(), other.kt2());
00303   double distance = min(_kt2, other._kt2);
00304   double dphi = abs(_phi - other._phi);
00305   if (dphi > pi) {dphi = twopi - dphi;}
00306   double drap = _rap - other._rap;
00307   distance = distance * (dphi*dphi + drap*drap);
00308   return distance;
00309 }
00310 
00311 
00312 //----------------------------------------------------------------------
00313 // return squared cylinder (eta-phi) distance between this jet and another one
00314 double PseudoJet::plain_distance(const PseudoJet & other) const {
00315   double dphi = abs(_phi - other._phi);
00316   if (dphi > pi) {dphi = twopi - dphi;}
00317   double drap = _rap - other._rap;
00318   return (dphi*dphi + drap*drap);
00319 }
00320 
00321 //----------------------------------------------------------------------
00322 // sort the indices so that values[indices[0..n-1]] is sorted
00323 // into increasing order 
00324 void sort_indices(vector<int> & indices, 
00325                          const vector<double> & values) {
00326   IndexedSortHelper index_sort_helper(&values);
00327   sort(indices.begin(), indices.end(), index_sort_helper);
00328 }
00329 
00330 //----------------------------------------------------------------------
00334 template<class T> vector<T>  objects_sorted_by_values(
00335                        const vector<T> & objects, 
00336                        const vector<double> & values) {
00337 
00338   assert(objects.size() == values.size());
00339 
00340   // get a vector of indices
00341   vector<int> indices(values.size());
00342   for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
00343   
00344   // sort the indices
00345   sort_indices(indices, values);
00346   
00347   // copy the objects 
00348   vector<T> objects_sorted(objects.size());
00349   
00350   // place the objects in the correct order
00351   for (size_t i = 0; i < indices.size(); i++) {
00352     objects_sorted[i] = objects[indices[i]];
00353   }
00354 
00355   return objects_sorted;
00356 }
00357 
00358 //----------------------------------------------------------------------
00360 vector<PseudoJet> sorted_by_pt(const vector<PseudoJet> & jets) {
00361   vector<double> minus_kt2(jets.size());
00362   for (size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
00363   return objects_sorted_by_values(jets, minus_kt2);
00364 }
00365 
00366 //----------------------------------------------------------------------
00368 vector<PseudoJet> sorted_by_rapidity(const vector<PseudoJet> & jets) {
00369   vector<double> rapidities(jets.size());
00370   for (size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
00371   return objects_sorted_by_values(jets, rapidities);
00372 }
00373 
00374 //----------------------------------------------------------------------
00376 vector<PseudoJet> sorted_by_E(const vector<PseudoJet> & jets) {
00377   vector<double> energies(jets.size());
00378   for (size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
00379   return objects_sorted_by_values(jets, energies);
00380 }
00381 
00382 
00383 FASTJET_END_NAMESPACE
00384 

Generated on Thu Jan 3 19:04:30 2008 for fastjet by  doxygen 1.5.2