PseudoJet.cc
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
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 #include<algorithm>
00039
00040 FASTJET_BEGIN_NAMESPACE
00041
00042 using namespace std;
00043
00044
00045
00046
00047 PseudoJet::PseudoJet(const double px, const double py, const double pz, const double E) {
00048
00049 _E = E ;
00050 _px = px;
00051 _py = py;
00052 _pz = pz;
00053
00054 this->_finish_init();
00055
00056
00057 _reset_indices();
00058
00059 }
00060
00061
00062
00064 void PseudoJet::_finish_init () {
00065 _kt2 = this->px()*this->px() + this->py()*this->py();
00066 if (_kt2 == 0.0) {
00067 _phi = 0.0; }
00068 else {
00069 _phi = atan2(this->py(),this->px());
00070 }
00071 if (_phi < 0.0) {_phi += twopi;}
00072 if (_phi >= twopi) {_phi -= twopi;}
00073 if (this->E() == abs(this->pz()) && _kt2 == 0) {
00074
00075
00076
00077
00078 double MaxRapHere = MaxRap + abs(this->pz());
00079 if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;}
00080 } else {
00081
00082
00083
00084 double effective_m2 = max(0.0,m2());
00085 double E_plus_pz = _E + abs(_pz);
00086
00087 _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
00088 if (_pz > 0) {_rap = - _rap;}
00089 }
00090
00092
00093
00094
00095
00096
00097
00098
00099
00100 }
00101
00102
00103
00104
00105 valarray<double> PseudoJet::four_mom() const {
00106 valarray<double> mom(4);
00107 mom[0] = _px;
00108 mom[1] = _py;
00109 mom[2] = _pz;
00110 mom[3] = _E ;
00111 return mom;
00112 }
00113
00114
00115
00116
00117 double PseudoJet::operator () (int i) const {
00118 switch(i) {
00119 case X:
00120 return px();
00121 case Y:
00122 return py();
00123 case Z:
00124 return pz();
00125 case T:
00126 return e();
00127 default:
00128 ostringstream err;
00129 err << "PseudoJet subscripting: bad index (" << i << ")";
00130 throw Error(err.str());
00131 }
00132 return 0.;
00133 }
00134
00135
00136
00137 double PseudoJet::pseudorapidity() const {
00138 if (px() == 0.0 && py() ==0.0) return MaxRap;
00139 if (pz() == 0.0) return 0.0;
00140
00141 double theta = atan(perp()/pz());
00142 if (theta < 0) theta += pi;
00143 return -log(tan(theta/2));
00144 }
00145
00146
00147
00148 PseudoJet operator+ (const PseudoJet & jet1, const PseudoJet & jet2) {
00149
00150 return PseudoJet(jet1.px()+jet2.px(),
00151 jet1.py()+jet2.py(),
00152 jet1.pz()+jet2.pz(),
00153 jet1.E() +jet2.E() );
00154 }
00155
00156
00157
00158 PseudoJet operator- (const PseudoJet & jet1, const PseudoJet & jet2) {
00159
00160 return PseudoJet(jet1.px()-jet2.px(),
00161 jet1.py()-jet2.py(),
00162 jet1.pz()-jet2.pz(),
00163 jet1.E() -jet2.E() );
00164 }
00165
00166
00167
00168 PseudoJet operator* (double coeff, const PseudoJet & jet) {
00169
00170
00171 PseudoJet coeff_times_jet(jet);
00172 coeff_times_jet *= coeff;
00173 return coeff_times_jet;
00174 }
00175
00176
00177
00178 PseudoJet operator* (const PseudoJet & jet, double coeff) {
00179 return coeff*jet;
00180 }
00181
00182
00183
00184 PseudoJet operator/ (const PseudoJet & jet, double coeff) {
00185 return (1.0/coeff)*jet;
00186 }
00187
00188
00190 void PseudoJet::operator*=(double coeff) {
00191 _px *= coeff;
00192 _py *= coeff;
00193 _pz *= coeff;
00194 _E *= coeff;
00195 _kt2*= coeff*coeff;
00196
00197 }
00198
00199
00201 void PseudoJet::operator/=(double coeff) {
00202 (*this) *= 1.0/coeff;
00203 }
00204
00205
00206
00208 void PseudoJet::operator+=(const PseudoJet & other_jet) {
00209 _px += other_jet._px;
00210 _py += other_jet._py;
00211 _pz += other_jet._pz;
00212 _E += other_jet._E ;
00213 _finish_init();
00214 }
00215
00216
00217
00219 void PseudoJet::operator-=(const PseudoJet & other_jet) {
00220 _px -= other_jet._px;
00221 _py -= other_jet._py;
00222 _pz -= other_jet._pz;
00223 _E -= other_jet._E ;
00224 _finish_init();
00225 }
00226
00227
00228
00231
00232
00233
00234 PseudoJet & PseudoJet::boost(const PseudoJet & prest) {
00235
00236 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
00237 return *this;
00238
00239 double m = prest.m();
00240 assert(m != 0);
00241
00242 double pf4 = ( px()*prest.px() + py()*prest.py()
00243 + pz()*prest.pz() + E()*prest.E() )/m;
00244 double fn = (pf4 + E()) / (prest.E() + m);
00245 _px += fn*prest.px();
00246 _py += fn*prest.py();
00247 _pz += fn*prest.pz();
00248 _E = pf4;
00249
00250 _finish_init();
00251 return *this;
00252 }
00253
00254
00255
00258
00259
00260
00261 PseudoJet & PseudoJet::unboost(const PseudoJet & prest) {
00262
00263 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
00264 return *this;
00265
00266 double m = prest.m();
00267 assert(m != 0);
00268
00269 double pf4 = ( -px()*prest.px() - py()*prest.py()
00270 - pz()*prest.pz() + E()*prest.E() )/m;
00271 double fn = (pf4 + E()) / (prest.E() + m);
00272 _px -= fn*prest.px();
00273 _py -= fn*prest.py();
00274 _pz -= fn*prest.pz();
00275 _E = pf4;
00276
00277 _finish_init();
00278 return *this;
00279 }
00280
00281
00282
00284 bool have_same_momentum(const PseudoJet & jeta, const PseudoJet & jetb) {
00285 return jeta.px() == jetb.px()
00286 && jeta.py() == jetb.py()
00287 && jeta.pz() == jetb.pz()
00288 && jeta.E() == jetb.E();
00289 }
00290
00291
00292
00294 PseudoJet PtYPhiM(double pt, double y, double phi, double m) {
00295 double ptm = sqrt(pt*pt+m*m);
00296 return PseudoJet(pt*cos(phi), pt*sin(phi), ptm*sinh(y), ptm*cosh(y));
00297 }
00298
00299
00300
00301
00302 double PseudoJet::kt_distance(const PseudoJet & other) const {
00303
00304 double distance = min(_kt2, other._kt2);
00305 double dphi = abs(_phi - other._phi);
00306 if (dphi > pi) {dphi = twopi - dphi;}
00307 double drap = _rap - other._rap;
00308 distance = distance * (dphi*dphi + drap*drap);
00309 return distance;
00310 }
00311
00312
00313
00314
00315 double PseudoJet::plain_distance(const PseudoJet & other) const {
00316 double dphi = abs(_phi - other._phi);
00317 if (dphi > pi) {dphi = twopi - dphi;}
00318 double drap = _rap - other._rap;
00319 return (dphi*dphi + drap*drap);
00320 }
00321
00322
00325 double PseudoJet::delta_phi_to(const PseudoJet & other) const {
00326 double dphi = abs(other._phi - _phi);
00327 if (dphi > pi) dphi -= twopi;
00328 if (dphi < -pi) dphi += twopi;
00329 return dphi;
00330 }
00331
00332
00333
00334
00335 void sort_indices(vector<int> & indices,
00336 const vector<double> & values) {
00337 IndexedSortHelper index_sort_helper(&values);
00338 sort(indices.begin(), indices.end(), index_sort_helper);
00339 }
00340
00341
00345 template<class T> vector<T> objects_sorted_by_values(
00346 const vector<T> & objects,
00347 const vector<double> & values) {
00348
00349 assert(objects.size() == values.size());
00350
00351
00352 vector<int> indices(values.size());
00353 for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
00354
00355
00356 sort_indices(indices, values);
00357
00358
00359 vector<T> objects_sorted(objects.size());
00360
00361
00362 for (size_t i = 0; i < indices.size(); i++) {
00363 objects_sorted[i] = objects[indices[i]];
00364 }
00365
00366 return objects_sorted;
00367 }
00368
00369
00371 vector<PseudoJet> sorted_by_pt(const vector<PseudoJet> & jets) {
00372 vector<double> minus_kt2(jets.size());
00373 for (size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
00374 return objects_sorted_by_values(jets, minus_kt2);
00375 }
00376
00377
00379 vector<PseudoJet> sorted_by_rapidity(const vector<PseudoJet> & jets) {
00380 vector<double> rapidities(jets.size());
00381 for (size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
00382 return objects_sorted_by_values(jets, rapidities);
00383 }
00384
00385
00387 vector<PseudoJet> sorted_by_E(const vector<PseudoJet> & jets) {
00388 vector<double> energies(jets.size());
00389 for (size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
00390 return objects_sorted_by_values(jets, energies);
00391 }
00392
00393
00395 vector<PseudoJet> sorted_by_pz(const vector<PseudoJet> & jets) {
00396 vector<double> pz(jets.size());
00397 for (size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();}
00398 return objects_sorted_by_values(jets, pz);
00399 }
00400
00401
00402 FASTJET_END_NAMESPACE
00403