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
00039 FASTJET_BEGIN_NAMESPACE
00040
00041 using namespace std;
00042
00043
00044
00045
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
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;}
00072 if (this->E() == abs(this->pz()) && _kt2 == 0) {
00073
00074
00075
00076
00077 double MaxRapHere = MaxRap + abs(this->pz());
00078 if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;}
00079 } else {
00080
00081
00082
00083 double effective_m2 = max(0.0,m2());
00084 double E_plus_pz = _E + abs(_pz);
00085
00086 _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
00087 if (_pz > 0) {_rap = - _rap;}
00088 }
00089
00091
00092
00093
00094
00095
00096
00097
00098
00099 }
00100
00101
00102
00103
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
00115
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
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
00147 PseudoJet operator+ (const PseudoJet & jet1, const PseudoJet & jet2) {
00148
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
00157 PseudoJet operator- (const PseudoJet & jet1, const PseudoJet & jet2) {
00158
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
00167 PseudoJet operator* (double coeff, const PseudoJet & jet) {
00168
00169
00170 PseudoJet coeff_times_jet(jet);
00171 coeff_times_jet *= coeff;
00172 return coeff_times_jet;
00173 }
00174
00175
00176
00177 PseudoJet operator* (const PseudoJet & jet, double coeff) {
00178 return coeff*jet;
00179 }
00180
00181
00182
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
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();
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();
00224 }
00225
00226
00227
00230
00231
00232
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();
00250 return *this;
00251 }
00252
00253
00254
00257
00258
00259
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();
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
00301 double PseudoJet::kt_distance(const PseudoJet & other) const {
00302
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
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
00323
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
00341 vector<int> indices(values.size());
00342 for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
00343
00344
00345 sort_indices(indices, values);
00346
00347
00348 vector<T> objects_sorted(objects.size());
00349
00350
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