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) 2020 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 // Forward reference to RotBstMatrix class; needed in Vec4 class.
23 class RotBstMatrix;
24 
25 //--------------------------------------------------------------------------
26 
27 // Vec4 class.
28 // This class implements four-vectors, in energy-momentum space.
29 // (But can equally well be used to hold space-time four-vectors.)
30 
31 class Vec4 {
32 
33 public:
34 
35  // Constructors.
36  Vec4(double xIn = 0., double yIn = 0., double zIn = 0., double tIn = 0.)
37  : xx(xIn), yy(yIn), zz(zIn), tt(tIn) { }
38  Vec4(const Vec4& v) : xx(v.xx), yy(v.yy), zz(v.zz), tt(v.tt) { }
39  Vec4& operator=(const Vec4& v) { if (this != &v) { xx = v.xx; yy = v.yy;
40  zz = v.zz; tt = v.tt; } return *this; }
41  Vec4& operator=(double value) { xx = value; yy = value; zz = value;
42  tt = value; return *this; }
43 
44  // Member functions for input.
45  void reset() {xx = 0.; yy = 0.; zz = 0.; tt = 0.;}
46  void p(double xIn, double yIn, double zIn, double tIn)
47  {xx = xIn; yy = yIn; zz = zIn; tt = tIn;}
48  void p(Vec4 pIn) {xx = pIn.xx; yy = pIn.yy; zz = pIn.zz; tt = pIn.tt;}
49  void px(double xIn) {xx = xIn;}
50  void py(double yIn) {yy = yIn;}
51  void pz(double zIn) {zz = zIn;}
52  void e(double tIn) {tt = tIn;}
53 
54  // Member functions for output.
55  double px() const {return xx;}
56  double py() const {return yy;}
57  double pz() const {return zz;}
58  double e() const {return tt;}
59  double& operator[](int i) {
60  if (i == 1) return xx;
61  else if (i == 2) return yy;
62  else if (i == 3) return zz;
63  else return tt;
64  }
65  double mCalc() const {double temp = tt*tt - xx*xx - yy*yy - zz*zz;
66  return (temp >= 0.) ? sqrt(temp) : -sqrt(-temp);}
67  double m2Calc() const {return tt*tt - xx*xx - yy*yy - zz*zz;}
68  double pT() const {return sqrt(xx*xx + yy*yy);}
69  double pT2() const {return xx*xx + yy*yy;}
70  double pAbs() const {return sqrt(xx*xx + yy*yy + zz*zz);}
71  double pAbs2() const {return xx*xx + yy*yy + zz*zz;}
72  double eT() const {double temp = xx*xx + yy*yy;
73  return tt * sqrt( temp / (temp + zz*zz) );}
74  double eT2() const {double temp = xx*xx + yy*yy;
75  return tt*tt * temp / (temp + zz*zz);}
76  double theta() const {return atan2(sqrt(xx*xx + yy*yy), zz);}
77  double phi() const {return atan2(yy,xx);}
78  double thetaXZ() const {return atan2(xx,zz);}
79  double pPos() const {return tt + zz;}
80  double pNeg() const {return tt - zz;}
81  double rap() const {return 0.5 * log( (tt + zz) / (tt - zz) );}
82  double eta() const {double xyz = sqrt(xx*xx + yy*yy + zz*zz);
83  return 0.5 * log( (xyz + zz) / (xyz - zz) );}
84 
85  // Member functions that perform operations.
86  void rescale3(double fac) {xx *= fac; yy *= fac; zz *= fac;}
87  void rescale4(double fac) {xx *= fac; yy *= fac; zz *= fac; tt *= fac;}
88  void flip3() {xx = -xx; yy = -yy; zz = -zz;}
89  void flip4() {xx = -xx; yy = -yy; zz = -zz; tt = -tt;}
90  void rot(double thetaIn, double phiIn);
91  void rotaxis(double phiIn, double nx, double ny, double nz);
92  void rotaxis(double phiIn, const Vec4& n);
93  void bst(double betaX, double betaY, double betaZ);
94  void bst(double betaX, double betaY, double betaZ, double gamma);
95  void bst(const Vec4& pIn);
96  void bst(const Vec4& pIn, double mIn);
97  void bstback(const Vec4& pIn);
98  void bstback(const Vec4& pIn, double mIn);
99  void rotbst(const RotBstMatrix& M);
100 
101  // Operator overloading with member functions
102  inline Vec4 operator-() const {Vec4 tmp; tmp.xx = -xx; tmp.yy = -yy;
103  tmp.zz = -zz; tmp.tt = -tt; return tmp;}
104  inline Vec4& operator+=(const Vec4& v) {xx += v.xx; yy += v.yy; zz += v.zz;
105  tt += v.tt; return *this;}
106  inline Vec4& operator-=(const Vec4& v) {xx -= v.xx; yy -= v.yy; zz -= v.zz;
107  tt -= v.tt; return *this;}
108  inline Vec4& operator*=(double f) {xx *= f; yy *= f; zz *= f;
109  tt *= f; return *this;}
110  inline Vec4& operator/=(double f) {xx /= f; yy /= f; zz /= f;
111  tt /= f; return *this;}
112  inline Vec4 operator+(const Vec4& v) const {
113  Vec4 tmp = *this; return tmp += v;}
114  inline Vec4 operator-(const Vec4& v) const {
115  Vec4 tmp = *this; return tmp -= v;}
116  inline Vec4 operator*(double f) const {
117  Vec4 tmp = *this; return tmp *= f;}
118  inline Vec4 operator/(double f) const {
119  Vec4 tmp = *this; return tmp /= f;}
120  inline double operator*(const Vec4& v) const {
121  return tt*v.tt - xx*v.xx - yy*v.yy - zz*v.zz;}
122 
123  // Operator overloading with friends.
124  friend Vec4 operator*(double f, const Vec4& v1);
125 
126  // Print a four-vector.
127  friend ostream& operator<<(ostream&, const Vec4& v) ;
128 
129  // Check if NaN, INF, or finite.
130  friend inline bool isnan(const Vec4 &v) {
131  return isnan(v.tt) || isnan(v.xx) || isnan(v.yy) || isnan(v.zz);}
132  friend inline bool isinf(const Vec4 &v) {
133  return isinf(v.tt) || isinf(v.xx) || isinf(v.yy) || isinf(v.zz);}
134  friend inline bool isfinite(const Vec4 &v) {
135  return isfinite(v.tt) && isfinite(v.xx) && isfinite(v.yy)
136  && isfinite(v.zz);}
137 
138  // Invariant mass of a pair and its square.
139  friend double m(const Vec4& v1, const Vec4& v2);
140  friend double m2(const Vec4& v1, const Vec4& v2);
141 
142  // Scalar and cross product of 3-vector parts.
143  friend double dot3(const Vec4& v1, const Vec4& v2);
144  friend Vec4 cross3(const Vec4& v1, const Vec4& v2);
145 
146  // Cross-product of three 4-vectors ( p_i = epsilon_{iabc} p_a p_b p_c).
147  friend Vec4 cross4(const Vec4& a, const Vec4& b, const Vec4& c);
148 
149  // theta is polar angle between v1 and v2.
150  friend double theta(const Vec4& v1, const Vec4& v2);
151  friend double costheta(const Vec4& v1, const Vec4& v2);
152 
153  // phi is azimuthal angle between v1 and v2 around z axis.
154  friend double phi(const Vec4& v1, const Vec4& v2);
155  friend double cosphi(const Vec4& v1, const Vec4& v2);
156 
157  // phi is azimuthal angle between v1 and v2 around n axis.
158  friend double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
159  friend double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
160 
161  // R is distance in cylindrical (y/eta, phi) coordinates.
162  friend double RRapPhi(const Vec4& v1, const Vec4& v2);
163  friend double REtaPhi(const Vec4& v1, const Vec4& v2);
164 
165  // Shift four-momenta within pair from old to new masses.
166  friend bool pShift( Vec4& p1Move, Vec4& p2Move, double m1New, double m2New);
167 
168  // Create two vectors that are perpendicular to both input vectors.
169  friend pair<Vec4,Vec4> getTwoPerpendicular(const Vec4& v1, const Vec4& v2);
170 
171 private:
172 
173  // Constants: could only be changed in the code itself.
174  static const double TINY;
175 
176  // The four-vector data members.
177  double xx, yy, zz, tt;
178 
179 };
180 
181 //--------------------------------------------------------------------------
182 
183 // Namespace function declarations; friends of Vec4 class.
184 
185 // Implementation of operator overloading with friends.
186 inline Vec4 operator*(double f, const Vec4& v1)
187  {Vec4 v = v1; return v *= f;}
188 
189 // Invariant mass of a pair and its square.
190 double m(const Vec4& v1, const Vec4& v2);
191 double m2(const Vec4& v1, const Vec4& v2);
192 
193 // Scalar and cross product of 3-vector parts.
194 double dot3(const Vec4& v1, const Vec4& v2);
195 Vec4 cross3(const Vec4& v1, const Vec4& v2);
196 
197 // Cross-product of three 4-vectors ( p_i = epsilon_{iabc} p_a p_b p_c).
198 Vec4 cross4(const Vec4& a, const Vec4& b, const Vec4& c);
199 
200 // theta is polar angle between v1 and v2.
201 double theta(const Vec4& v1, const Vec4& v2);
202 double costheta(const Vec4& v1, const Vec4& v2);
203 
204 // phi is azimuthal angle between v1 and v2 around z axis.
205 double phi(const Vec4& v1, const Vec4& v2);
206 double cosphi(const Vec4& v1, const Vec4& v2);
207 
208 // phi is azimuthal angle between v1 and v2 around n axis.
209 double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
210 double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
211 
212 // R is distance in cylindrical (y/eta, phi) coordinates.
213 double RRapPhi(const Vec4& v1, const Vec4& v2);
214 double REtaPhi(const Vec4& v1, const Vec4& v2);
215 
216 // Print a four-vector.
217 ostream& operator<<(ostream&, const Vec4& v) ;
218 
219 // Shift four-momenta within pair from old to new masses.
220 bool pShift( Vec4& p1Move, Vec4& p2Move, double m1New, double m2New);
221 
222 // Create two vectors that are perpendicular to both input vectors.
223 pair<Vec4,Vec4> getTwoPerpendicular(const Vec4& v1, const Vec4& v2);
224 
225 //==========================================================================
226 
227 // RotBstMatrix class.
228 // This class implements 4 * 4 matrices that encode an arbitrary combination
229 // of rotations and boosts, that can be applied to Vec4 four-vectors.
230 
231 class RotBstMatrix {
232 
233 public:
234 
235  // Constructors.
236  RotBstMatrix() : M() {for (int i = 0; i < 4; ++i) {
237  for (int j = 0; j < 4; ++j) { M[i][j] = (i==j) ? 1. : 0.; } } }
238  RotBstMatrix(const RotBstMatrix& Min) : M() {
239  for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
240  M[i][j] = Min.M[i][j]; } } }
241  RotBstMatrix& operator=(const RotBstMatrix& Min) {if (this != &Min) {
242  for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
243  M[i][j] = Min.M[i][j]; } } } return *this; }
244 
245  // Member functions.
246  void rot(double = 0., double = 0.);
247  void rot(const Vec4& p);
248  void bst(double = 0., double = 0., double = 0.);
249  void bst(const Vec4&);
250  void bstback(const Vec4&);
251  void bst(const Vec4&, const Vec4&);
252  void toCMframe(const Vec4&, const Vec4&);
253  void fromCMframe(const Vec4&, const Vec4&);
254  void toSameVframe(const Vec4&, const Vec4&);
255  void fromSameVframe(const Vec4&, const Vec4&);
256  void rotbst(const RotBstMatrix&);
257  void invert();
258  RotBstMatrix inverse() const { RotBstMatrix tmp = *this;
259  tmp.invert(); return tmp; }
260  void reset();
261 
262  // Return value of matrix element.
263  double value(int i, int j) { return M[i][j];}
264 
265  // Crude estimate deviation from unit matrix.
266  double deviation() const;
267 
268  // Print a transformation matrix.
269  friend ostream& operator<<(ostream&, const RotBstMatrix&) ;
270 
271  // Private members to be accessible from Vec4.
272  friend class Vec4;
273 
274  // Multiplication.
275  Vec4 operator*(Vec4 p) const { p.rotbst(*this); return p; }
276  RotBstMatrix operator*(RotBstMatrix R) const { R.rotbst(*this); return R; }
277 
278 private:
279 
280  // Constants: could only be changed in the code itself.
281  static const double TINY;
282 
283  // The rotation-and-boost matrix data members.
284  double M[4][4];
285 
286 };
287 
288 //--------------------------------------------------------------------------
289 
290 // Namespace function declaration; friend of RotBstMatrix class.
291 
292 // Print a transformation matrix.
293 ostream& operator<<(ostream&, const RotBstMatrix&) ;
294 
295 // Get a RotBstMatrix to rest frame of p.
296 inline RotBstMatrix toCMframe(const Vec4& p) {
297  RotBstMatrix tmp; tmp.bstback(p); return tmp; }
298 
299 // Get a RotBstMatrix from rest frame of p.
300 inline RotBstMatrix fromCMframe(const Vec4& p) {
301  RotBstMatrix tmp; tmp.bst(p); return tmp; }
302 
303 // Get a RotBstMatrix to rest frame of p1 and p2, where p1 is along
304 // the z-axis.
305 inline RotBstMatrix toCMframe(const Vec4& p1, const Vec4& p2) {
306  RotBstMatrix tmp; tmp.toCMframe(p1, p2); return tmp; }
307 
308 // Get a RotBstMatrix from rest frame of p1 and p2, where p1 is along
309 // the z-axis.
310 inline RotBstMatrix fromCMframe(const Vec4& p1, const Vec4& p2) {
311  RotBstMatrix tmp; tmp.fromCMframe(p1, p2); return tmp; }
312 
313 // Get a RotBstMatrix to rest frame of ptot where pz is along the
314 // z-axis and pxz is in the xz-plane with positive x.
315 inline RotBstMatrix toCMframe(const Vec4& ptot, const Vec4& pz,
316  const Vec4 & pxz) { RotBstMatrix tmp = toCMframe(ptot);
317  Vec4 pzp = tmp*pz; tmp.rot(0.0, -pzp.phi()); tmp.rot(-pzp.theta());
318  tmp.rot(0.0, -(tmp*pxz).phi()); return tmp; }
319 
320 // Get a RotBstMatrix from rest frame of ptot where pz is along the
321 // z-axis and pxz is in the xz-plane with positive x.
322 inline RotBstMatrix fromCMframe(const Vec4& ptot, const Vec4& pz,
323  const Vec4 & pxz) { return toCMframe(ptot, pz, pxz).inverse(); }
324 
325 //==========================================================================
326 
327 // RndmEngine is the base class for external random number generators.
328 // There is only one pure virtual method, that should do the generation.
329 
330 class RndmEngine {
331 
332 public:
333 
334  // Destructor.
335  virtual ~RndmEngine() {}
336 
337  // A pure virtual method, wherein the derived class method
338  // generates a random number uniformly distributed between 1 and 1.
339  virtual double flat() = 0;
340 
341 };
342 
343 //==========================================================================
344 
345 // Rndm class.
346 // This class handles random number generation according to the
347 // Marsaglia-Zaman-Tsang algorithm.
348 
349 class Rndm {
350 
351 public:
352 
353  // Constructors.
354  Rndm() : initRndm(false), i97(), j97(), seedSave(0), sequence(0), u(), c(),
355  cd(), cm(), useExternalRndm(false), rndmEngPtr(0) { }
356  Rndm(int seedIn) : initRndm(false), i97(), j97(), seedSave(0), sequence(0),
357  u(), c(), cd(), cm(), useExternalRndm(false), rndmEngPtr(0) {init(seedIn);}
358 
359  // Possibility to pass in pointer for external random number generation.
360  bool rndmEnginePtr( RndmEngine* rndmEngPtrIn);
361 
362  // Initialize, normally at construction or in first call.
363  void init(int seedIn = 0) ;
364 
365  // Generate next random number uniformly between 0 and 1.
366  double flat() ;
367 
368  // Generate random numbers according to exp(-x).
369  double exp() { return -log(flat()) ;}
370 
371  // Generate random numbers according to x * exp(-x).
372  double xexp() { return -log(flat() * flat()) ;}
373 
374  // Generate random numbers according to exp(-x^2/2).
375  double gauss() {return sqrt(-2. * log(flat())) * cos(M_PI * flat());}
376 
377  // Generate two random numbers according to exp(-x^2/2-y^2/2).
378  pair<double, double> gauss2() {double r = sqrt(-2. * log(flat()));
379  double phi = 2. * M_PI * flat();
380  return { r * sin(phi), r * cos(phi) };}
381 
382  // Generate two random vectors according to the phase space distribution
383  pair<Vec4, Vec4> phaseSpace2(double eCM, double m1, double m2);
384 
385  // Pick one option among vector of (positive) probabilities.
386  int pick(const vector<double>& prob) ;
387 
388  // Save or read current state to or from a binary file.
389  bool dumpState(string fileName);
390  bool readState(string fileName);
391 
392 private:
393 
394  // Default random number sequence.
395  static const int DEFAULTSEED;
396 
397  // State of the random number generator.
398  bool initRndm;
399  int i97, j97, seedSave;
400  long sequence;
401  double u[97], c, cd, cm;
402 
403  // Pointer for external random number generation.
404  bool useExternalRndm;
405  RndmEngine* rndmEngPtr;
406 
407 };
408 
409 //==========================================================================
410 
411 // Hist class.
412 // This class handles a single histogram at a time.
413 
414 class Hist {
415 
416 public:
417 
418  // Constructors, including copy constructors.
419  Hist() : titleSave(""), nBin(), nFill(), nNonFinite(),
420  xMin(), xMax(), linX(), dx(), under(), inside(), over(), sumxw()
421  { }
422  Hist(string titleIn, int nBinIn = 100, double xMinIn = 0.,
423  double xMaxIn = 1., bool logXIn = false) : nBin(), nFill(), nNonFinite(),
424  xMin(), xMax(), linX(), dx(), under(), inside(), over(), sumxw()
425  { book(titleIn, nBinIn, xMinIn, xMaxIn, logXIn); }
426  Hist(const Hist& h)
427  : titleSave(h.titleSave), nBin(h.nBin), nFill(h.nFill),
428  nNonFinite(h.nNonFinite), xMin(h.xMin), xMax(h.xMax), linX(h.linX),
429  dx(h.dx), under(h.under), inside(h.inside), over(h.over),
430  sumxw(h.sumxw), res(h.res) { }
431  Hist(string titleIn, const Hist& h)
432  : titleSave(titleIn), nBin(h.nBin), nFill(h.nFill),
433  nNonFinite(h.nNonFinite), xMin(h.xMin), xMax(h.xMax), linX(h.linX),
434  dx(h.dx), under(h.under), inside(h.inside), over(h.over),
435  sumxw(h.sumxw), res(h.res) { }
436  Hist& operator=(const Hist& h) { if(this != &h) {
437  nBin = h.nBin; nFill = h.nFill; nNonFinite = h.nNonFinite; xMin = h.xMin;
438  xMax = h.xMax; linX = h.linX; dx = h.dx; under = h.under;
439  inside = h.inside; over = h.over; sumxw = h.sumxw;
440  res = h.res; } return *this; }
441 
442  // Create a histogram that is the plot of the given function.
443  static Hist plotFunc(function<double(double)> f, string titleIn,
444  int nBinIn, double xMinIn, double xMaxIn, bool logXIn = false);
445 
446  // Book a histogram.
447  void book(string titleIn = " ", int nBinIn = 100, double xMinIn = 0.,
448  double xMaxIn = 1., bool logXIn = false) ;
449 
450  // Set title of a histogram.
451  void title(string titleIn = " ") {titleSave = titleIn; }
452 
453  // Reset bin contents.
454  void null() ;
455 
456  // Fill bin with weight.
457  void fill(double x, double w = 1.) ;
458 
459  // Print a histogram with overloaded << operator.
460  friend ostream& operator<<(ostream& os, const Hist& h) ;
461 
462  // Print histogram contents as a table (e.g. for Gnuplot, Rivet or Pyplot).
463  void table(ostream& os = cout, bool printOverUnder = false,
464  bool xMidBin = true) const ;
465  void table(string fileName, bool printOverUnder = false,
466  bool xMidBin = true) const { ofstream streamName(fileName.c_str());
467  table(streamName, printOverUnder, xMidBin);}
468  void rivetTable(ostream& os = cout, bool printError = false) const ;
469  void rivetTable(string fileName, bool printError = false) const {
470  ofstream streamName(fileName.c_str()); rivetTable(streamName, printError);}
471  void pyplotTable(ostream& os = cout, bool isHist = true) const ;
472  void pyplotTable(string fileName, bool isHist = true) const {
473  ofstream streamName(fileName.c_str()); pyplotTable(streamName, isHist);}
474 
475  // Print a table out of two histograms with same x axis.
476  friend void table(const Hist& h1, const Hist& h2, ostream& os,
477  bool printOverUnder, bool xMidBin) ;
478  friend void table(const Hist& h1, const Hist& h2, string fileName,
479  bool printOverUnder, bool xMidBin) ;
480 
481  // Return title and size of histogram. Also if logarithmic x scale.
482  string getTitle() const {return titleSave;}
483  int getBinNumber() const {return nBin;}
484  int getNonFinite() const {return nNonFinite;}
485  bool getLinX() const {return linX;}
486 
487  // Return min and max in x and y directions.
488  double getXMin() const {return xMin;}
489  double getXMax() const {return xMax;}
490  double getYMin() const { double yMin = res[0];
491  for (int ix = 1; ix < nBin; ++ix)
492  if (res[ix] < yMin ) yMin = res[ix];
493  return yMin;}
494  double getYMax() const {double yMax = res[0];
495  for (int ix = 1; ix < nBin; ++ix)
496  if (res[ix] > yMax ) yMax = res[ix];
497  return yMax;}
498  double getYAbsMin() const { double yAbsMin = 1e20; double yAbs;
499  for (int ix = 0; ix < nBin; ++ix) { yAbs = abs(res[ix]);
500  if (yAbs > 1e-20 && yAbs < yAbsMin) yAbsMin = yAbs; }
501  return yAbsMin;}
502 
503  double getXMean() const { return sumxw / inside; }
504  double getYMean() const { return inside / nFill; }
505 
506  // Return content of specific bin: 0 gives underflow and nBin+1 overflow.
507  double getBinContent(int iBin) const;
508 
509  // Return number of entries.
510  int getEntries(bool alsoNonFinite = true) const {
511  return alsoNonFinite ? nNonFinite + nFill : nFill; }
512 
513  // Check whether another histogram has same size and limits.
514  bool sameSize(const Hist& h) const ;
515 
516  // Take logarithm (base 10 or e) of bin contents.
517  void takeLog(bool tenLog = true) ;
518 
519  // Take square root of bin contents.
520  void takeSqrt() ;
521 
522  // Normalize bin contents to given sum, by default including overflow bins.
523  void normalize( double sum = 1., bool alsoOverflow = true) ;
524 
525  // Operator overloading with member functions
526  Hist& operator+=(const Hist& h) ;
527  Hist& operator-=(const Hist& h) ;
528  Hist& operator*=(const Hist& h) ;
529  Hist& operator/=(const Hist& h) ;
530  Hist& operator+=(double f) ;
531  Hist& operator-=(double f) ;
532  Hist& operator*=(double f) ;
533  Hist& operator/=(double f) ;
534  Hist operator+(double f) const;
535  Hist operator+(const Hist& h2) const;
536  Hist operator-(double f) const;
537  Hist operator-(const Hist& h2) const;
538  Hist operator*(double f) const;
539  Hist operator*(const Hist& h2) const;
540  Hist operator/(double f) const;
541  Hist operator/(const Hist& h2) const;
542 
543  // Operator overloading with friends
544  friend Hist operator+(double f, const Hist& h1);
545  friend Hist operator-(double f, const Hist& h1);
546  friend Hist operator*(double f, const Hist& h1);
547  friend Hist operator/(double f, const Hist& h1);
548 
549 private:
550 
551  // Constants: could only be changed in the code itself.
552  static const int NBINMAX, NCOLMAX, NLINES;
553  static const double TOLERANCE, TINY, LARGE, SMALLFRAC, DYAC[];
554  static const char NUMBER[];
555 
556  // Properties and contents of a histogram.
557  string titleSave;
558  int nBin, nFill, nNonFinite;
559  double xMin, xMax;
560  bool linX;
561  double dx, under, inside, over, sumxw;
562  vector<double> res;
563 
564 };
565 
566 //--------------------------------------------------------------------------
567 
568 // Namespace function declarations; friends of Hist class.
569 
570 // Print a histogram with overloaded << operator.
571 ostream& operator<<(ostream& os, const Hist& h) ;
572 
573 // Print a table out of two histograms with same x axis.
574 void table(const Hist& h1, const Hist& h2, ostream& os = cout,
575  bool printOverUnder = false, bool xMidBin = true) ;
576 void table(const Hist& h1, const Hist& h2, string fileName,
577  bool printOverUnder = false, bool xMidBin = true) ;
578 
579 // Operator overloading with friends
580 Hist operator+(double f, const Hist& h1);
581 Hist operator-(double f, const Hist& h1);
582 Hist operator*(double f, const Hist& h1);
583 Hist operator/(double f, const Hist& h1);
584 
585 //==========================================================================
586 
587 // HistPlot class.
588 // Writes a Python program that can generate PDF plots from Hist histograms.
589 
590 class HistPlot {
591 
592 public:
593 
594  // Constructor requires name of Python program (and adds .py).
595  HistPlot(string pythonName) : nFrame(), nTable() {
596  toPython.open( (pythonName + ".py").c_str() );
597  toPython << "from matplotlib import pyplot as plt" << endl
598  << "from matplotlib.backends.backend_pdf import PdfPages" << endl;
599  nPDF = 0; }
600 
601  // Destructor should do final close.
602  ~HistPlot() { toPython << "pp.close()" << endl; }
603 
604  // New plot frame, with title, x and y labels, x and y sizes..
605  void frame( string frameIn, string titleIn = "", string xLabIn = "",
606  string yLabIn = "", double xSizeIn = 8., double ySizeIn = 6.) {
607  framePrevious = frameName; frameName = frameIn; title = titleIn;
608  xLabel = xLabIn; yLabel = yLabIn; xSize = xSizeIn; ySize = ySizeIn;
609  histos.clear(); styles.clear(); legends.clear(); files.clear();
610  fileStyles.clear(); fileLegends.clear(); filexyerr.clear();}
611 
612  // Add a histogram to the current plot, with optional style and legend.
613  void add( const Hist& histIn, string styleIn = "h",
614  string legendIn = "void") { histos.push_back(histIn);
615  styles.push_back(styleIn); legends.push_back(legendIn); }
616 
617  // Add a file of (x, y) values not from a histogram, e.g. data points.
618  void addFile( string fileIn, string styleIn = "o",
619  string legendIn = "void", string xyerrIn="") { files.push_back(fileIn);
620  fileStyles.push_back(styleIn); fileLegends.push_back(legendIn);
621  filexyerr.push_back(xyerrIn);}
622 
623  // Plot a frame given the information from the new and add calls.
624  void plot( bool logY = false, bool logX = false, bool userBorders = false);
625  void plot( double xMinUserIn, double xMaxUserIn, double yMinUserIn,
626  double yMaxUserIn, bool logY = false, bool logX = false) {
627  xMinUser = xMinUserIn; xMaxUser = xMaxUserIn; yMinUser = yMinUserIn;
628  yMaxUser = yMaxUserIn; plot( logY, logX, true);}
629 
630  // Omnibus single call when only one histogram in the frame.
631  void plotFrame( string frameIn, const Hist& histIn, string titleIn = "",
632  string xLabIn = "", string yLabIn = "", string styleIn = "h",
633  string legendIn = "void", bool logY = false) {
634  frame( frameIn, titleIn, xLabIn, yLabIn);
635  add( histIn, styleIn, legendIn); plot( logY); }
636 
637 private:
638 
639  // Initialization code.
640  void init( string pythonName);
641 
642  // Stored quantities.
643  ofstream toPython;
644  int nPDF, nFrame, nTable;
645  double xSize, ySize, xMinUser, xMaxUser, yMinUser, yMaxUser;
646  string frameName, framePrevious, title, xLabel, yLabel, fileName, tmpFig;
647  vector<Hist> histos;
648  vector<string> styles, legends, files, fileStyles, fileLegends, filexyerr;
649 
650 };
651 
652 //==========================================================================
653 
654 } // end namespace Pythia8
655 
656 #endif // end Pythia8_Basics_H
Generated by   doxygen 1.8.5