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 set_cluster_hist_index(-1);
00057 set_user_index(-1);
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 PseudoJet operator+ (const PseudoJet & jet1, const PseudoJet & jet2) {
00138
00139 return PseudoJet(jet1.px()+jet2.px(),
00140 jet1.py()+jet2.py(),
00141 jet1.pz()+jet2.pz(),
00142 jet1.E() +jet2.E() );
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* (double coeff, const PseudoJet & jet) {
00158
00159
00160 PseudoJet coeff_times_jet(jet);
00161 coeff_times_jet *= coeff;
00162 return coeff_times_jet;
00163 }
00164
00165
00166
00167 PseudoJet operator* (const PseudoJet & jet, double coeff) {
00168 return coeff*jet;
00169 }
00170
00171
00172
00173 PseudoJet operator/ (const PseudoJet & jet, double coeff) {
00174 return (1.0/coeff)*jet;
00175 }
00176
00177
00179 void PseudoJet::operator*=(double coeff) {
00180 _px *= coeff;
00181 _py *= coeff;
00182 _pz *= coeff;
00183 _E *= coeff;
00184 _kt2*= coeff*coeff;
00185
00186 }
00187
00188
00190 void PseudoJet::operator/=(double coeff) {
00191 (*this) *= 1.0/coeff;
00192 }
00193
00194
00195
00197 void PseudoJet::operator+=(const PseudoJet & other_jet) {
00198 _px += other_jet._px;
00199 _py += other_jet._py;
00200 _pz += other_jet._pz;
00201 _E += other_jet._E ;
00202 _finish_init();
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
00219
00220
00221
00222 PseudoJet & PseudoJet::boost(const PseudoJet & prest) {
00223
00224 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
00225 return *this;
00226
00227 double m = prest.m();
00228 assert(m != 0);
00229
00230 double pf4 = ( px()*prest.px() + py()*prest.py()
00231 + pz()*prest.pz() + E()*prest.E() )/m;
00232 double fn = (pf4 + E()) / (prest.E() + m);
00233 _px += fn*prest.px();
00234 _py += fn*prest.py();
00235 _pz += fn*prest.pz();
00236 _E = pf4;
00237
00238 _finish_init();
00239 return *this;
00240 }
00241
00242
00243
00246
00247
00248
00249 PseudoJet & PseudoJet::unboost(const PseudoJet & prest) {
00250
00251 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
00252 return *this;
00253
00254 double m = prest.m();
00255 assert(m != 0);
00256
00257 double pf4 = ( -px()*prest.px() - py()*prest.py()
00258 - pz()*prest.pz() + E()*prest.E() )/m;
00259 double fn = (pf4 + E()) / (prest.E() + m);
00260 _px -= fn*prest.px();
00261 _py -= fn*prest.py();
00262 _pz -= fn*prest.pz();
00263 _E = pf4;
00264
00265 _finish_init();
00266 return *this;
00267 }
00268
00269
00270
00272 bool have_same_momentum(const PseudoJet & jeta, const PseudoJet & jetb) {
00273 return jeta.px() == jetb.px()
00274 && jeta.py() == jetb.py()
00275 && jeta.pz() == jetb.pz()
00276 && jeta.E() == jetb.E();
00277 }
00278
00279
00280
00281
00282 double PseudoJet::kt_distance(const PseudoJet & other) const {
00283
00284 double distance = min(_kt2, other._kt2);
00285 double dphi = abs(_phi - other._phi);
00286 if (dphi > pi) {dphi = twopi - dphi;}
00287 double drap = _rap - other._rap;
00288 distance = distance * (dphi*dphi + drap*drap);
00289 return distance;
00290 }
00291
00292
00293
00294
00295 double PseudoJet::plain_distance(const PseudoJet & other) const {
00296 double dphi = abs(_phi - other._phi);
00297 if (dphi > pi) {dphi = twopi - dphi;}
00298 double drap = _rap - other._rap;
00299 return (dphi*dphi + drap*drap);
00300 }
00301
00302
00303
00304
00305 void sort_indices(vector<int> & indices,
00306 const vector<double> & values) {
00307 IndexedSortHelper index_sort_helper(&values);
00308 sort(indices.begin(), indices.end(), index_sort_helper);
00309 }
00310
00311
00315 template<class T> vector<T> objects_sorted_by_values(
00316 const vector<T> & objects,
00317 const vector<double> & values) {
00318
00319 assert(objects.size() == values.size());
00320
00321
00322 vector<int> indices(values.size());
00323 for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
00324
00325
00326 sort_indices(indices, values);
00327
00328
00329 vector<T> objects_sorted(objects.size());
00330
00331
00332 for (size_t i = 0; i < indices.size(); i++) {
00333 objects_sorted[i] = objects[indices[i]];
00334 }
00335
00336 return objects_sorted;
00337 }
00338
00339
00341 vector<PseudoJet> sorted_by_pt(const vector<PseudoJet> & jets) {
00342 vector<double> minus_kt2(jets.size());
00343 for (size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
00344 return objects_sorted_by_values(jets, minus_kt2);
00345 }
00346
00347
00349 vector<PseudoJet> sorted_by_rapidity(const vector<PseudoJet> & jets) {
00350 vector<double> rapidities(jets.size());
00351 for (size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
00352 return objects_sorted_by_values(jets, rapidities);
00353 }
00354
00355
00357 vector<PseudoJet> sorted_by_E(const vector<PseudoJet> & jets) {
00358 vector<double> energies(jets.size());
00359 for (size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
00360 return objects_sorted_by_values(jets, energies);
00361 }
00362
00363
00364 FASTJET_END_NAMESPACE
00365