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) 2018 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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);
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  : 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& operator[](int i) {
141  if (i == 1) return xx;
142  else if (i == 2) return yy;
143  else if (i == 3) return zz;
144  else return tt;
145  }
146  double mCalc() const {double temp = tt*tt - xx*xx - yy*yy - zz*zz;
147  return (temp >= 0.) ? sqrt(temp) : -sqrt(-temp);}
148  double m2Calc() const {return tt*tt - xx*xx - yy*yy - zz*zz;}
149  double pT() const {return sqrt(xx*xx + yy*yy);}
150  double pT2() const {return xx*xx + yy*yy;}
151  double pAbs() const {return sqrt(xx*xx + yy*yy + zz*zz);}
152  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  double theta() const {return atan2(sqrt(xx*xx + yy*yy), zz);}
158  double phi() const {return atan2(yy,xx);}
159  double thetaXZ() const {return atan2(xx,zz);}
160  double pPos() const {return tt + zz;}
161  double pNeg() const {return tt - zz;}
162  double rap() const {return 0.5 * log( (tt + zz) / (tt - zz) );}
163  double eta() const {double xyz = sqrt(xx*xx + yy*yy + zz*zz);
164  return 0.5 * log( (xyz + zz) / (xyz - zz) );}
165 
166  // Member functions that perform operations.
167  void rescale3(double fac) {xx *= fac; yy *= fac; zz *= fac;}
168  void rescale4(double fac) {xx *= fac; yy *= fac; zz *= fac; tt *= fac;}
169  void flip3() {xx = -xx; yy = -yy; zz = -zz;}
170  void flip4() {xx = -xx; yy = -yy; zz = -zz; tt = -tt;}
171  void rot(double thetaIn, double phiIn);
172  void rotaxis(double phiIn, double nx, double ny, double nz);
173  void rotaxis(double phiIn, const Vec4& n);
174  void bst(double betaX, double betaY, double betaZ);
175  void bst(double betaX, double betaY, double betaZ, double gamma);
176  void bst(const Vec4& pIn);
177  void bst(const Vec4& pIn, double mIn);
178  void bstback(const Vec4& pIn);
179  void bstback(const Vec4& pIn, double mIn);
180  void rotbst(const RotBstMatrix& M);
181 
182  // Operator overloading with member functions
183  inline Vec4 operator-() const {Vec4 tmp; tmp.xx = -xx; tmp.yy = -yy;
184  tmp.zz = -zz; tmp.tt = -tt; return tmp;}
185  inline Vec4& operator+=(const Vec4& v) {xx += v.xx; yy += v.yy; zz += v.zz;
186  tt += v.tt; return *this;}
187  inline Vec4& operator-=(const Vec4& v) {xx -= v.xx; yy -= v.yy; zz -= v.zz;
188  tt -= v.tt; return *this;}
189  inline Vec4& operator*=(double f) {xx *= f; yy *= f; zz *= f;
190  tt *= f; return *this;}
191  inline Vec4& operator/=(double f) {xx /= f; yy /= f; zz /= f;
192  tt /= f; return *this;}
193  inline Vec4 operator+(const Vec4& v) const {
194  Vec4 tmp = *this; return tmp += v;}
195  inline Vec4 operator-(const Vec4& v) const {
196  Vec4 tmp = *this; return tmp -= v;}
197  inline Vec4 operator*(double f) const {
198  Vec4 tmp = *this; return tmp *= f;}
199  inline Vec4 operator/(double f) const {
200  Vec4 tmp = *this; return tmp /= f;}
201  inline double operator*(const Vec4& v) const {
202  return tt*v.tt - xx*v.xx - yy*v.yy - zz*v.zz;}
203 
204  // Operator overloading with friends
205  friend Vec4 operator*(double f, const Vec4& v1);
206 
207  // Print a four-vector.
208  friend ostream& operator<<(ostream&, const Vec4& v) ;
209 
210  // Invariant mass of a pair and its square.
211  friend double m(const Vec4& v1, const Vec4& v2);
212  friend double m2(const Vec4& v1, const Vec4& v2);
213 
214  // Scalar and cross product of 3-vector parts.
215  friend double dot3(const Vec4& v1, const Vec4& v2);
216  friend Vec4 cross3(const Vec4& v1, const Vec4& v2);
217 
218  // Cross-product of three 4-vectors ( p_i = epsilon_{iabc} p_a p_b p_c).
219  friend Vec4 cross4(const Vec4& a, const Vec4& b, const Vec4& c);
220 
221  // theta is polar angle between v1 and v2.
222  friend double theta(const Vec4& v1, const Vec4& v2);
223  friend double costheta(const Vec4& v1, const Vec4& v2);
224 
225  // phi is azimuthal angle between v1 and v2 around z axis.
226  friend double phi(const Vec4& v1, const Vec4& v2);
227  friend double cosphi(const Vec4& v1, const Vec4& v2);
228 
229  // phi is azimuthal angle between v1 and v2 around n axis.
230  friend double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
231  friend double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
232 
233  // R is distance in cylindrical (y/eta, phi) coordinates.
234  friend double RRapPhi(const Vec4& v1, const Vec4& v2);
235  friend double REtaPhi(const Vec4& v1, const Vec4& v2);
236 
237  // Shift four-momenta within pair from old to new masses.
238  friend bool pShift( Vec4& p1Move, Vec4& p2Move, double m1New, double m2New);
239 
240  // Create two vectors that are perpendicular to both input vectors.
241  friend pair<Vec4,Vec4> getTwoPerpendicular(const Vec4& v1, const Vec4& v2);
242 
243 private:
244 
245  // Constants: could only be changed in the code itself.
246  static const double TINY;
247 
248  // The four-vector data members.
249  double xx, yy, zz, tt;
250 
251 };
252 
253 //--------------------------------------------------------------------------
254 
255 // Namespace function declarations; friends of Vec4 class.
256 
257 // Implementation of operator overloading with friends.
258 inline Vec4 operator*(double f, const Vec4& v1)
259  {Vec4 v = v1; return v *= f;}
260 
261 // Invariant mass of a pair and its square.
262 double m(const Vec4& v1, const Vec4& v2);
263 double m2(const Vec4& v1, const Vec4& v2);
264 
265 // Scalar and cross product of 3-vector parts.
266 double dot3(const Vec4& v1, const Vec4& v2);
267 Vec4 cross3(const Vec4& v1, const Vec4& v2);
268 
269 // Cross-product of three 4-vectors ( p_i = epsilon_{iabc} p_a p_b p_c).
270 Vec4 cross4(const Vec4& a, const Vec4& b, const Vec4& c);
271 
272 // theta is polar angle between v1 and v2.
273 double theta(const Vec4& v1, const Vec4& v2);
274 double costheta(const Vec4& v1, const Vec4& v2);
275 
276 // phi is azimuthal angle between v1 and v2 around z axis.
277 double phi(const Vec4& v1, const Vec4& v2);
278 double cosphi(const Vec4& v1, const Vec4& v2);
279 
280 // phi is azimuthal angle between v1 and v2 around n axis.
281 double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
282 double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
283 
284 // R is distance in cylindrical (y/eta, phi) coordinates.
285 double RRapPhi(const Vec4& v1, const Vec4& v2);
286 double REtaPhi(const Vec4& v1, const Vec4& v2);
287 
288 // Print a four-vector.
289 ostream& operator<<(ostream&, const Vec4& v) ;
290 
291 // Shift four-momenta within pair from old to new masses.
292 bool pShift( Vec4& p1Move, Vec4& p2Move, double m1New, double m2New);
293 
294 // Create two vectors that are perpendicular to both input vectors.
295 pair<Vec4,Vec4> getTwoPerpendicular(const Vec4& v1, const Vec4& v2);
296 
297 //==========================================================================
298 
299 // RotBstMatrix class.
300 // This class implements 4 * 4 matrices that encode an arbitrary combination
301 // of rotations and boosts, that can be applied to Vec4 four-vectors.
302 
303 class RotBstMatrix {
304 
305 public:
306 
307  // Constructors.
308  RotBstMatrix() {for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j)
309  { M[i][j] = (i==j) ? 1. : 0.; } } }
310  RotBstMatrix(const RotBstMatrix& Min) {
311  for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
312  M[i][j] = Min.M[i][j]; } } }
313  RotBstMatrix& operator=(const RotBstMatrix& Min) {if (this != &Min) {
314  for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
315  M[i][j] = Min.M[i][j]; } } } return *this; }
316 
317  // Member functions.
318  void rot(double = 0., double = 0.);
319  void rot(const Vec4& p);
320  void bst(double = 0., double = 0., double = 0.);
321  void bst(const Vec4&);
322  void bstback(const Vec4&);
323  void bst(const Vec4&, const Vec4&);
324  void toCMframe(const Vec4&, const Vec4&);
325  void fromCMframe(const Vec4&, const Vec4&);
326  void rotbst(const RotBstMatrix&);
327  void invert();
328  void reset();
329 
330  // Return value of matrix element.
331  double value(int i, int j) { return M[i][j];}
332 
333  // Crude estimate deviation from unit matrix.
334  double deviation() const;
335 
336  // Print a transformation matrix.
337  friend ostream& operator<<(ostream&, const RotBstMatrix&) ;
338 
339  // Private members to be accessible from Vec4.
340  friend class Vec4;
341 
342 private:
343 
344  // Constants: could only be changed in the code itself.
345  static const double TINY;
346 
347  // The rotation-and-boost matrix data members.
348  double M[4][4];
349 
350 };
351 
352 //--------------------------------------------------------------------------
353 
354 // Namespace function declaration; friend of RotBstMatrix class.
355 
356 // Print a transformation matrix.
357 ostream& operator<<(ostream&, const RotBstMatrix&) ;
358 
359 //==========================================================================
360 
361 // Hist class.
362 // This class handles a single histogram at a time.
363 
364 class Hist {
365 
366 public:
367 
368  // Constructors, including copy constructors.
369  Hist() {titleSave = "";}
370  Hist(string titleIn, int nBinIn = 100, double xMinIn = 0.,
371  double xMaxIn = 1., bool logXIn = false) {
372  book(titleIn, nBinIn, xMinIn, xMaxIn, logXIn);}
373  Hist(const Hist& h)
374  : titleSave(h.titleSave), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
375  xMax(h.xMax), linX(h.linX), dx(h.dx), under(h.under), inside(h.inside),
376  over(h.over), res(h.res) { }
377  Hist(string titleIn, const Hist& h)
378  : titleSave(titleIn), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
379  xMax(h.xMax), linX(h.linX), dx(h.dx), under(h.under), inside(h.inside),
380  over(h.over), res(h.res) { }
381  Hist& operator=(const Hist& h) { if(this != &h) {
382  nBin = h.nBin; nFill = h.nFill; xMin = h.xMin; xMax = h.xMax;
383  linX = h.linX; dx = h.dx; under = h.under; inside = h.inside;
384  over = h.over; res = h.res; } return *this; }
385 
386  // Book a histogram.
387  void book(string titleIn = " ", int nBinIn = 100, double xMinIn = 0.,
388  double xMaxIn = 1., bool logXIn = false) ;
389 
390  // Set title of a histogram.
391  void title(string titleIn = " ") {titleSave = titleIn; }
392 
393  // Reset bin contents.
394  void null() ;
395 
396  // Fill bin with weight.
397  void fill(double x, double w = 1.) ;
398 
399  // Print a histogram with overloaded << operator.
400  friend ostream& operator<<(ostream& os, const Hist& h) ;
401 
402  // Print histogram contents as a table (e.g. for Gnuplot, Rivet or Pyplot).
403  void table(ostream& os = cout, bool printOverUnder = false,
404  bool xMidBin = true) const ;
405  void table(string fileName, bool printOverUnder = false,
406  bool xMidBin = true) const { ofstream streamName(fileName.c_str());
407  table(streamName, printOverUnder, xMidBin);}
408  void rivetTable(ostream& os = cout, bool printError = false) const ;
409  void rivetTable(string fileName, bool printError = false) const {
410  ofstream streamName(fileName.c_str()); rivetTable(streamName, printError);}
411  void pyplotTable(ostream& os = cout, bool isHist = true) const ;
412  void pyplotTable(string fileName, bool isHist = true) const {
413  ofstream streamName(fileName.c_str()); pyplotTable(streamName, isHist);}
414 
415  // Print a table out of two histograms with same x axis.
416  friend void table(const Hist& h1, const Hist& h2, ostream& os,
417  bool printOverUnder, bool xMidBin) ;
418  friend void table(const Hist& h1, const Hist& h2, string fileName,
419  bool printOverUnder, bool xMidBin) ;
420 
421  // Return title and size of histogram. Also if logarithmic x scale.
422  string getTitle() const {return titleSave;}
423  int getBinNumber() const {return nBin;}
424  bool getLinX() const {return linX;}
425 
426  // Return content of specific bin: 0 gives underflow and nBin+1 overflow.
427  double getBinContent(int iBin) const;
428 
429  // Return number of entries.
430  int getEntries() const {return nFill; }
431 
432  // Check whether another histogram has same size and limits.
433  bool sameSize(const Hist& h) const ;
434 
435  // Take logarithm (base 10 or e) of bin contents.
436  void takeLog(bool tenLog = true) ;
437 
438  // Take square root of bin contents.
439  void takeSqrt() ;
440 
441  // Find smallest nonzero absolute value of bin contents.
442  double smallestAbsValue() const ;
443 
444  // Operator overloading with member functions
445  Hist& operator+=(const Hist& h) ;
446  Hist& operator-=(const Hist& h) ;
447  Hist& operator*=(const Hist& h) ;
448  Hist& operator/=(const Hist& h) ;
449  Hist& operator+=(double f) ;
450  Hist& operator-=(double f) ;
451  Hist& operator*=(double f) ;
452  Hist& operator/=(double f) ;
453  Hist operator+(double f) const;
454  Hist operator+(const Hist& h2) const;
455  Hist operator-(double f) const;
456  Hist operator-(const Hist& h2) const;
457  Hist operator*(double f) const;
458  Hist operator*(const Hist& h2) const;
459  Hist operator/(double f) const;
460  Hist operator/(const Hist& h2) const;
461 
462  // Operator overloading with friends
463  friend Hist operator+(double f, const Hist& h1);
464  friend Hist operator-(double f, const Hist& h1);
465  friend Hist operator*(double f, const Hist& h1);
466  friend Hist operator/(double f, const Hist& h1);
467 
468 private:
469 
470  // Constants: could only be changed in the code itself.
471  static const int NBINMAX, NCOLMAX, NLINES;
472  static const double TOLERANCE, TINY, LARGE, SMALLFRAC, DYAC[];
473  static const char NUMBER[];
474 
475  // Properties and contents of a histogram.
476  string titleSave;
477  int nBin, nFill;
478  double xMin, xMax;
479  bool linX;
480  double dx, under, inside, over;
481  vector<double> res;
482 
483 };
484 
485 //--------------------------------------------------------------------------
486 
487 // Namespace function declarations; friends of Hist class.
488 
489 // Print a histogram with overloaded << operator.
490 ostream& operator<<(ostream& os, const Hist& h) ;
491 
492 // Print a table out of two histograms with same x axis.
493 void table(const Hist& h1, const Hist& h2, ostream& os = cout,
494  bool printOverUnder = false, bool xMidBin = true) ;
495 void table(const Hist& h1, const Hist& h2, string fileName,
496  bool printOverUnder = false, bool xMidBin = true) ;
497 
498 // Operator overloading with friends
499 Hist operator+(double f, const Hist& h1);
500 Hist operator-(double f, const Hist& h1);
501 Hist operator*(double f, const Hist& h1);
502 Hist operator/(double f, const Hist& h1);
503 
504 //==========================================================================
505 
506 // HistPlot class.
507 // Writes a Python program that can generate PDF plots from Hist histograms.
508 
509 class HistPlot {
510 
511 public:
512 
513  // Constructor requires name of Python program (and adds .py).
514  HistPlot(string pythonName) { toPython.open( (pythonName + ".py").c_str() );
515  toPython << "from matplotlib import pyplot as plt" << endl
516  << "from matplotlib.backends.backend_pdf import PdfPages" << endl;
517  nPDF = 0; }
518 
519  // Destructor should do final close.
520  ~HistPlot() { toPython << "pp.close()" << endl; }
521 
522  // New plot frame, with title, x and y labels.
523  void frame( string frameIn, string titleIn = "", string xLabIn = "",
524  string yLabIn = "") {frameName = frameIn; title = titleIn; xLabel = xLabIn;
525  yLabel = yLabIn; histos.clear(); styles.clear(); legends.clear(); }
526 
527  // Add a histogram to the current plot, with optional style and legend.
528  void add( const Hist& histIn, string styleIn = "h",
529  string legendIn = "void") { histos.push_back(&histIn);
530  styles.push_back(styleIn); legends.push_back(legendIn); }
531 
532  // Plot a frame given the information from the new and add calls.
533  void plot( bool logY = false);
534 
535  // Omnibus single call when only one histogram in the frame.
536  void plotFrame( string frameIn, const Hist& histIn, string titleIn = "",
537  string xLabIn = "", string yLabIn = "", string styleIn = "h",
538  string legendIn = "void", bool logY = false) {
539  frame( frameIn, titleIn, xLabIn, yLabIn);
540  add( histIn, styleIn, legendIn); plot( logY); }
541 
542 private:
543 
544  // Initialization code.
545  void init( string pythonName);
546 
547  // Stored quantities.
548  ofstream toPython;
549  int nPDF, nFrame, nTable;
550  string frameName, title, xLabel, yLabel, fileName, tmpFig;
551  vector<const Hist*> histos;
552  vector<string> styles, legends;
553 
554 };
555 
556 //==========================================================================
557 
558 } // end namespace Pythia8
559 
560 #endif // end Pythia8_Basics_H
Generated by   doxygen 1.8.5