StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Basics.h
1 // Basics.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2012 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 "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  // A pure virtual method, wherein the derived class method
30  // generates a random number uniformly distributed between 1 and 1.
31  virtual double flat() = 0;
32 
33 protected:
34 
35  // Destructor.
36  virtual ~RndmEngine() {}
37 
38 };
39 
40 //==========================================================================
41 
42 // Rndm class.
43 // This class handles random number generation according to the
44 // Marsaglia-Zaman-Tsang algorithm.
45 
46 class Rndm {
47 
48 public:
49 
50  // Constructors.
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);}
55 
56  // Possibility to pass in pointer for external random number generation.
57  bool rndmEnginePtr( RndmEngine* rndmEngPtrIn);
58 
59  // Initialize, normally at construction or in first call.
60  void init(int seedIn = 0) ;
61 
62  // Generate next random number uniformly between 0 and 1.
63  double flat() ;
64 
65  // Generate random numbers according to exp(-x).
66  double exp() { return -log(flat()) ;}
67 
68  // Generate random numbers according to x * exp(-x).
69  double xexp() { return -log(flat() * flat()) ;}
70 
71  // Generate random numbers according to exp(-x^2/2).
72  double gauss() {return sqrt(-2. * log(flat())) * cos(M_PI * flat());}
73 
74  // Generate two random numbers according to exp(-x^2/2-y^2/2).
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));}
78 
79  // Pick one option among vector of (positive) probabilities.
80  int pick(const vector<double>& prob) ;
81 
82  // Save or read current state to or from a binary file.
83  bool dumpState(string fileName);
84  bool readState(string fileName);
85 
86 private:
87 
88  // Default random number sequence.
89  static const int DEFAULTSEED;
90 
91  // State of the random number generator.
92  bool initRndm;
93  int i97, j97, defaultSeed, seedSave;
94  long sequence;
95  double u[97], c, cd, cm;
96 
97  // Pointer for external random number generation.
98  bool useExternalRndm;
99  RndmEngine* rndmEngPtr;
100 
101 };
102 
103 //==========================================================================
104 
105 // Forward reference to RotBstMatrix class; needed in Vec4 class.
106 class RotBstMatrix;
107 
108 //--------------------------------------------------------------------------
109 
110 // Vec4 class.
111 // This class implements four-vectors, in energy-momentum space.
112 // (But can equally well be used to hold space-time four-vectors.)
113 
114 class Vec4 {
115 
116 public:
117 
118  // Constructors.
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; }
126 
127  // Member functions for input.
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;}
136 
137  // Member functions for output.
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;}
158 
159  // Member functions that perform operations.
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);
173  void rotbst(const RotBstMatrix& M);
174 
175  // Operator overloading with member functions
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;}
186 
187  // Operator overloading with friends
188  friend Vec4 operator+(const Vec4& v1, const Vec4& v2);
189  friend Vec4 operator-(const Vec4& v1, const Vec4& v2);
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);
194 
195  // Invariant mass of a pair and its square.
196  friend double m(const Vec4& v1, const Vec4& v2);
197  friend double m2(const Vec4& v1, const Vec4& v2);
198 
199  // Scalar and cross product of 3-vector parts.
200  friend double dot3(const Vec4& v1, const Vec4& v2);
201  friend Vec4 cross3(const Vec4& v1, const Vec4& v2);
202 
203  // theta is polar angle between v1 and v2.
204  friend double theta(const Vec4& v1, const Vec4& v2);
205  friend double costheta(const Vec4& v1, const Vec4& v2);
206 
207  // phi is azimuthal angle between v1 and v2 around z axis.
208  friend double phi(const Vec4& v1, const Vec4& v2);
209  friend double cosphi(const Vec4& v1, const Vec4& v2);
210 
211  // phi is azimuthal angle between v1 and v2 around n axis.
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);
214 
215  // Print a four-vector
216  friend ostream& operator<<(ostream&, const Vec4& v) ;
217 
218 private:
219 
220  // Constants: could only be changed in the code itself.
221  static const double TINY;
222 
223  // The four-vector data members.
224  double xx, yy, zz, tt;
225 
226 };
227 
228 //--------------------------------------------------------------------------
229 
230 // Namespace function declarations; friends of Vec4 class.
231 
232 // Implementation of operator overloading with friends.
233 
234 inline Vec4 operator+(const Vec4& v1, const Vec4& v2)
235  {Vec4 v = v1 ; return v += v2;}
236 
237 inline Vec4 operator-(const Vec4& v1, const Vec4& v2)
238  {Vec4 v = v1 ; return v -= v2;}
239 
240 inline Vec4 operator*(double f, const Vec4& v1)
241  {Vec4 v = v1; return v *= f;}
242 
243 inline Vec4 operator*(const Vec4& v1, double f)
244  {Vec4 v = v1; return v *= f;}
245 
246 inline Vec4 operator/(const Vec4& v1, double f)
247  {Vec4 v = v1; return v /= f;}
248 
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;}
251 
252 // Invariant mass of a pair and its square.
253 double m(const Vec4& v1, const Vec4& v2);
254 double m2(const Vec4& v1, const Vec4& v2);
255 
256 // Scalar and cross product of 3-vector parts.
257 double dot3(const Vec4& v1, const Vec4& v2);
258 Vec4 cross3(const Vec4& v1, const Vec4& v2);
259 
260 // theta is polar angle between v1 and v2.
261 double theta(const Vec4& v1, const Vec4& v2);
262 double costheta(const Vec4& v1, const Vec4& v2);
263 
264 // phi is azimuthal angle between v1 and v2 around z axis.
265 double phi(const Vec4& v1, const Vec4& v2);
266 double cosphi(const Vec4& v1, const Vec4& v2);
267 
268 // phi is azimuthal angle between v1 and v2 around n axis.
269 double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
270 double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
271 
272 // Print a four-vector.
273 ostream& operator<<(ostream&, const Vec4& v) ;
274 
275 //==========================================================================
276 
277 // RotBstMatrix class.
278 // This class implements 4 * 4 matrices that encode an arbitrary combination
279 // of rotations and boosts, that can be applied to Vec4 four-vectors.
280 
281 class RotBstMatrix {
282 
283 public:
284 
285  // Constructors.
286  RotBstMatrix() {for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j)
287  { M[i][j] = (i==j) ? 1. : 0.; } } }
288  RotBstMatrix(const RotBstMatrix& Min) {
289  for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
290  M[i][j] = Min.M[i][j]; } } }
291  RotBstMatrix& operator=(const RotBstMatrix& Min) {if (this != &Min) {
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; }
294 
295  // Member functions.
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&);
301  void bst(const Vec4&, const Vec4&);
302  void toCMframe(const Vec4&, const Vec4&);
303  void fromCMframe(const Vec4&, const Vec4&);
304  void rotbst(const RotBstMatrix&);
305  void invert();
306  void reset();
307 
308  // Crude estimate deviation from unit matrix.
309  double deviation() const;
310 
311  // Print a transformation matrix.
312  friend ostream& operator<<(ostream&, const RotBstMatrix&) ;
313 
314  // Private members to be accessible from Vec4.
315  friend class Vec4;
316 
317 private:
318 
319  // Constants: could only be changed in the code itself.
320  static const double TINY;
321 
322  // The rotation-and-boost matrix data members.
323  double M[4][4];
324 
325 };
326 
327 //--------------------------------------------------------------------------
328 
329 // Namespace function declaration; friend of RotBstMatrix class.
330 
331 // Print a transformation matrix.
332 ostream& operator<<(ostream&, const RotBstMatrix&) ;
333 
334 //==========================================================================
335 
336 // Hist class.
337 // This class handles a single histogram at a time.
338 
339 class Hist{
340 
341 public:
342 
343  // Constructors, including copy constructors.
344  Hist() {;}
345  Hist(string titleIn, int nBinIn = 100, double xMinIn = 0.,
346  double xMaxIn = 1.) {
347  book(titleIn, nBinIn, xMinIn, xMaxIn);}
348  Hist(const Hist& h)
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) { }
352  Hist(string titleIn, const Hist& h)
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; }
360 
361  // Book a histogram.
362  void book(string titleIn = " ", int nBinIn = 100, double xMinIn = 0.,
363  double xMaxIn = 1.) ;
364 
365  // Set title of a histogram.
366  void name(string titleIn = " ") {title = titleIn; }
367 
368  // Reset bin contents.
369  void null() ;
370 
371  // Fill bin with weight.
372  void fill(double x, double w = 1.) ;
373 
374  // Print a histogram with overloaded << operator.
375  friend ostream& operator<<(ostream& os, const Hist& h) ;
376 
377  // Print histogram contents as a table (e.g. for Gnuplot).
378  void table(ostream& os = cout) const ;
379  void table(string fileName) const {
380  ofstream streamName(fileName.c_str()); table(streamName); }
381 
382  // Print a table out of two histograms with same x axis.
383  friend void table(const Hist& h1, const Hist& h2, ostream& os) ;
384  friend void table(const Hist& h1, const Hist& h2, string fileName) ;
385 
386  // Return content of specific bin: -1 gives underflow and nBin overflow.
387  double getBinContent(int iBin) ;
388 
389  // Return number of entries
390  int getEntries() {return nFill; }
391 
392  // Check whether another histogram has same size and limits.
393  bool sameSize(const Hist& h) const ;
394 
395  // Take logarithm (base 10 or e) of bin contents.
396  void takeLog(bool tenLog = true) ;
397 
398  // Take square root of bin contents.
399  void takeSqrt() ;
400 
401  // Operator overloading with member functions
402  Hist& operator+=(const Hist& h) ;
403  Hist& operator-=(const Hist& h) ;
404  Hist& operator*=(const Hist& h) ;
405  Hist& operator/=(const Hist& h) ;
406  Hist& operator+=(double f) ;
407  Hist& operator-=(double f) ;
408  Hist& operator*=(double f) ;
409  Hist& operator/=(double f) ;
410 
411  // Operator overloading with friends
412  friend Hist operator+(double f, const Hist& h1);
413  friend Hist operator+(const Hist& h1, double f);
414  friend Hist operator+(const Hist& h1, const Hist& h2);
415  friend Hist operator-(double f, const Hist& h1);
416  friend Hist operator-(const Hist& h1, double f);
417  friend Hist operator-(const Hist& h1, const Hist& h2);
418  friend Hist operator*(double f, const Hist& h1);
419  friend Hist operator*(const Hist& h1, double f);
420  friend Hist operator*(const Hist& h1, const Hist& h2);
421  friend Hist operator/(double f, const Hist& h1);
422  friend Hist operator/(const Hist& h1, double f);
423  friend Hist operator/(const Hist& h1, const Hist& h2);
424 
425 private:
426 
427  // Constants: could only be changed in the code itself.
428  static const int NBINMAX, NCOLMAX, NLINES;
429  static const double TOLERANCE, TINY, SMALLFRAC, DYAC[];
430  static const char NUMBER[];
431 
432  // Properties and contents of a histogram.
433  string title;
434  int nBin, nFill;
435  double xMin, xMax, dx, under, inside, over;
436  vector<double> res;
437 
438 };
439 
440 //--------------------------------------------------------------------------
441 
442 // Namespace function declarations; friends of Hist class.
443 
444 // Print a histogram with overloaded << operator.
445 ostream& operator<<(ostream& os, const Hist& h) ;
446 
447 // Print a table out of two histograms with same x axis.
448 void table(const Hist& h1, const Hist& h2, ostream& os = cout) ;
449 void table(const Hist& h1, const Hist& h2, string fileName) ;
450 
451 // Operator overloading with friends
452 Hist operator+(double f, const Hist& h1);
453 Hist operator+(const Hist& h1, double f);
454 Hist operator+(const Hist& h1, const Hist& h2);
455 Hist operator-(double f, const Hist& h1);
456 Hist operator-(const Hist& h1, double f);
457 Hist operator-(const Hist& h1, const Hist& h2);
458 Hist operator*(double f, const Hist& h1);
459 Hist operator*(const Hist& h1, double f);
460 Hist operator*(const Hist& h1, const Hist& h2);
461 Hist operator/(double f, const Hist& h1);
462 Hist operator/(const Hist& h1, double f);
463 Hist operator/(const Hist& h1, const Hist& h2);
464 
465 //==========================================================================
466 
467 } // end namespace Pythia8
468 
469 #endif // end Pythia8_Basics_H
Generated by   doxygen 1.8.5