StRoot  1
Basics.h
1 // Basics.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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  Rndm() : initRndm(false), seedSave(0), sequence(0),
50  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  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  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  pair<double, double> gauss2() {double r = sqrt(-2. * log(flat()));
74  double phi = 2. * M_PI * flat();
75  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);
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  : xx(xIn), yy(yIn), zz(zIn), tt(tIn) { }
119  Vec4(const Vec4& v) : xx(v.xx), yy(v.yy), zz(v.zz), tt(v.tt) { }
120  Vec4& operator=(const Vec4& v) { if (this != &v) { xx = v.xx; yy = v.yy;
121  zz = v.zz; tt = v.tt; } return *this; }
122  Vec4& operator=(double value) { xx = value; yy = value; zz = value;
123  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  {xx = xIn; yy = yIn; zz = zIn; tt = tIn;}
129  void p(Vec4 pIn) {xx = pIn.xx; yy = pIn.yy; zz = pIn.zz; tt = pIn.tt;}
130  void px(double xIn) {xx = xIn;}
131  void py(double yIn) {yy = yIn;}
132  void pz(double zIn) {zz = zIn;}
133  void e(double tIn) {tt = tIn;}
134
135  // Member functions for output.
136  double px() const {return xx;}
137  double py() const {return yy;}
138  double pz() const {return zz;}
139  double e() const {return tt;}
140  double mCalc() const {double temp = tt*tt - xx*xx - yy*yy - zz*zz;
141  return (temp >= 0.) ? sqrt(temp) : -sqrt(-temp);}
142  double m2Calc() const {return tt*tt - xx*xx - yy*yy - zz*zz;}
143  double pT() const {return sqrt(xx*xx + yy*yy);}
144  double pT2() const {return xx*xx + yy*yy;}
145  double pAbs() const {return sqrt(xx*xx + yy*yy + zz*zz);}
146  double pAbs2() const {return xx*xx + yy*yy + zz*zz;}
147  double eT() const {double temp = xx*xx + yy*yy;
148  return tt * sqrt( temp / (temp + zz*zz) );}
149  double eT2() const {double temp = xx*xx + yy*yy;
150  return tt*tt * temp / (temp + zz*zz);}
151  double theta() const {return atan2(sqrt(xx*xx + yy*yy), zz);}
152  double phi() const {return atan2(yy,xx);}
153  double thetaXZ() const {return atan2(xx,zz);}
154  double pPos() const {return tt + zz;}
155  double pNeg() const {return tt - zz;}
156  double rap() const {return 0.5 * log( (tt + zz) / (tt - zz) );}
157  double eta() const {double xyz = sqrt(xx*xx + yy*yy + zz*zz);
158  return 0.5 * log( (xyz + zz) / (xyz - zz) );}
159
160  // Member functions that perform operations.
161  void rescale3(double fac) {xx *= fac; yy *= fac; zz *= fac;}
162  void rescale4(double fac) {xx *= fac; yy *= fac; zz *= fac; tt *= fac;}
163  void flip3() {xx = -xx; yy = -yy; zz = -zz;}
164  void flip4() {xx = -xx; yy = -yy; zz = -zz; tt = -tt;}
165  void rot(double thetaIn, double phiIn);
166  void rotaxis(double phiIn, double nx, double ny, double nz);
167  void rotaxis(double phiIn, const Vec4& n);
168  void bst(double betaX, double betaY, double betaZ);
169  void bst(double betaX, double betaY, double betaZ, double gamma);
170  void bst(const Vec4& pIn);
171  void bst(const Vec4& pIn, double mIn);
172  void bstback(const Vec4& pIn);
173  void bstback(const Vec4& pIn, double mIn);
174  void rotbst(const RotBstMatrix& M);
175
177  Vec4 operator-() {Vec4 tmp; tmp.xx = -xx; tmp.yy = -yy; tmp.zz = -zz;
178  tmp.tt = -tt; return tmp;}
179  Vec4& operator+=(const Vec4& v) {xx += v.xx; yy += v.yy; zz += v.zz;
180  tt += v.tt; return *this;}
181  Vec4& operator-=(const Vec4& v) {xx -= v.xx; yy -= v.yy; zz -= v.zz;
182  tt -= v.tt; return *this;}
183  Vec4& operator*=(double f) {xx *= f; yy *= f; zz *= f;
184  tt *= f; return *this;}
185  Vec4& operator/=(double f) {xx /= f; yy /= f; zz /= f;
186  tt /= f; return *this;}
187
189  friend Vec4 operator+(const Vec4& v1, const Vec4& v2);
190  friend Vec4 operator-(const Vec4& v1, const Vec4& v2);
191  friend Vec4 operator*(double f, const Vec4& v1);
192  friend Vec4 operator*(const Vec4& v1, double f);
193  friend Vec4 operator/(const Vec4& v1, double f);
194  friend double operator*(const Vec4& v1, const Vec4& v2);
195
196  // Invariant mass of a pair and its square.
197  friend double m(const Vec4& v1, const Vec4& v2);
198  friend double m2(const Vec4& v1, const Vec4& v2);
199
200  // Scalar and cross product of 3-vector parts.
201  friend double dot3(const Vec4& v1, const Vec4& v2);
202  friend Vec4 cross3(const Vec4& v1, const Vec4& v2);
203
204  // theta is polar angle between v1 and v2.
205  friend double theta(const Vec4& v1, const Vec4& v2);
206  friend double costheta(const Vec4& v1, const Vec4& v2);
207
208  // phi is azimuthal angle between v1 and v2 around z axis.
209  friend double phi(const Vec4& v1, const Vec4& v2);
210  friend double cosphi(const Vec4& v1, const Vec4& v2);
211
212  // phi is azimuthal angle between v1 and v2 around n axis.
213  friend double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
214  friend double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
215
216  // R is distance in cylindrical (y/eta, phi) coordinates.
217  friend double RRapPhi(const Vec4& v1, const Vec4& v2);
218  friend double REtaPhi(const Vec4& v1, const Vec4& v2);
219
220  // Print a four-vector.
221  friend ostream& operator<<(ostream&, const Vec4& v) ;
222
223 private:
224
225  // Constants: could only be changed in the code itself.
226  static const double TINY;
227
228  // The four-vector data members.
229  double xx, yy, zz, tt;
230
231 };
232
233 //--------------------------------------------------------------------------
234
235 // Namespace function declarations; friends of Vec4 class.
236
238
239 inline Vec4 operator+(const Vec4& v1, const Vec4& v2)
240  {Vec4 v = v1 ; return v += v2;}
241
242 inline Vec4 operator-(const Vec4& v1, const Vec4& v2)
243  {Vec4 v = v1 ; return v -= v2;}
244
245 inline Vec4 operator*(double f, const Vec4& v1)
246  {Vec4 v = v1; return v *= f;}
247
248 inline Vec4 operator*(const Vec4& v1, double f)
249  {Vec4 v = v1; return v *= f;}
250
251 inline Vec4 operator/(const Vec4& v1, double f)
252  {Vec4 v = v1; return v /= f;}
253
254 inline double operator*(const Vec4& v1, const Vec4& v2)
255  {return v1.tt*v2.tt - v1.xx*v2.xx - v1.yy*v2.yy - v1.zz*v2.zz;}
256
257 // Invariant mass of a pair and its square.
258 double m(const Vec4& v1, const Vec4& v2);
259 double m2(const Vec4& v1, const Vec4& v2);
260
261 // Scalar and cross product of 3-vector parts.
262 double dot3(const Vec4& v1, const Vec4& v2);
263 Vec4 cross3(const Vec4& v1, const Vec4& v2);
264
265 // theta is polar angle between v1 and v2.
266 double theta(const Vec4& v1, const Vec4& v2);
267 double costheta(const Vec4& v1, const Vec4& v2);
268
269 // phi is azimuthal angle between v1 and v2 around z axis.
270 double phi(const Vec4& v1, const Vec4& v2);
271 double cosphi(const Vec4& v1, const Vec4& v2);
272
273 // phi is azimuthal angle between v1 and v2 around n axis.
274 double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
275 double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
276
277 // R is distance in cylindrical (y/eta, phi) coordinates.
278 double RRapPhi(const Vec4& v1, const Vec4& v2);
279 double REtaPhi(const Vec4& v1, const Vec4& v2);
280
281 // Print a four-vector.
282 ostream& operator<<(ostream&, const Vec4& v) ;
283
284 //==========================================================================
285
286 // RotBstMatrix class.
287 // This class implements 4 * 4 matrices that encode an arbitrary combination
288 // of rotations and boosts, that can be applied to Vec4 four-vectors.
289
290 class RotBstMatrix {
291
292 public:
293
294  // Constructors.
295  RotBstMatrix() {for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j)
296  { M[i][j] = (i==j) ? 1. : 0.; } } }
297  RotBstMatrix(const RotBstMatrix& Min) {
298  for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
299  M[i][j] = Min.M[i][j]; } } }
300  RotBstMatrix& operator=(const RotBstMatrix& Min) {if (this != &Min) {
301  for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
302  M[i][j] = Min.M[i][j]; } } } return *this; }
303
304  // Member functions.
305  void rot(double = 0., double = 0.);
306  void rot(const Vec4& p);
307  void bst(double = 0., double = 0., double = 0.);
308  void bst(const Vec4&);
309  void bstback(const Vec4&);
310  void bst(const Vec4&, const Vec4&);
311  void toCMframe(const Vec4&, const Vec4&);
312  void fromCMframe(const Vec4&, const Vec4&);
313  void rotbst(const RotBstMatrix&);
314  void invert();
315  void reset();
316
317  // Crude estimate deviation from unit matrix.
318  double deviation() const;
319
320  // Print a transformation matrix.
321  friend ostream& operator<<(ostream&, const RotBstMatrix&) ;
322
323  // Private members to be accessible from Vec4.
324  friend class Vec4;
325
326 private:
327
328  // Constants: could only be changed in the code itself.
329  static const double TINY;
330
331  // The rotation-and-boost matrix data members.
332  double M[4][4];
333
334 };
335
336 //--------------------------------------------------------------------------
337
338 // Namespace function declaration; friend of RotBstMatrix class.
339
340 // Print a transformation matrix.
341 ostream& operator<<(ostream&, const RotBstMatrix&) ;
342
343 //==========================================================================
344
345 // Hist class.
346 // This class handles a single histogram at a time.
347
348 class Hist{
349
350 public:
351
352  // Constructors, including copy constructors.
353  Hist() {;}
354  Hist(string titleIn, int nBinIn = 100, double xMinIn = 0.,
355  double xMaxIn = 1.) {
356  book(titleIn, nBinIn, xMinIn, xMaxIn);}
357  Hist(const Hist& h)
358  : title(h.title), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
359  xMax(h.xMax), dx(h.dx), under(h.under), inside(h.inside),
360  over(h.over), res(h.res) { }
361  Hist(string titleIn, const Hist& h)
362  : title(titleIn), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
363  xMax(h.xMax), dx(h.dx), under(h.under), inside(h.inside),
364  over(h.over), res(h.res) { }
365  Hist& operator=(const Hist& h) { if(this != &h) {
366  nBin = h.nBin; nFill = h.nFill; xMin = h.xMin; xMax = h.xMax;
367  dx = h.dx; under = h.under; inside = h.inside; over = h.over;
368  res = h.res; } return *this; }
369
370  // Book a histogram.
371  void book(string titleIn = " ", int nBinIn = 100, double xMinIn = 0.,
372  double xMaxIn = 1.) ;
373
374  // Set title of a histogram.
375  void name(string titleIn = " ") {title = titleIn; }
376
377  // Reset bin contents.
378  void null() ;
379
380  // Fill bin with weight.
381  void fill(double x, double w = 1.) ;
382
383  // Print a histogram with overloaded << operator.
384  friend ostream& operator<<(ostream& os, const Hist& h) ;
385
386  // Print histogram contents as a table (e.g. for Gnuplot).
387  void table(ostream& os = cout, bool printOverUnder = false,
388  bool xMidBin = true) const ;
389  void table(string fileName, bool printOverUnder = false,
390  bool xMidBin = true) const { ofstream streamName(fileName.c_str());
391  table(streamName, printOverUnder, xMidBin);}
392
393  // Print a table out of two histograms with same x axis.
394  friend void table(const Hist& h1, const Hist& h2, ostream& os,
395  bool printOverUnder, bool xMidBin) ;
396  friend void table(const Hist& h1, const Hist& h2, string fileName,
397  bool printOverUnder, bool xMidBin) ;
398
399  // Return content of specific bin: 0 gives underflow and nBin+1 overflow.
400  double getBinContent(int iBin) const;
401
402  // Return number of entries.
403  int getEntries() const {return nFill; }
404
405  // Check whether another histogram has same size and limits.
406  bool sameSize(const Hist& h) const ;
407
408  // Take logarithm (base 10 or e) of bin contents.
409  void takeLog(bool tenLog = true) ;
410
411  // Take square root of bin contents.
412  void takeSqrt() ;
413
415  Hist& operator+=(const Hist& h) ;
416  Hist& operator-=(const Hist& h) ;
417  Hist& operator*=(const Hist& h) ;
418  Hist& operator/=(const Hist& h) ;
419  Hist& operator+=(double f) ;
420  Hist& operator-=(double f) ;
421  Hist& operator*=(double f) ;
422  Hist& operator/=(double f) ;
423
425  friend Hist operator+(double f, const Hist& h1);
426  friend Hist operator+(const Hist& h1, double f);
427  friend Hist operator+(const Hist& h1, const Hist& h2);
428  friend Hist operator-(double f, const Hist& h1);
429  friend Hist operator-(const Hist& h1, double f);
430  friend Hist operator-(const Hist& h1, const Hist& h2);
431  friend Hist operator*(double f, const Hist& h1);
432  friend Hist operator*(const Hist& h1, double f);
433  friend Hist operator*(const Hist& h1, const Hist& h2);
434  friend Hist operator/(double f, const Hist& h1);
435  friend Hist operator/(const Hist& h1, double f);
436  friend Hist operator/(const Hist& h1, const Hist& h2);
437
438 private:
439
440  // Constants: could only be changed in the code itself.
441  static const int NBINMAX, NCOLMAX, NLINES;
442  static const double TOLERANCE, TINY, LARGE, SMALLFRAC, DYAC[];
443  static const char NUMBER[];
444
445  // Properties and contents of a histogram.
446  string title;
447  int nBin, nFill;
448  double xMin, xMax, dx, under, inside, over;
449  vector<double> res;
450
451 };
452
453 //--------------------------------------------------------------------------
454
455 // Namespace function declarations; friends of Hist class.
456
457 // Print a histogram with overloaded << operator.
458 ostream& operator<<(ostream& os, const Hist& h) ;
459
460 // Print a table out of two histograms with same x axis.
461 void table(const Hist& h1, const Hist& h2, ostream& os = cout,
462  bool printOverUnder = false, bool xMidBin = true) ;
463 void table(const Hist& h1, const Hist& h2, string fileName,
464  bool printOverUnder = false, bool xMidBin = true) ;
465
467 Hist operator+(double f, const Hist& h1);
468 Hist operator+(const Hist& h1, double f);
469 Hist operator+(const Hist& h1, const Hist& h2);
470 Hist operator-(double f, const Hist& h1);
471 Hist operator-(const Hist& h1, double f);
472 Hist operator-(const Hist& h1, const Hist& h2);
473 Hist operator*(double f, const Hist& h1);
474 Hist operator*(const Hist& h1, double f);
475 Hist operator*(const Hist& h1, const Hist& h2);
476 Hist operator/(double f, const Hist& h1);
477 Hist operator/(const Hist& h1, double f);
478 Hist operator/(const Hist& h1, const Hist& h2);
479
480 //==========================================================================
481
482 } // end namespace Pythia8
483
484 #endif // end Pythia8_Basics_H