herwig is hosted by Hepforge, IPPP Durham
Herwig++  2.7.0
EventShapes.h
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 */