PYTHIA8/pythia8130/include/Basics.h

   1 // Basics.h is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 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 (static member functions).
   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() {}
  52   Rndm(int seedIn) { init(seedIn);}
  53
  54   // Possibility to pass in pointer for external random number generation.
  55   static bool rndmEnginePtr( RndmEngine* rndmPtrIn);
  56
  57   // Initialize, normally at construction or in first call.
  58   static void init(int seedIn = 0) ;
  59
  60   // Generate next random number uniformly between 0 and 1.
  61   static double flat() ;
  62
  63   // Generate random numbers according to exp(-x).
  64   static double exp() { return -log(flat()) ;}
  65
  66   // Generate random numbers according to x * exp(-x).
  67   static double xexp() { return -log(flat() * flat()) ;}
  68
  69   // Generate random numbers according to exp(-x^2/2).
  70   static double gauss() ;
  71
  72   // Pick one option among  vector of (positive) probabilities.
  73   static int pick(const vector<double>& prob) ;
  74
  75 private:
  76
  77   // State of the random number generator.
  78   static bool   initRndm, saveGauss;
  79   static int    i97, j97, defaultSeed;
  80   static double u[97], c, cd, cm, save;
  81
  82   // Pointer for external random number generation.
  83   static bool   useExternalRndm;
  84   static RndmEngine* rndmPtr;
  85
  86 };
  87
  88 //**************************************************************************
  89
  90 // Forward reference to RotBstMatrix class.
  91 class RotBstMatrix;
  92
  93 //**************************************************************************
  94
  95 // Vec4 class.
  96 // This class implements four-vectors, in energy-momentum space.
  97 // (But can equally well be used to hold space-time four-vectors.)
  98
  99 class Vec4 {
 100
 101 public:
 102
 103   // Constructors.
 104   Vec4(double xIn = 0., double yIn = 0., double zIn = 0., double tIn = 0.)
 105     : xx(xIn), yy(yIn), zz(zIn), tt(tIn) { }
 106   Vec4(const Vec4& v) : xx(v.xx), yy(v.yy), zz(v.zz), tt(v.tt) { }
 107   Vec4& operator=(const Vec4& v) { if (this != &v) { xx = v.xx; yy = v.yy;
 108     zz = v.zz; tt = v.tt; } return *this; }
 109   Vec4& operator=(double value) { xx = value; yy = value; zz = value;
 110     tt = value; return *this; }
 111
 112   // Member functions for input.
 113   void reset() {xx = 0.; yy = 0.; zz = 0.; tt = 0.;}
 114   void p(double xIn, double yIn, double zIn, double tIn)
 115     {xx = xIn; yy = yIn; zz = zIn; tt = tIn;}
 116   void p(Vec4 pIn) {xx = pIn.xx; yy = pIn.yy; zz = pIn.zz; tt = pIn.tt;}
 117   void px(double xIn) {xx = xIn;}
 118   void py(double yIn) {yy = yIn;}
 119   void pz(double zIn) {zz = zIn;}
 120   void e(double tIn) {tt = tIn;}
 121
 122   // Member functions for output.
 123   double px() const {return xx;}
 124   double py() const {return yy;}
 125   double pz() const {return zz;}
 126   double e() const {return tt;}
 127   double mCalc() const {double temp = tt*tt - xx*xx - yy*yy - zz*zz;
 128     return (temp >= 0) ? sqrt(temp) : -sqrt(-temp);}
 129   double m2Calc() const {return tt*tt - xx*xx - yy*yy - zz*zz;}
 130   double pT() const {return sqrt(xx*xx + yy*yy);}
 131   double pT2() const {return xx*xx + yy*yy;}
 132   double pAbs() const {return sqrt(xx*xx + yy*yy + zz*zz);}
 133   double pAbs2() const {return xx*xx + yy*yy + zz*zz;}
 134   double theta() const {return atan2(sqrt(xx*xx + yy*yy), zz);}
 135   double phi() const {return atan2(yy,xx);}
 136   double thetaXZ() const {return atan2(xx,zz);}
 137   double pPlus() const {return tt + zz;}
 138   double pMinus() const {return tt - zz;}
 139
 140   // Member functions that perform operations.
 141   void rescale3(double fac) {xx *= fac; yy *= fac; zz *= fac;}
 142   void rescale4(double fac) {xx *= fac; yy *= fac; zz *= fac; tt *= fac;}
 143   void flip3() {xx = -xx; yy = -yy; zz = -zz;}
 144   void flip4() {xx = -xx; yy = -yy; zz = -zz; tt = -tt;}
 145   void rot(double thetaIn, double phiIn);
 146   void rotaxis(double phiIn, double nx, double ny, double nz);
 147   void rotaxis(double phiIn, const Vec4& n);
 148   void bst(double betaX, double betaY, double betaZ);
 149   void bst(double betaX, double betaY, double betaZ, double gamma);
 150   void bst(const Vec4& pIn);
 151   void bst(const Vec4& pIn, double mIn);
 152   void bstback(const Vec4& pIn);
 153   void bstback(const Vec4& pIn, double mIn);
 154   void rotbst(const RotBstMatrix& M);
 155
 156   // Operator overloading with member functions
 157   Vec4& operator-() {xx = -xx; yy = -yy; zz = -zz; tt = -tt; return *this;}
 158   Vec4& operator+=(const Vec4& v) {xx += v.xx; yy += v.yy; zz += v.zz;
 159     tt += v.tt; return *this;}
 160   Vec4& operator-=(const Vec4& v) {xx -= v.xx; yy -= v.yy; zz -= v.zz;
 161     tt -= v.tt; return *this;}
 162   Vec4& operator*=(double f) {xx *= f; yy *= f; zz *= f;
 163     tt *= f; return *this;}
 164   Vec4& operator/=(double f) {xx /= f; yy /= f; zz /= f;
 165     tt /= f; return *this;}
 166
 167   // Operator overloading with friends
 168   friend Vec4 operator+(const Vec4& v1, const Vec4& v2);
 169   friend Vec4 operator-(const Vec4& v1, const Vec4& v2);
 170   friend Vec4 operator*(double f, const Vec4& v1);
 171   friend Vec4 operator*(const Vec4& v1, double f);
 172   friend Vec4 operator/(const Vec4& v1, double f);
 173   friend double operator*(const Vec4& v1, const Vec4& v2);
 174
 175   // Invariant mass of a pair and its square.
 176   friend double m(const Vec4& v1, const Vec4& v2);
 177   friend double m2(const Vec4& v1, const Vec4& v2);
 178
 179   // Scalar and cross product of 3-vector parts.
 180   friend double dot3(const Vec4& v1, const Vec4& v2);
 181   friend Vec4 cross3(const Vec4& v1, const Vec4& v2);
 182
 183   // theta is polar angle between v1 and v2.
 184   friend double theta(const Vec4& v1, const Vec4& v2);
 185   friend double costheta(const Vec4& v1, const Vec4& v2);
 186
 187   // phi is azimuthal angle between v1 and v2 around z axis.
 188   friend double phi(const Vec4& v1, const Vec4& v2);
 189   friend double cosphi(const Vec4& v1, const Vec4& v2);
 190
 191   // phi is azimuthal angle between v1 and v2 around n axis.
 192   friend double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
 193   friend double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
 194
 195   // Print a four-vector
 196   friend ostream& operator<<(ostream&, const Vec4& v) ;
 197
 198 private:
 199
 200   // Constants: could only be changed in the code itself.
 201   static const double TINY;
 202
 203   // The four-vector data members.
 204   double xx, yy, zz, tt;
 205
 206 };
 207
 208 // Implementation of operator overloading with friends.
 209
 210 inline Vec4 operator+(const Vec4& v1, const Vec4& v2)
 211   {Vec4 v = v1 ; return v += v2;}
 212
 213 inline Vec4 operator-(const Vec4& v1, const Vec4& v2)
 214   {Vec4 v = v1 ; return v -= v2;}
 215
 216 inline Vec4 operator*(double f, const Vec4& v1)
 217   {Vec4 v = v1; return v *= f;}
 218
 219 inline Vec4 operator*(const Vec4& v1, double f)
 220   {Vec4 v = v1; return v *= f;}
 221
 222 inline Vec4 operator/(const Vec4& v1, double f)
 223   {Vec4 v = v1; return v /= f;}
 224
 225 inline double operator*(const Vec4& v1, const Vec4& v2)
 226   {return v1.tt*v2.tt - v1.xx*v2.xx - v1.yy*v2.yy - v1.zz*v2.zz;}
 227
 228 //**************************************************************************
 229
 230 // RotBstMatrix class.
 231 // This class implements 4 * 4 matrices that encode an arbitrary combination
 232 // of rotations and boosts, that can be applied to Vec4 four-vectors.
 233
 234 class RotBstMatrix {
 235
 236 public:
 237
 238   // Constructors.
 239   RotBstMatrix() {for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j)
 240     { M[i][j] = (i==j) ? 1. : 0.; } } }
 241   RotBstMatrix(const RotBstMatrix& Min) {
 242     for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
 243     M[i][j] = Min.M[i][j]; } } }
 244   RotBstMatrix& operator=(const RotBstMatrix& Min) {if (this != &Min) {
 245     for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
 246     M[i][j] = Min.M[i][j]; } } } return *this; }
 247
 248   // Member functions.
 249   void rot(double = 0., double = 0.);
 250   void rot(const Vec4& p);
 251   void bst(double = 0., double = 0., double = 0.);
 252   void bst(const Vec4&);
 253   void bstback(const Vec4&);
 254   void bst(const Vec4&, const Vec4&);
 255   void toCMframe(const Vec4&, const Vec4&);
 256   void fromCMframe(const Vec4&, const Vec4&);
 257   void rotbst(const RotBstMatrix&);
 258   void invert();
 259   void reset();
 260
 261   // Crude estimate deviation from unit matrix.
 262   double deviation() const;
 263
 264   // Print a transformation matrix.
 265   friend ostream& operator<<(ostream&, const RotBstMatrix&) ;
 266
 267   // Private members to be accessible from Vec4.
 268   friend class Vec4;
 269
 270 private:
 271
 272   // Constants: could only be changed in the code itself.
 273   static const double TINY;
 274
 275   // The rotation-and-boost matrix data members.
 276   double M[4][4];
 277
 278 };
 279
 280 //**************************************************************************
 281
 282 // Hist class.
 283 // This class handles a single histogram at a time.
 284
 285 class Hist{
 286
 287 public:
 288
 289   // Constructors, including copy constructors.
 290   Hist() {;}
 291   Hist(string titleIn, int nBinIn = 100, double xMinIn = 0.,
 292     double xMaxIn = 1.) {
 293     book(titleIn, nBinIn, xMinIn, xMaxIn);}
 294   Hist(const Hist& h)
 295     : title(h.title), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
 296     xMax(h.xMax), dx(h.dx), under(h.under), inside(h.inside),
 297     over(h.over), res(h.res) { }
 298   Hist(string titleIn, const Hist& h)
 299     : title(titleIn), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
 300     xMax(h.xMax), dx(h.dx), under(h.under), inside(h.inside),
 301     over(h.over), res(h.res) { }
 302   Hist& operator=(const Hist& h) { if(this != &h) {
 303     nBin = h.nBin; nFill = h.nFill; xMin = h.xMin; xMax = h.xMax;
 304     dx = h.dx;  under = h.under; inside = h.inside; over = h.over;
 305     res = h.res; } return *this; }
 306
 307   // Book a histogram.
 308   void book(string titleIn = "  ", int nBinIn = 100, double xMinIn = 0.,
 309     double xMaxIn = 1.) ;
 310
 311   // Set title of a histogram.
 312   void name(string titleIn = "  ") {title = titleIn; }
 313
 314   // Reset bin contents.
 315   void null() ;
 316
 317   // Fill bin with weight.
 318   void fill(double x, double w = 1.) ;
 319
 320   // Return content of specific bin: -1 gives underflow and nBin overflow.
 321   double getBinContent(int iBin) ;
 322
 323   // Return number of entries
 324   double getEntries() {return nFill; }
 325
 326   // Print histogram contents as a table (e.g. for Gnuplot).
 327   void table(ostream& os = cout) const ;
 328
 329   // Check whether another histogram has same size and limits.
 330   bool sameSize(const Hist& h) const ;
 331
 332   // Take logarithm (base 10 or e) of bin contents.
 333   void takeLog(bool tenLog = true) ;
 334
 335   // Operator overloading with member functions
 336   Hist& operator+=(const Hist& h) ;
 337   Hist& operator-=(const Hist& h) ;
 338   Hist& operator*=(const Hist& h) ;
 339   Hist& operator/=(const Hist& h) ;
 340   Hist& operator+=(double f) ;
 341   Hist& operator-=(double f) ;
 342   Hist& operator*=(double f) ;
 343   Hist& operator/=(double f) ;
 344
 345   // Operator overloading with friends
 346   friend Hist operator+(double f, const Hist& h1);
 347   friend Hist operator+(const Hist& h1, double f);
 348   friend Hist operator+(const Hist& h1, const Hist& h2);
 349   friend Hist operator-(double f, const Hist& h1);
 350   friend Hist operator-(const Hist& h1, double f);
 351   friend Hist operator-(const Hist& h1, const Hist& h2);
 352   friend Hist operator*(double f, const Hist& h1);
 353   friend Hist operator*(const Hist& h1, double f);
 354   friend Hist operator*(const Hist& h1, const Hist& h2);
 355   friend Hist operator/(double f, const Hist& h1);
 356   friend Hist operator/(const Hist& h1, double f);
 357   friend Hist operator/(const Hist& h1, const Hist& h2);
 358
 359   // Print a histogram with overloaded << operator.
 360   friend ostream& operator<<(ostream& os, const Hist& h) ;
 361
 362 private:
 363
 364   // Constants: could only be changed in the code itself.
 365   static const int    NBINMAX, NLINES;
 366   static const double TOLERANCE, TINY, SMALLFRAC, DYAC[];
 367   static const char   NUMBER[];
 368
 369   // Properties and contents of a histogram.
 370   string title;
 371   int    nBin, nFill;
 372   double xMin, xMax, dx, under, inside, over;
 373   vector<double> res;
 374
 375 };
 376
 377 //**************************************************************************
 378
 379 } // end namespace Pythia8
 380
 381 #endif // end Pythia8_Basics_H