00001
00002
00003
00004 #ifndef Basics_H
00005 #define Basics_H
00006
00007 #include <iostream>
00008 #include <vector>
00009 #include <cmath>
00010 #include <ctime>
00011 #include "Pythia7/Config/Pythia7.h"
00012
00013 namespace Pythia7 {
00014
00020 namespace Shower {
00021
00022 using namespace std;
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00038 inline double sqrtpos(double x) {if (x < 0.) return 0.; else return sqrt(x);}
00039
00045 inline double sqrtsgn(double x) {if (x < 0.) return -sqrt(-x);
00046 else return sqrt(x);}
00047
00051 inline double lambdaK(double a, double b, double c) {
00052 return max(0., pow(a - b - c, 2) - 4. * b * c);
00053 }
00054
00059 inline double lambdaKRoot(double a, double b, double c) {
00060 return sqrt(lambdaK(a, b, c));
00061 }
00062
00069 inline double pAbsLambdaK(double a, double b, double c) {
00070 return 0.5*sqrt(lambdaK(a, b, c)/a);
00071 }
00072
00078 class Rndm {
00079
00080 public:
00081
00083 Rndm() {;}
00090 Rndm(long inseed) {init(inseed);}
00091
00095 static void init(long = 0);
00096
00100 static double flat() ;
00101
00105 static double exp() { return -log(flat()) ;}
00106
00110 static double xexp() { return -log(flat() * flat()) ;}
00111
00116 static double gauss() ;
00117
00122 static long pick(const vector<double>&) ;
00123
00124 private:
00125
00129 static bool initrndm;
00130
00135 static bool savgauss;
00136
00140 static long i97;
00141
00145 static long j97;
00146
00150 static double u[97];
00151
00155 static double c;
00156
00160 static double cd;
00161
00165 static double cm;
00166
00170 static double save;
00171 };
00172
00173
00174
00179 class Hist{
00180
00181 public:
00182
00184 static long NLINES;
00185
00187 static double TOLERANCE;
00188
00190 static double TINY;
00191
00195 Hist() {}
00196
00204 Hist(string titlein, long nbinin = 100, double xminin = 0.,
00205 double xmaxin = 1.) {
00206 book(titlein, nbinin, xminin, xmaxin);
00207 }
00208
00212 Hist(const Hist& h) {
00213 title = h.title; nbin = h.nbin; xmin = h.xmin; xmax = h.xmax;
00214 dx = h.dx; nfill = h.nfill; under = h.under; inside = h.inside;
00215 over = h.over; res = h.res;
00216 }
00217
00221 Hist(string titlein, const Hist& h) {
00222 title = titlein; nbin = h.nbin; xmin = h.xmin; xmax = h.xmax;
00223 dx = h.dx; nfill = h.nfill; under = h.under; inside = h.inside;
00224 over = h.over;res = h.res;
00225 }
00226
00230 Hist& operator=(const Hist& h) {
00231 if(this != &h) {
00232 nbin = h.nbin; xmin = h.xmin; xmax = h.xmax;
00233 dx = h.dx; nfill = h.nfill; under = h.under; inside = h.inside;
00234 over = h.over; res = h.res;
00235 return *this;
00236 }
00237 }
00238
00243 void book(string = " ", long = 100, double = 0., double = 1.);
00244
00248 void name(string = " ");
00249
00253 void null() ;
00254
00258 void fill(double, double = 1.);
00259
00264 void table(ostream& = cout) const ;
00265
00269 bool sameSize(const Hist&) const ;
00270
00272 Hist& operator+=(const Hist&) ;
00274 Hist& operator-=(const Hist&) ;
00276 Hist& operator*=(const Hist&) ;
00278 Hist& operator/=(const Hist&) ;
00280 Hist& operator+=(double) ;
00282 Hist& operator-=(double) ;
00284 Hist& operator*=(double) ;
00286 Hist& operator/=(double) ;
00287
00289 friend Hist operator+(double, const Hist&);
00291 friend Hist operator+(const Hist&, double);
00293 friend Hist operator+(const Hist&, const Hist&);
00295 friend Hist operator-(double, const Hist&);
00297 friend Hist operator-(const Hist&, double);
00299 friend Hist operator-(const Hist&, const Hist&);
00301 friend Hist operator*(double, const Hist&);
00303 friend Hist operator*(const Hist&, double);
00305 friend Hist operator*(const Hist&, const Hist&);
00307 friend Hist operator/(double, const Hist&);
00309 friend Hist operator/(const Hist&, double);
00311 friend Hist operator/(const Hist&, const Hist&);
00312
00314 friend ostream& operator<<(ostream&, const Hist&) ;
00315
00317 friend void table(const vector<Hist>&, ostream& = cout) ;
00318
00319 private:
00320
00322 string title;
00324 long nbin;
00326 long nfill;
00328 double xmin;
00330 double xmax;
00332 double dx;
00334 double under;
00336 double inside;
00338 double over;
00340 vector<double> res;
00341
00342 };
00343
00344
00345
00346
00347 class RotBstMatrix;
00348
00349
00350
00351
00358 class Vec4 {
00359
00360 public:
00364 Vec4(double xin = 0., double yin = 0., double zin = 0., double tin = 0.)
00365 : xx(xin), yy(yin), zz(zin), tt(tin) {}
00367 Vec4(const Vec4& v)
00368 : xx(v.xx), yy(v.yy), zz(v.zz), tt(v.tt) {}
00370 Vec4& operator=(const Vec4& v) {
00371 if (this != &v) {xx = v.xx; yy = v.yy; zz = v.zz; tt = v.tt; }
00372 return *this; }
00373
00374
00378 void p(double xin, double yin, double zin, double tin)
00379 {xx = xin; yy = yin; zz = zin; tt = tin;}
00381 void px(double xin) {xx = xin;}
00383 void py(double yin) {yy = yin;}
00385 void pz(double zin) {zz = zin;}
00387 void e(double tin) {tt = tin;}
00388
00389
00393 double px() const {return xx;}
00395 double py() const {return yy;}
00397 double pz() const {return zz;}
00399 double e() const {return tt;}
00401 double m2calc() const {return tt*tt - xx*xx - yy*yy - zz*zz;}
00403 double mcalc() const {return sqrtsgn(tt*tt - xx*xx - yy*yy - zz*zz);}
00405 double pT() const {return sqrt(xx*xx + yy*yy);}
00407 double pT2() const {return xx*xx + yy*yy;}
00409 double pAbs() const {return sqrt(xx*xx + yy*yy + zz*zz);}
00411 double p2() const {return xx*xx + yy*yy + zz*zz;}
00413 double theta() const {return atan2(sqrt(xx*xx + yy*yy), zz);}
00415 double phi() const {return atan2(yy,xx);}
00416
00417
00421 void rescalep(double fac) {xx = fac * xx; yy = fac * yy; zz = fac * zz;}
00423 void flip(){xx = -xx; yy = -yy; zz = -zz;}
00425 void rot(double, double);
00427 void rotaxis(double, double, double, double);
00429 void rotaxis(double, const Vec4&);
00433 void bst(double, double, double);
00435 void bst(const Vec4&);
00437 void rotbst(const RotBstMatrix&);
00438
00439
00443 Vec4& operator+=(const Vec4& v) {xx += v.xx; yy += v.yy; zz += v.zz;
00444 tt += v.tt; return *this;}
00446 Vec4& operator-=(const Vec4& v) {xx -= v.xx; yy -= v.yy; zz -= v.zz;
00447 tt -= v.tt; return *this;}
00449 Vec4& operator*=(double f) {xx *= f; yy *= f; zz *= f;
00450 tt *= f; return *this;}
00452 Vec4& operator/=(double f) {xx /= f; yy /= f; zz /= f;
00453 tt /= f; return *this;}
00454
00455
00459 friend Vec4 operator+(const Vec4&, const Vec4&);
00463 friend Vec4 operator-(const Vec4&, const Vec4&);
00467 friend Vec4 operator*(double, const Vec4&);
00471 friend Vec4 operator*(const Vec4&, double);
00475 friend Vec4 operator/(const Vec4&, double);
00479 friend double operator*(const Vec4&, const Vec4&);
00480
00484 friend double dot3(const Vec4&, const Vec4&);
00485
00489 friend Vec4 cross3(const Vec4&, const Vec4&);
00490
00494 friend double costheta(const Vec4 & v1, const Vec4 & v2);
00495
00499 friend double theta(const Vec4& v1, const Vec4& v2) {
00500 return acos(costheta(v1, v2)); }
00501
00506 friend double cosphi(const Vec4 & v1, const Vec4 & v2, const Vec4 & n);
00507
00511 friend double phi(const Vec4 & v1, const Vec4 & v2, const Vec4 & n) {
00512 return acos(cosphi(v1, v2, n)); }
00513
00514
00518 friend ostream& operator<<(ostream&, const Vec4&) ;
00519
00520 private:
00521
00523 double xx;
00524
00526 double yy;
00527
00529 double zz;
00530
00532 double tt;
00533
00534 };
00535
00536
00537
00539 inline Vec4 operator+(const Vec4& v1, const Vec4& v2)
00540 {Vec4 v = v1 ; return v += v2;}
00541
00543 inline Vec4 operator-(const Vec4& v1, const Vec4& v2)
00544 {Vec4 v = v1 ; return v -= v2;}
00545
00547 inline Vec4 operator*(double f, const Vec4& v1)
00548 {Vec4 v = v1; return v *= f;}
00549
00551 inline Vec4 operator*(const Vec4& v1, double f)
00552 {Vec4 v = v1; return v *= f;}
00553
00555 inline Vec4 operator/(const Vec4& v1, double f)
00556 {Vec4 v = v1; return v /= f;}
00557
00559 inline double operator*(const Vec4& v1, const Vec4& v2)
00560 {return v1.tt*v2.tt - v1.xx*v2.xx - v1.yy*v2.yy - v1.zz*v2.zz;}
00561
00562
00563
00564
00571 class RotBstMatrix {
00572
00573 public:
00577 RotBstMatrix() {for (long i = 0; i < 4; ++i) { for (long j = 0; j < 4; ++j) {
00578 M[i][j] = (i==j) ? 1. : 0.; } } }
00580 RotBstMatrix(const RotBstMatrix& Min) {
00581 for (long i = 0; i < 4; ++i) { for (long j = 0; j < 4; ++j) {
00582 M[i][j] = Min.M[i][j]; } } }
00584 RotBstMatrix& operator=(const RotBstMatrix& Min) {if (this != &Min) {
00585 for (long i = 0; i < 4; ++i) { for (long j = 0; j < 4; ++j) {
00586 M[i][j] = Min.M[i][j]; } }} return *this; }
00587
00588
00592 void rot(double = 0., double = 0.);
00596 void rot(const Vec4& p);
00600 void bst(double = 0., double = 0., double = 0.);
00604 void bst(const Vec4&);
00608 void bstback(const Vec4&);
00612 void bst(const Vec4&, const Vec4&);
00616 void rotbst(const RotBstMatrix&);
00620 void invert();
00624 void toCMframe(const Vec4&, const Vec4&);
00628 void fromCMframe(const Vec4&, const Vec4&);
00632 void reset();
00633
00637 friend ostream& operator<<(ostream&, const RotBstMatrix&) ;
00638
00642 friend class Vec4;
00643
00644 private:
00646 double M[4][4];
00647
00648 };
00649
00650
00651
00652
00658 class Particle {
00659
00660 public:
00664 Particle() {idp = 0; statusp = 0; mother1p = -1; mother2p = -1; prevp = -1;
00665 colp = 0; anticolp = 0; pp = Vec4(0.,0.,0.,0.); mp = 0.; scalep = 0.;}
00667 Particle(long idin, long statusin = 0, long mother1in = -1,
00668 long mother2in = -1, long colin = 0, long anticolin = 0,
00669 Vec4 pin = Vec4(0.,0.,0.,0.), double min = 0., double scalein = 0.) {
00670 idp = idin; statusp = statusin; mother1p = mother1in;
00671 mother2p = mother2in; prevp = -1; colp = colin; anticolp = anticolin;
00672 pp = pin; mp = min; scalep = scalein;}
00674 Particle(long idin, long statusin = 0, long mother1in = -1,
00675 long mother2in = -1, long colin = 0, long anticolin = 0,
00676 double pxin = 0., double pyin = 0., double pzin = 0.,
00677 double ein = 0., double min = 0., double scalein = 0.) {
00678 idp = idin; statusp = statusin; mother1p = mother1in;
00679 mother2p = mother2in; prevp = -1; colp = colin; anticolp = anticolin;
00680 pp = Vec4(pxin, pyin, pzin, ein); mp = min; scalep = scalein;}
00682 Particle(const Particle& pt) {
00683 idp = pt.idp; statusp = pt.statusp; mother1p = pt.mother1p;
00684 mother2p = pt.mother2p; prevp = pt.prevp;
00688 colp = pt.colp; anticolp = pt.anticolp;
00689 pp = pt.pp; mp = pt.mp; scalep = pt.scalep;}
00691 Particle& operator=(const Particle& pt) {if (this != &pt) {
00692 idp = pt.idp; statusp = pt.statusp; mother1p = pt.mother1p;
00693 mother2p = pt.mother2p; prevp = pt.prevp;
00694 colp = pt.colp; anticolp = pt.anticolp;
00695 pp = pt.pp; mp = pt.mp; scalep = pt.scalep;} return *this; }
00696
00697
00701 void id(long idin) {idp = idin;}
00703 void status(long statusin) {statusp = statusin;}
00705 void addstatus(long change) {statusp += change;}
00707 void mother1(long mother1in) {mother1p = mother1in;}
00709 void mother2(long mother2in) {mother2p = mother2in;}
00711 void mothers(long mother1in = -1, long mother2in = -1)
00712 {mother1p = mother1in; mother2p = mother2in;}
00714 void prev(long previn) {prevp = previn;}
00716 void col(long colin) {colp = colin;}
00718 void anticol(long anticolin) {anticolp = anticolin;}
00720 void cols(long colin = 0,long anticolin = 0)
00721 {colp = colin; anticolp = anticolin;}
00723 void p(Vec4 pin) {pp = pin;}
00725 void p(double pxin, double pyin, double pzin, double ein)
00726 {pp.p(pxin, pyin, pzin, ein);}
00728 void px(double pxin) {pp.px(pxin);}
00730 void py(double pyin) {pp.py(pyin);}
00732 void pz(double pzin) {pp.pz(pzin);}
00734 void e(double ein) {pp.e(ein);}
00736 void m(double min) {mp = min;}
00738 void scale(double scalein) {scalep = scalein;}
00739
00740
00744 long id() const {return idp;}
00746 long status() const {return statusp;}
00748 long mother1() const {return mother1p;}
00750 long mother2() const {return mother2p;}
00752 long prev() const {return prevp;}
00754 long col() const {return colp;}
00756 long anticol() const {return anticolp;}
00758 Vec4 p() const {return pp;}
00760 double px() const {return pp.px();}
00762 double py() const {return pp.py();}
00764 double pz() const {return pp.pz();}
00766 double e() const {return pp.e();}
00768 double m() const {return mp;}
00770 double scale() const {return scalep;}
00772 double m2() const {return mp*mp;}
00774 double mcalc() const {return pp.mcalc();}
00776 double m2calc() const {return pp.m2calc();}
00778 double pT() const {return pp.pT();}
00780 double pT2() const {return pp.pT2();}
00782 double mT() const {return sqrt(mp*mp + pp.pT2());}
00784 double mT2() const {return mp*mp + pp.pT2();}
00786 double pAbs() const {return pp.pAbs();}
00788 double p2() const {return pp.p2();}
00790 double theta() const {return pp.theta();}
00792 double phi() const {return pp.phi();}
00793
00794
00798 void rescalep(double fac) {pp.rescalep(fac);}
00800 void rot(double theta, double phi) {pp.rot(theta, phi);}
00802 void bst(double betaX, double betaY, double betaZ)
00803 {pp.bst(betaX, betaY, betaZ);}
00805 void bst(const Vec4& vec) {pp.bst(vec);}
00807 void rotbst(const RotBstMatrix& M) {pp.rotbst(M);}
00808
00809
00813 friend ostream& operator<<(ostream&, const Particle&) ;
00814
00815 private:
00817 long idp;
00819 long statusp;
00821 long mother1p;
00823 long mother2p;
00825 long prevp;
00827 long colp;
00829 long anticolp;
00831 Vec4 pp;
00833 double mp;
00835 double scalep;
00836 };
00837
00838
00839
00841 double Mass(long);
00842
00844 long iSpin(long);
00845
00847 long iCharge(long);
00848
00850 long iColour(long);
00851
00852
00853 }
00854 }
00855
00856 #endif // Basics_H
