PYTHIA8/pythia8170/include/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