Line data Source code
1 : // Basics.h is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 Torbjorn Sjostrand.
3 : // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 : // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 :
6 : // Header file for basic, often-used helper classes.
7 : // RndmEngine: base class for external random number generators.
8 : // Rndm: random number generator.
9 : // Vec4: simple four-vectors.
10 : // RotBstMatrix: matrices encoding rotations and boosts of Vec4 objects.
11 : // Hist: simple one-dimensional histograms.
12 :
13 : #ifndef Pythia8_Basics_H
14 : #define Pythia8_Basics_H
15 :
16 : #include "Pythia8/PythiaStdlib.h"
17 :
18 : namespace Pythia8 {
19 :
20 : //==========================================================================
21 :
22 : // RndmEngine is the base class for external random number generators.
23 : // There is only one pure virtual method, that should do the generation.
24 :
25 : class RndmEngine {
26 :
27 : public:
28 :
29 : // Destructor.
30 : virtual ~RndmEngine() {}
31 :
32 : // A pure virtual method, wherein the derived class method
33 : // generates a random number uniformly distributed between 1 and 1.
34 : virtual double flat() = 0;
35 :
36 : };
37 :
38 : //==========================================================================
39 :
40 : // Rndm class.
41 : // This class handles random number generation according to the
42 : // Marsaglia-Zaman-Tsang algorithm.
43 :
44 : class Rndm {
45 :
46 : public:
47 :
48 : // Constructors.
49 0 : Rndm() : initRndm(false), seedSave(0), sequence(0),
50 0 : useExternalRndm(false), rndmEngPtr(0) { }
51 : Rndm(int seedIn) : initRndm(false), seedSave(0), sequence(0),
52 : useExternalRndm(false), rndmEngPtr(0) { init(seedIn);}
53 :
54 : // Possibility to pass in pointer for external random number generation.
55 : bool rndmEnginePtr( RndmEngine* rndmEngPtrIn);
56 :
57 : // Initialize, normally at construction or in first call.
58 : void init(int seedIn = 0) ;
59 :
60 : // Generate next random number uniformly between 0 and 1.
61 : double flat() ;
62 :
63 : // Generate random numbers according to exp(-x).
64 0 : double exp() { return -log(flat()) ;}
65 :
66 : // Generate random numbers according to x * exp(-x).
67 : double xexp() { return -log(flat() * flat()) ;}
68 :
69 : // Generate random numbers according to exp(-x^2/2).
70 0 : double gauss() {return sqrt(-2. * log(flat())) * cos(M_PI * flat());}
71 :
72 : // Generate two random numbers according to exp(-x^2/2-y^2/2).
73 0 : pair<double, double> gauss2() {double r = sqrt(-2. * log(flat()));
74 0 : double phi = 2. * M_PI * flat();
75 0 : return pair<double, double>(r * sin(phi), r * cos(phi));}
76 :
77 : // Pick one option among vector of (positive) probabilities.
78 : int pick(const vector<double>& prob) ;
79 :
80 : // Save or read current state to or from a binary file.
81 : bool dumpState(string fileName);
82 : bool readState(string fileName);
83 :
84 : private:
85 :
86 : // Default random number sequence.
87 : static const int DEFAULTSEED;
88 :
89 : // State of the random number generator.
90 : bool initRndm;
91 : int i97, j97, seedSave;
92 : long sequence;
93 : double u[97], c, cd, cm;
94 :
95 : // Pointer for external random number generation.
96 : bool useExternalRndm;
97 : RndmEngine* rndmEngPtr;
98 :
99 : };
100 :
101 : //==========================================================================
102 :
103 : // Forward reference to RotBstMatrix class; needed in Vec4 class.
104 : class RotBstMatrix;
105 :
106 : //--------------------------------------------------------------------------
107 :
108 : // Vec4 class.
109 : // This class implements four-vectors, in energy-momentum space.
110 : // (But can equally well be used to hold space-time four-vectors.)
111 :
112 : class Vec4 {
113 :
114 : public:
115 :
116 : // Constructors.
117 : Vec4(double xIn = 0., double yIn = 0., double zIn = 0., double tIn = 0.)
118 0 : : xx(xIn), yy(yIn), zz(zIn), tt(tIn) { }
119 0 : Vec4(const Vec4& v) : xx(v.xx), yy(v.yy), zz(v.zz), tt(v.tt) { }
120 0 : Vec4& operator=(const Vec4& v) { if (this != &v) { xx = v.xx; yy = v.yy;
121 0 : zz = v.zz; tt = v.tt; } return *this; }
122 0 : Vec4& operator=(double value) { xx = value; yy = value; zz = value;
123 0 : tt = value; return *this; }
124 :
125 : // Member functions for input.
126 : void reset() {xx = 0.; yy = 0.; zz = 0.; tt = 0.;}
127 : void p(double xIn, double yIn, double zIn, double tIn)
128 0 : {xx = xIn; yy = yIn; zz = zIn; tt = tIn;}
129 0 : void p(Vec4 pIn) {xx = pIn.xx; yy = pIn.yy; zz = pIn.zz; tt = pIn.tt;}
130 0 : void px(double xIn) {xx = xIn;}
131 0 : void py(double yIn) {yy = yIn;}
132 0 : void pz(double zIn) {zz = zIn;}
133 0 : void e(double tIn) {tt = tIn;}
134 :
135 : // Member functions for output.
136 0 : double px() const {return xx;}
137 0 : double py() const {return yy;}
138 0 : double pz() const {return zz;}
139 0 : double e() const {return tt;}
140 : double& operator[](int i) {
141 0 : if (i == 1) return xx;
142 0 : else if (i == 2) return yy;
143 0 : else if (i == 3) return zz;
144 0 : else return tt;
145 0 : }
146 0 : double mCalc() const {double temp = tt*tt - xx*xx - yy*yy - zz*zz;
147 0 : return (temp >= 0.) ? sqrt(temp) : -sqrt(-temp);}
148 0 : double m2Calc() const {return tt*tt - xx*xx - yy*yy - zz*zz;}
149 0 : double pT() const {return sqrt(xx*xx + yy*yy);}
150 0 : double pT2() const {return xx*xx + yy*yy;}
151 0 : double pAbs() const {return sqrt(xx*xx + yy*yy + zz*zz);}
152 0 : double pAbs2() const {return xx*xx + yy*yy + zz*zz;}
153 : double eT() const {double temp = xx*xx + yy*yy;
154 : return tt * sqrt( temp / (temp + zz*zz) );}
155 : double eT2() const {double temp = xx*xx + yy*yy;
156 : return tt*tt * temp / (temp + zz*zz);}
157 0 : double theta() const {return atan2(sqrt(xx*xx + yy*yy), zz);}
158 0 : double phi() const {return atan2(yy,xx);}
159 : double thetaXZ() const {return atan2(xx,zz);}
160 0 : double pPos() const {return tt + zz;}
161 0 : double pNeg() const {return tt - zz;}
162 0 : double rap() const {return 0.5 * log( (tt + zz) / (tt - zz) );}
163 0 : double eta() const {double xyz = sqrt(xx*xx + yy*yy + zz*zz);
164 0 : if (xx != 0 || yy != 0) return 0.5 * log( (xyz + zz) / (xyz - zz) );
165 0 : else return 1e6;
166 0 : }
167 :
168 : // Member functions that perform operations.
169 0 : void rescale3(double fac) {xx *= fac; yy *= fac; zz *= fac;}
170 0 : void rescale4(double fac) {xx *= fac; yy *= fac; zz *= fac; tt *= fac;}
171 0 : void flip3() {xx = -xx; yy = -yy; zz = -zz;}
172 : void flip4() {xx = -xx; yy = -yy; zz = -zz; tt = -tt;}
173 : void rot(double thetaIn, double phiIn);
174 : void rotaxis(double phiIn, double nx, double ny, double nz);
175 : void rotaxis(double phiIn, const Vec4& n);
176 : void bst(double betaX, double betaY, double betaZ);
177 : void bst(double betaX, double betaY, double betaZ, double gamma);
178 : void bst(const Vec4& pIn);
179 : void bst(const Vec4& pIn, double mIn);
180 : void bstback(const Vec4& pIn);
181 : void bstback(const Vec4& pIn, double mIn);
182 : void rotbst(const RotBstMatrix& M);
183 :
184 : // Operator overloading with member functions
185 0 : Vec4 operator-() {Vec4 tmp; tmp.xx = -xx; tmp.yy = -yy; tmp.zz = -zz;
186 0 : tmp.tt = -tt; return tmp;}
187 0 : Vec4& operator+=(const Vec4& v) {xx += v.xx; yy += v.yy; zz += v.zz;
188 0 : tt += v.tt; return *this;}
189 0 : Vec4& operator-=(const Vec4& v) {xx -= v.xx; yy -= v.yy; zz -= v.zz;
190 0 : tt -= v.tt; return *this;}
191 0 : Vec4& operator*=(double f) {xx *= f; yy *= f; zz *= f;
192 0 : tt *= f; return *this;}
193 0 : Vec4& operator/=(double f) {xx /= f; yy /= f; zz /= f;
194 0 : tt /= f; return *this;}
195 :
196 : // Operator overloading with friends
197 : friend Vec4 operator+(const Vec4& v1, const Vec4& v2);
198 : friend Vec4 operator-(const Vec4& v1, const Vec4& v2);
199 : friend Vec4 operator*(double f, const Vec4& v1);
200 : friend Vec4 operator*(const Vec4& v1, double f);
201 : friend Vec4 operator/(const Vec4& v1, double f);
202 : friend double operator*(const Vec4& v1, const Vec4& v2);
203 :
204 : // Invariant mass of a pair and its square.
205 : friend double m(const Vec4& v1, const Vec4& v2);
206 : friend double m2(const Vec4& v1, const Vec4& v2);
207 :
208 : // Scalar and cross product of 3-vector parts.
209 : friend double dot3(const Vec4& v1, const Vec4& v2);
210 : friend Vec4 cross3(const Vec4& v1, const Vec4& v2);
211 :
212 : // theta is polar angle between v1 and v2.
213 : friend double theta(const Vec4& v1, const Vec4& v2);
214 : friend double costheta(const Vec4& v1, const Vec4& v2);
215 :
216 : // phi is azimuthal angle between v1 and v2 around z axis.
217 : friend double phi(const Vec4& v1, const Vec4& v2);
218 : friend double cosphi(const Vec4& v1, const Vec4& v2);
219 :
220 : // phi is azimuthal angle between v1 and v2 around n axis.
221 : friend double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
222 : friend double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
223 :
224 : // R is distance in cylindrical (y/eta, phi) coordinates.
225 : friend double RRapPhi(const Vec4& v1, const Vec4& v2);
226 : friend double REtaPhi(const Vec4& v1, const Vec4& v2);
227 :
228 : // Print a four-vector.
229 : friend ostream& operator<<(ostream&, const Vec4& v) ;
230 :
231 : private:
232 :
233 : // Constants: could only be changed in the code itself.
234 : static const double TINY;
235 :
236 : // The four-vector data members.
237 : double xx, yy, zz, tt;
238 :
239 : };
240 :
241 : //--------------------------------------------------------------------------
242 :
243 : // Namespace function declarations; friends of Vec4 class.
244 :
245 : // Implementation of operator overloading with friends.
246 :
247 : inline Vec4 operator+(const Vec4& v1, const Vec4& v2)
248 0 : {Vec4 v = v1 ; return v += v2;}
249 :
250 : inline Vec4 operator-(const Vec4& v1, const Vec4& v2)
251 0 : {Vec4 v = v1 ; return v -= v2;}
252 :
253 : inline Vec4 operator*(double f, const Vec4& v1)
254 0 : {Vec4 v = v1; return v *= f;}
255 :
256 : inline Vec4 operator*(const Vec4& v1, double f)
257 0 : {Vec4 v = v1; return v *= f;}
258 :
259 : inline Vec4 operator/(const Vec4& v1, double f)
260 0 : {Vec4 v = v1; return v /= f;}
261 :
262 : inline double operator*(const Vec4& v1, const Vec4& v2)
263 0 : {return v1.tt*v2.tt - v1.xx*v2.xx - v1.yy*v2.yy - v1.zz*v2.zz;}
264 :
265 : // Invariant mass of a pair and its square.
266 : double m(const Vec4& v1, const Vec4& v2);
267 : double m2(const Vec4& v1, const Vec4& v2);
268 :
269 : // Scalar and cross product of 3-vector parts.
270 : double dot3(const Vec4& v1, const Vec4& v2);
271 : Vec4 cross3(const Vec4& v1, const Vec4& v2);
272 :
273 : // theta is polar angle between v1 and v2.
274 : double theta(const Vec4& v1, const Vec4& v2);
275 : double costheta(const Vec4& v1, const Vec4& v2);
276 :
277 : // phi is azimuthal angle between v1 and v2 around z axis.
278 : double phi(const Vec4& v1, const Vec4& v2);
279 : double cosphi(const Vec4& v1, const Vec4& v2);
280 :
281 : // phi is azimuthal angle between v1 and v2 around n axis.
282 : double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
283 : double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
284 :
285 : // R is distance in cylindrical (y/eta, phi) coordinates.
286 : double RRapPhi(const Vec4& v1, const Vec4& v2);
287 : double REtaPhi(const Vec4& v1, const Vec4& v2);
288 :
289 : // Print a four-vector.
290 : ostream& operator<<(ostream&, const Vec4& v) ;
291 :
292 : //==========================================================================
293 :
294 : // RotBstMatrix class.
295 : // This class implements 4 * 4 matrices that encode an arbitrary combination
296 : // of rotations and boosts, that can be applied to Vec4 four-vectors.
297 :
298 : class RotBstMatrix {
299 :
300 : public:
301 :
302 : // Constructors.
303 0 : RotBstMatrix() {for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j)
304 0 : { M[i][j] = (i==j) ? 1. : 0.; } } }
305 0 : RotBstMatrix(const RotBstMatrix& Min) {
306 0 : for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
307 0 : M[i][j] = Min.M[i][j]; } } }
308 0 : RotBstMatrix& operator=(const RotBstMatrix& Min) {if (this != &Min) {
309 0 : for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
310 0 : M[i][j] = Min.M[i][j]; } } } return *this; }
311 :
312 : // Member functions.
313 : void rot(double = 0., double = 0.);
314 : void rot(const Vec4& p);
315 : void bst(double = 0., double = 0., double = 0.);
316 : void bst(const Vec4&);
317 : void bstback(const Vec4&);
318 : void bst(const Vec4&, const Vec4&);
319 : void toCMframe(const Vec4&, const Vec4&);
320 : void fromCMframe(const Vec4&, const Vec4&);
321 : void rotbst(const RotBstMatrix&);
322 : void invert();
323 : void reset();
324 :
325 : // Crude estimate deviation from unit matrix.
326 : double deviation() const;
327 :
328 : // Print a transformation matrix.
329 : friend ostream& operator<<(ostream&, const RotBstMatrix&) ;
330 :
331 : // Private members to be accessible from Vec4.
332 : friend class Vec4;
333 :
334 : private:
335 :
336 : // Constants: could only be changed in the code itself.
337 : static const double TINY;
338 :
339 : // The rotation-and-boost matrix data members.
340 : double M[4][4];
341 :
342 : };
343 :
344 : //--------------------------------------------------------------------------
345 :
346 : // Namespace function declaration; friend of RotBstMatrix class.
347 :
348 : // Print a transformation matrix.
349 : ostream& operator<<(ostream&, const RotBstMatrix&) ;
350 :
351 : //==========================================================================
352 :
353 : // Hist class.
354 : // This class handles a single histogram at a time.
355 :
356 0 : class Hist{
357 :
358 : public:
359 :
360 : // Constructors, including copy constructors.
361 : Hist() {;}
362 : Hist(string titleIn, int nBinIn = 100, double xMinIn = 0.,
363 : double xMaxIn = 1.) {
364 : book(titleIn, nBinIn, xMinIn, xMaxIn);}
365 : Hist(const Hist& h)
366 0 : : title(h.title), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
367 0 : xMax(h.xMax), dx(h.dx), under(h.under), inside(h.inside),
368 0 : over(h.over), res(h.res) { }
369 : Hist(string titleIn, const Hist& h)
370 : : title(titleIn), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
371 : xMax(h.xMax), dx(h.dx), under(h.under), inside(h.inside),
372 : over(h.over), res(h.res) { }
373 : Hist& operator=(const Hist& h) { if(this != &h) {
374 : nBin = h.nBin; nFill = h.nFill; xMin = h.xMin; xMax = h.xMax;
375 : dx = h.dx; under = h.under; inside = h.inside; over = h.over;
376 : res = h.res; } return *this; }
377 :
378 : // Book a histogram.
379 : void book(string titleIn = " ", int nBinIn = 100, double xMinIn = 0.,
380 : double xMaxIn = 1.) ;
381 :
382 : // Set title of a histogram.
383 : void name(string titleIn = " ") {title = titleIn; }
384 :
385 : // Reset bin contents.
386 : void null() ;
387 :
388 : // Fill bin with weight.
389 : void fill(double x, double w = 1.) ;
390 :
391 : // Print a histogram with overloaded << operator.
392 : friend ostream& operator<<(ostream& os, const Hist& h) ;
393 :
394 : // Print histogram contents as a table (e.g. for Gnuplot).
395 : void table(ostream& os = cout, bool printOverUnder = false,
396 : bool xMidBin = true) const ;
397 : void table(string fileName, bool printOverUnder = false,
398 : bool xMidBin = true) const { ofstream streamName(fileName.c_str());
399 : table(streamName, printOverUnder, xMidBin);}
400 :
401 : // Print a table out of two histograms with same x axis.
402 : friend void table(const Hist& h1, const Hist& h2, ostream& os,
403 : bool printOverUnder, bool xMidBin) ;
404 : friend void table(const Hist& h1, const Hist& h2, string fileName,
405 : bool printOverUnder, bool xMidBin) ;
406 :
407 : // Return content of specific bin: 0 gives underflow and nBin+1 overflow.
408 : double getBinContent(int iBin) const;
409 :
410 : // Return number of entries.
411 : int getEntries() const {return nFill; }
412 :
413 : // Check whether another histogram has same size and limits.
414 : bool sameSize(const Hist& h) const ;
415 :
416 : // Take logarithm (base 10 or e) of bin contents.
417 : void takeLog(bool tenLog = true) ;
418 :
419 : // Take square root of bin contents.
420 : void takeSqrt() ;
421 :
422 : // Operator overloading with member functions
423 : Hist& operator+=(const Hist& h) ;
424 : Hist& operator-=(const Hist& h) ;
425 : Hist& operator*=(const Hist& h) ;
426 : Hist& operator/=(const Hist& h) ;
427 : Hist& operator+=(double f) ;
428 : Hist& operator-=(double f) ;
429 : Hist& operator*=(double f) ;
430 : Hist& operator/=(double f) ;
431 :
432 : // Operator overloading with friends
433 : friend Hist operator+(double f, const Hist& h1);
434 : friend Hist operator+(const Hist& h1, double f);
435 : friend Hist operator+(const Hist& h1, const Hist& h2);
436 : friend Hist operator-(double f, const Hist& h1);
437 : friend Hist operator-(const Hist& h1, double f);
438 : friend Hist operator-(const Hist& h1, const Hist& h2);
439 : friend Hist operator*(double f, const Hist& h1);
440 : friend Hist operator*(const Hist& h1, double f);
441 : friend Hist operator*(const Hist& h1, const Hist& h2);
442 : friend Hist operator/(double f, const Hist& h1);
443 : friend Hist operator/(const Hist& h1, double f);
444 : friend Hist operator/(const Hist& h1, const Hist& h2);
445 :
446 : private:
447 :
448 : // Constants: could only be changed in the code itself.
449 : static const int NBINMAX, NCOLMAX, NLINES;
450 : static const double TOLERANCE, TINY, LARGE, SMALLFRAC, DYAC[];
451 : static const char NUMBER[];
452 :
453 : // Properties and contents of a histogram.
454 : string title;
455 : int nBin, nFill;
456 : double xMin, xMax, dx, under, inside, over;
457 : vector<double> res;
458 :
459 : };
460 :
461 : //--------------------------------------------------------------------------
462 :
463 : // Namespace function declarations; friends of Hist class.
464 :
465 : // Print a histogram with overloaded << operator.
466 : ostream& operator<<(ostream& os, const Hist& h) ;
467 :
468 : // Print a table out of two histograms with same x axis.
469 : void table(const Hist& h1, const Hist& h2, ostream& os = cout,
470 : bool printOverUnder = false, bool xMidBin = true) ;
471 : void table(const Hist& h1, const Hist& h2, string fileName,
472 : bool printOverUnder = false, bool xMidBin = true) ;
473 :
474 : // Operator overloading with friends
475 : Hist operator+(double f, const Hist& h1);
476 : Hist operator+(const Hist& h1, double f);
477 : Hist operator+(const Hist& h1, const Hist& h2);
478 : Hist operator-(double f, const Hist& h1);
479 : Hist operator-(const Hist& h1, double f);
480 : Hist operator-(const Hist& h1, const Hist& h2);
481 : Hist operator*(double f, const Hist& h1);
482 : Hist operator*(const Hist& h1, double f);
483 : Hist operator*(const Hist& h1, const Hist& h2);
484 : Hist operator/(double f, const Hist& h1);
485 : Hist operator/(const Hist& h1, double f);
486 : Hist operator/(const Hist& h1, const Hist& h2);
487 :
488 : //==========================================================================
489 :
490 : } // end namespace Pythia8
491 :
492 : #endif // end Pythia8_Basics_H
|
