Herwig++
2.7.0
|
00001 // -*- C++ -*- 00002 // 00003 // EventShapes.h is a part of Herwig++ - A multi-purpose Monte Carlo 00004 // event generator Copyright (C) 2002-2011 The Herwig Collaboration 00005 // 00006 // Herwig++ is licenced under version 2 of the GPL, see COPYING for 00007 // details. Please respect the MCnet academic guidelines, see 00008 // GUIDELINES for details. 00009 // 00010 #ifndef HERWIG_EventShapes_H 00011 #define HERWIG_EventShapes_H 00012 // 00013 // This is the declaration of the EventShapes class. 00014 // 00015 00016 #include "ThePEG/Interface/Interfaced.h" 00017 #include "ThePEG/Handlers/AnalysisHandler.h" 00018 #include "ThePEG/Vectors/Lorentz5Vector.h" 00019 #include "ThePEG/Vectors/ThreeVector.h" 00020 #include "ThePEG/EventRecord/Particle.h" 00021 #include "EventShapes.fh" 00022 00023 namespace Herwig { 00024 00025 using namespace ThePEG; 00026 00036 class EventShapes: public Interfaced { 00037 00038 public: 00039 00043 EventShapes() : _thrustDone(false), _spherDone(false), _linTenDone(false), 00044 _hemDone(false), _useCmBoost(false), 00045 _mPlus(), _mMinus(), _bPlus(), _bMinus() 00046 {} 00047 00051 void reset(const tPVector &part) { 00052 _pv.resize(part.size()); 00053 for(unsigned int ix=0;ix<part.size();++ix) _pv[ix]=part[ix]->momentum(); 00054 _thrustDone = false; 00055 _spherDone = false; 00056 _linTenDone = false; 00057 _hemDone = false; 00058 _useCmBoost = false; 00059 } 00060 00061 00062 public: 00063 00071 double thrust() { 00072 checkThrust(); 00073 return _thrust[0]; 00074 } 00075 00079 double thrustMajor() { 00080 checkThrust(); 00081 return _thrust[1]; 00082 } 00083 00087 double thrustMinor() { 00088 checkThrust(); 00089 return _thrust[2]; 00090 } 00091 00095 double oblateness() { 00096 checkThrust(); 00097 return _thrust[1]-_thrust[2]; 00098 } 00099 00103 Axis thrustAxis() { 00104 checkThrust(); 00105 return _thrustAxis[0]; 00106 } 00107 00111 Axis majorAxis() { 00112 checkThrust(); 00113 return _thrustAxis[1]; 00114 } 00115 00119 Axis minorAxis() { 00120 checkThrust(); 00121 return _thrustAxis[2]; 00122 } 00124 00132 double CParameter() { 00133 checkLinTen(); 00134 return 3.*(_linTen[0]*_linTen[1]+_linTen[1]*_linTen[2] 00135 +_linTen[2]*_linTen[0]); 00136 } 00137 00141 double DParameter() { 00142 checkLinTen(); 00143 return 27.*(_linTen[0]*_linTen[1]*_linTen[2]); 00144 } 00145 00149 vector<double> linTenEigenValues() { 00150 checkLinTen(); 00151 return _linTen; 00152 } 00153 00154 00158 vector<Axis> linTenEigenVectors() { 00159 checkLinTen(); 00160 return _linTenAxis; 00161 } 00162 00164 00172 double sphericity() { 00173 checkSphericity(); 00174 return 3./2.*(_spher[1]+_spher[2]); 00175 } 00176 00180 double aplanarity() { 00181 checkSphericity(); 00182 return 3./2.*_spher[2]; 00183 } 00184 00185 00189 double planarity() { 00190 checkSphericity(); 00191 return _spher[1]-_spher[2]; 00192 } 00193 00197 Axis sphericityAxis() { 00198 checkSphericity(); 00199 return _spherAxis[0]; 00200 } 00201 00202 00206 vector<double> sphericityEigenValues() { 00207 checkSphericity(); 00208 return _spher; 00209 } 00210 00214 vector<Axis> sphericityEigenVectors() { 00215 checkSphericity(); 00216 return _spherAxis; 00217 } //@} 00218 00227 double Mhigh2() { 00228 checkHemispheres(); 00229 return _mPlus; 00230 } 00231 00236 double Mlow2() { 00237 checkHemispheres(); 00238 return _mMinus; 00239 } 00240 00245 double Mdiff2() { 00246 checkHemispheres(); 00247 return _mPlus-_mMinus; 00248 } 00249 00251 00259 double Bmax() { 00260 checkHemispheres(); 00261 return _bPlus; 00262 } 00263 00267 double Bmin() { 00268 checkHemispheres(); 00269 return _bMinus; 00270 } 00271 00275 double Bsum() { 00276 checkHemispheres(); 00277 return _bPlus+_bMinus; 00278 } 00279 00280 00284 double Bdiff() { 00285 checkHemispheres(); 00286 return _bPlus-_bMinus; 00287 } 00289 00294 00298 double getXi(const Lorentz5Momentum & p, 00299 const Energy & Ebeam) { 00300 return((Ebeam > 0*MeV && p.vect().mag() > 0*MeV) ? 00301 log(Ebeam/p.vect().mag()) : -1.); 00302 } 00303 00307 Energy getPt(const Lorentz5Momentum & p) { 00308 return p.perp(); 00309 } 00310 00314 double getRapidity(const Lorentz5Momentum & p) { 00315 return (p.t() > p.z() ? p.rapidity() : 1e99); 00316 } 00318 00326 Energy ptInT(const Lorentz5Momentum & p) { 00327 checkThrust(); 00328 return p.vect()*_thrustAxis[1]; 00329 } 00330 00335 Energy ptOutT(const Lorentz5Momentum & p) { 00336 checkThrust(); 00337 return p.vect()*_thrustAxis[2]; 00338 } 00339 00343 double yT(const Lorentz5Momentum & p) { 00344 checkThrust(); 00345 return (p.t() > p.vect()*_thrustAxis[0] ? 00346 p.rapidity(_thrustAxis[0]) : 1e99); 00347 } 00348 00353 Energy ptInS(const Lorentz5Momentum & p) { 00354 checkSphericity(); 00355 return p.vect()*_spherAxis[1]; 00356 } 00357 00362 Energy ptOutS(const Lorentz5Momentum & p) { 00363 checkSphericity(); 00364 return p.vect()*_spherAxis[2]; 00365 } 00366 00370 double yS(const Lorentz5Momentum & p) { 00371 checkSphericity(); 00372 return (p.t() > p.vect()*_spherAxis[0] ? 00373 p.rapidity(_spherAxis[0]) : 1e99); 00374 } 00376 00377 00386 void bookEEC(vector<double> & hi); 00387 00392 void normalizeEEC(vector<double> & hi, long evts) { 00393 for (unsigned int bin = 0; bin < hi.size(); bin++) bin /= (hi.size()*evts); 00394 } 00395 00400 double AEEC(vector<double> & hi, double& coschi) { 00401 if (coschi > 0. && coschi <= 1.) { 00402 int i = static_cast<int>( floor((-coschi+1.)/2.*hi.size()) ); 00403 int j = static_cast<int>( floor(( coschi+1.)/2.*hi.size()) ); 00404 return hi[i]-hi[j]; 00405 } else { 00406 return 1e99; 00407 } 00408 } 00409 00410 public: 00411 00418 static void Init(); 00419 00420 protected: 00421 00428 virtual IBPtr clone() const {return new_ptr(*this);} 00429 00433 virtual IBPtr fullclone() const {return new_ptr(*this);} 00435 00436 private: 00437 00447 void checkThrust() { 00448 if (!_thrustDone) { 00449 _thrustDone = true; 00450 calculateThrust(); 00451 } 00452 } 00453 00458 void checkLinTen() { 00459 if (!_linTenDone) { 00460 _linTenDone = true; 00461 diagonalizeTensors(true, _useCmBoost); 00462 } 00463 } 00464 00469 void checkSphericity() { 00470 if (!_spherDone) { 00471 _spherDone = true; 00472 diagonalizeTensors(false, _useCmBoost); 00473 } 00474 } 00475 00480 void checkHemispheres() { 00481 if (!_hemDone) { 00482 _hemDone = true; 00483 calcHemisphereMasses(); 00484 } 00485 } 00487 00495 void calcHemisphereMasses(); 00496 00500 void calculateThrust(); 00501 00508 void diagonalizeTensors(bool linear, bool cmboost); 00509 00517 vector<double> eigenvalues(const double T[3][3]); 00518 00522 Axis eigenvector(const double T[3][3], const double &lam); 00523 00529 vector<Axis> eigenvectors(const double T[3][3], const vector<double> &lam); 00530 00537 void calcT(const vector<Momentum3> &p, Energy2 &t, Axis &taxis); 00538 00545 void calcM(const vector<Momentum3> &p, Energy2 &m, Axis &maxis); 00547 00548 private: 00549 00554 static NoPIOClassDescription<EventShapes> initEventShapes; 00555 00560 EventShapes & operator=(const EventShapes &); 00561 00562 private: 00563 00567 vector<Lorentz5Momentum> _pv; 00568 00576 vector<Axis> _thrustAxis; 00577 00581 vector<Axis> _spherAxis; 00582 00586 vector<Axis> _linTenAxis; 00588 00596 vector<double> _thrust; 00597 00601 vector<double> _spher; 00602 00606 vector<double> _linTen; 00608 00616 bool _thrustDone; 00617 00621 bool _spherDone; 00622 00626 bool _linTenDone; 00627 00631 bool _hemDone; 00633 00637 bool _useCmBoost; 00638 00646 double _mPlus; 00647 00651 double _mMinus; 00653 00661 double _bPlus; 00662 00666 double _bMinus; 00668 }; 00669 00670 } 00671 00672 #include "ThePEG/Utilities/ClassTraits.h" 00673 00674 namespace ThePEG { 00675 00680 template <> 00681 struct BaseClassTrait<Herwig::EventShapes,1> { 00683 typedef Interfaced NthBase; 00684 }; 00685 00688 template <> 00689 struct ClassTraits<Herwig::EventShapes> 00690 : public ClassTraitsBase<Herwig::EventShapes> { 00692 static string className() { return "Herwig::EventShapes"; } 00696 static string library() { return "HwAnalysis.so"; } 00697 }; 00698 00701 } 00702 00703 #endif /* HERWIG_EventShapes_H */