13 #ifndef Pythia8_Basics_H
14 #define Pythia8_Basics_H
16 #include "PythiaStdlib.h"
31 virtual double flat() = 0;
51 Rndm() : initRndm(
false), seedSave(0), sequence(0),
52 useExternalRndm(
false), rndmEngPtr(0) { }
53 Rndm(
int seedIn) : initRndm(
false), seedSave(0), sequence(0),
54 useExternalRndm(
false), rndmEngPtr(0) { init(seedIn);}
60 void init(
int seedIn = 0) ;
66 double exp() {
return -log(flat()) ;}
69 double xexp() {
return -log(flat() * flat()) ;}
72 double gauss() {
return sqrt(-2. * log(flat())) * cos(M_PI * flat());}
75 pair<double, double> gauss2() {
double r = sqrt(-2. * log(flat()));
76 double phi = 2. * M_PI * flat();
77 return pair<double, double>(r * sin(phi), r * cos(phi));}
80 int pick(
const vector<double>& prob) ;
83 bool dumpState(
string fileName);
84 bool readState(
string fileName);
89 static const int DEFAULTSEED;
93 int i97, j97, defaultSeed, seedSave;
95 double u[97], c, cd, cm;
119 Vec4(
double xIn = 0.,
double yIn = 0.,
double zIn = 0.,
double tIn = 0.)
120 : xx(xIn), yy(yIn), zz(zIn), tt(tIn) { }
121 Vec4(
const Vec4& v) : xx(v.xx), yy(v.yy), zz(v.zz), tt(v.tt) { }
122 Vec4& operator=(
const Vec4& v) {
if (
this != &v) { xx = v.xx; yy = v.yy;
123 zz = v.zz; tt = v.tt; }
return *
this; }
124 Vec4& operator=(
double value) { xx = value; yy = value; zz = value;
125 tt = value;
return *
this; }
128 void reset() {xx = 0.; yy = 0.; zz = 0.; tt = 0.;}
129 void p(
double xIn,
double yIn,
double zIn,
double tIn)
130 {xx = xIn; yy = yIn; zz = zIn; tt = tIn;}
131 void p(
Vec4 pIn) {xx = pIn.xx; yy = pIn.yy; zz = pIn.zz; tt = pIn.tt;}
132 void px(
double xIn) {xx = xIn;}
133 void py(
double yIn) {yy = yIn;}
134 void pz(
double zIn) {zz = zIn;}
135 void e(
double tIn) {tt = tIn;}
138 double px()
const {
return xx;}
139 double py()
const {
return yy;}
140 double pz()
const {
return zz;}
141 double e()
const {
return tt;}
142 double mCalc()
const {
double temp = tt*tt - xx*xx - yy*yy - zz*zz;
143 return (temp >= 0.) ? sqrt(temp) : -sqrt(-temp);}
144 double m2Calc()
const {
return tt*tt - xx*xx - yy*yy - zz*zz;}
145 double pT()
const {
return sqrt(xx*xx + yy*yy);}
146 double pT2()
const {
return xx*xx + yy*yy;}
147 double pAbs()
const {
return sqrt(xx*xx + yy*yy + zz*zz);}
148 double pAbs2()
const {
return xx*xx + yy*yy + zz*zz;}
149 double eT()
const {
double temp = xx*xx + yy*yy;
150 return tt * sqrt( temp / (temp + zz*zz) );}
151 double eT2()
const {
double temp = xx*xx + yy*yy;
152 return tt*tt * temp / (temp + zz*zz);}
153 double theta()
const {
return atan2(sqrt(xx*xx + yy*yy), zz);}
154 double phi()
const {
return atan2(yy,xx);}
155 double thetaXZ()
const {
return atan2(xx,zz);}
156 double pPos()
const {
return tt + zz;}
157 double pNeg()
const {
return tt - zz;}
160 void rescale3(
double fac) {xx *= fac; yy *= fac; zz *= fac;}
161 void rescale4(
double fac) {xx *= fac; yy *= fac; zz *= fac; tt *= fac;}
162 void flip3() {xx = -xx; yy = -yy; zz = -zz;}
163 void flip4() {xx = -xx; yy = -yy; zz = -zz; tt = -tt;}
164 void rot(
double thetaIn,
double phiIn);
165 void rotaxis(
double phiIn,
double nx,
double ny,
double nz);
166 void rotaxis(
double phiIn,
const Vec4& n);
167 void bst(
double betaX,
double betaY,
double betaZ);
168 void bst(
double betaX,
double betaY,
double betaZ,
double gamma);
169 void bst(
const Vec4& pIn);
170 void bst(
const Vec4& pIn,
double mIn);
171 void bstback(
const Vec4& pIn);
172 void bstback(
const Vec4& pIn,
double mIn);
176 Vec4 operator-() {
Vec4 tmp; tmp.xx = -xx; tmp.yy = -yy; tmp.zz = -zz;
177 tmp.tt = -tt;
return tmp;}
178 Vec4& operator+=(
const Vec4& v) {xx += v.xx; yy += v.yy; zz += v.zz;
179 tt += v.tt;
return *
this;}
180 Vec4& operator-=(
const Vec4& v) {xx -= v.xx; yy -= v.yy; zz -= v.zz;
181 tt -= v.tt;
return *
this;}
182 Vec4& operator*=(
double f) {xx *= f; yy *= f; zz *= f;
183 tt *= f;
return *
this;}
184 Vec4& operator/=(
double f) {xx /= f; yy /= f; zz /= f;
185 tt /= f;
return *
this;}
190 friend Vec4 operator*(
double f,
const Vec4& v1);
191 friend Vec4 operator*(
const Vec4& v1,
double f);
192 friend Vec4 operator/(
const Vec4& v1,
double f);
193 friend double operator*(
const Vec4& v1,
const Vec4& v2);
196 friend double m(
const Vec4& v1,
const Vec4& v2);
197 friend double m2(
const Vec4& v1,
const Vec4& v2);
200 friend double dot3(
const Vec4& v1,
const Vec4& v2);
204 friend double theta(
const Vec4& v1,
const Vec4& v2);
205 friend double costheta(
const Vec4& v1,
const Vec4& v2);
208 friend double phi(
const Vec4& v1,
const Vec4& v2);
209 friend double cosphi(
const Vec4& v1,
const Vec4& v2);
212 friend double phi(
const Vec4& v1,
const Vec4& v2,
const Vec4& n);
213 friend double cosphi(
const Vec4& v1,
const Vec4& v2,
const Vec4& n);
216 friend ostream& operator<<(ostream&,
const Vec4& v) ;
221 static const double TINY;
224 double xx, yy, zz, tt;
235 {
Vec4 v = v1 ;
return v += v2;}
237 inline Vec4 operator-(
const Vec4& v1,
const Vec4& v2)
238 {Vec4 v = v1 ;
return v -= v2;}
240 inline Vec4 operator*(
double f,
const Vec4& v1)
241 {Vec4 v = v1;
return v *= f;}
243 inline Vec4 operator*(
const Vec4& v1,
double f)
244 {Vec4 v = v1;
return v *= f;}
246 inline Vec4 operator/(
const Vec4& v1,
double f)
247 {Vec4 v = v1;
return v /= f;}
249 inline double operator*(
const Vec4& v1,
const Vec4& v2)
250 {
return v1.tt*v2.tt - v1.xx*v2.xx - v1.yy*v2.yy - v1.zz*v2.zz;}
253 double m(
const Vec4& v1,
const Vec4& v2);
254 double m2(
const Vec4& v1,
const Vec4& v2);
257 double dot3(
const Vec4& v1,
const Vec4& v2);
258 Vec4 cross3(
const Vec4& v1,
const Vec4& v2);
261 double theta(
const Vec4& v1,
const Vec4& v2);
262 double costheta(
const Vec4& v1,
const Vec4& v2);
265 double phi(
const Vec4& v1,
const Vec4& v2);
266 double cosphi(
const Vec4& v1,
const Vec4& v2);
269 double phi(
const Vec4& v1,
const Vec4& v2,
const Vec4& n);
270 double cosphi(
const Vec4& v1,
const Vec4& v2,
const Vec4& n);
273 ostream& operator<<(ostream&,
const Vec4& v) ;
286 RotBstMatrix() {
for (
int i = 0; i < 4; ++i) {
for (
int j = 0; j < 4; ++j)
287 { M[i][j] = (i==j) ? 1. : 0.; } } }
289 for (
int i = 0; i < 4; ++i) {
for (
int j = 0; j < 4; ++j) {
290 M[i][j] = Min.M[i][j]; } } }
292 for (
int i = 0; i < 4; ++i) {
for (
int j = 0; j < 4; ++j) {
293 M[i][j] = Min.M[i][j]; } } }
return *
this; }
296 void rot(
double = 0.,
double = 0.);
297 void rot(
const Vec4& p);
298 void bst(
double = 0.,
double = 0.,
double = 0.);
299 void bst(
const Vec4&);
300 void bstback(
const Vec4&);
302 void toCMframe(
const Vec4&,
const Vec4&);
303 void fromCMframe(
const Vec4&,
const Vec4&);
309 double deviation()
const;
312 friend ostream& operator<<(ostream&,
const RotBstMatrix&) ;
320 static const double TINY;
345 Hist(
string titleIn,
int nBinIn = 100,
double xMinIn = 0.,
346 double xMaxIn = 1.) {
347 book(titleIn, nBinIn, xMinIn, xMaxIn);}
349 : title(h.title), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
350 xMax(h.xMax), dx(h.dx), under(h.under), inside(h.inside),
351 over(h.over), res(h.res) { }
353 : title(titleIn), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
354 xMax(h.xMax), dx(h.dx), under(h.under), inside(h.inside),
355 over(h.over), res(h.res) { }
356 Hist& operator=(
const Hist& h) {
if(
this != &h) {
357 nBin = h.nBin; nFill = h.nFill; xMin = h.xMin; xMax = h.xMax;
358 dx = h.dx; under = h.under; inside = h.inside; over = h.over;
359 res = h.res; }
return *
this; }
362 void book(
string titleIn =
" ",
int nBinIn = 100,
double xMinIn = 0.,
363 double xMaxIn = 1.) ;
366 void name(
string titleIn =
" ") {title = titleIn; }
372 void fill(
double x,
double w = 1.) ;
375 friend ostream& operator<<(ostream& os,
const Hist& h) ;
378 void table(ostream& os = cout)
const ;
379 void table(
string fileName)
const {
380 ofstream streamName(fileName.c_str()); table(streamName); }
383 friend void table(
const Hist& h1,
const Hist& h2, ostream& os) ;
384 friend void table(
const Hist& h1,
const Hist& h2,
string fileName) ;
387 double getBinContent(
int iBin) ;
390 int getEntries() {
return nFill; }
393 bool sameSize(
const Hist& h)
const ;
396 void takeLog(
bool tenLog =
true) ;
406 Hist& operator+=(
double f) ;
407 Hist& operator-=(
double f) ;
408 Hist& operator*=(
double f) ;
409 Hist& operator/=(
double f) ;
412 friend Hist operator+(
double f,
const Hist& h1);
413 friend Hist operator+(
const Hist& h1,
double f);
415 friend Hist operator-(
double f,
const Hist& h1);
416 friend Hist operator-(
const Hist& h1,
double f);
418 friend Hist operator*(
double f,
const Hist& h1);
419 friend Hist operator*(
const Hist& h1,
double f);
421 friend Hist operator/(
double f,
const Hist& h1);
422 friend Hist operator/(
const Hist& h1,
double f);
428 static const int NBINMAX, NCOLMAX, NLINES;
429 static const double TOLERANCE, TINY, SMALLFRAC, DYAC[];
430 static const char NUMBER[];
435 double xMin, xMax, dx, under, inside, over;
445 ostream& operator<<(ostream& os,
const Hist& h) ;
448 void table(
const Hist& h1,
const Hist& h2, ostream& os = cout) ;
449 void table(
const Hist& h1,
const Hist& h2,
string fileName) ;
452 Hist operator+(
double f,
const Hist& h1);
453 Hist operator+(
const Hist& h1,
double f);
455 Hist operator-(
double f,
const Hist& h1);
456 Hist operator-(
const Hist& h1,
double f);
458 Hist operator*(
double f,
const Hist& h1);
459 Hist operator*(
const Hist& h1,
double f);
461 Hist operator/(
double f,
const Hist& h1);
462 Hist operator/(
const Hist& h1,
double f);
469 #endif // end Pythia8_Basics_H