FastJet 3.0.2
|
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