PYTHIA8/pythia8130/src/Basics.cxx

   1 // Basics.cc 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 // Function definitions (not found in the header) for the Rndm, Vec4,
   7 // RotBstMatrix and Hist classes, and some related global functions.
   8
   9 #include "Basics.h"
  10
  11 // Access time information.
  12 #include <ctime>
  13
  14 namespace Pythia8 {
  15
  16 //**************************************************************************
  17
  18 // Rndm class.
  19 // This class handles random number generation according to the
  20 // Marsaglia-Zaman-Tsang algorithm
  21
  22 //*********
  23
  24 // Definitions of static variables.
  25
  26 bool   Rndm::initRndm        = false;
  27 bool   Rndm::saveGauss       = false;
  28 int    Rndm::i97, Rndm::j97;
  29 int    Rndm::defaultSeed     = 19780503;
  30 double Rndm::u[97], Rndm::c, Rndm::cd, Rndm::cm, Rndm::save;
  31 bool   Rndm::useExternalRndm = false;
  32 RndmEngine* Rndm::rndmPtr    = 0;
  33
  34 //*********
  35
  36 // Method to pass in pointer for external random number generation.
  37
  38 bool Rndm::rndmEnginePtr( RndmEngine* rndmPtrIn) {
  39
  40   // Save pointer.
  41   if (rndmPtrIn == 0) return false;
  42   rndmPtr         = rndmPtrIn;
  43   useExternalRndm = true;
  44
  45   // Done.
  46   return true;
  47
  48 }
  49
  50 //*********
  51
  52 // Initialize, normally at construction or in first call.
  53
  54 void Rndm::init(int seedIn) {
  55
  56   // Pick seed in convenient way. Assure it to be non-negative.
  57   int seed = seedIn;
  58   if (seedIn < 0) seed = defaultSeed;
  59   else if (seedIn == 0) seed = int(time(0));
  60   if (seed < 0) seed = -seed;
  61
  62   // Unpack seed.
  63   int ij = (seed/30082) % 31329;
  64   int kl = seed % 30082;
  65   int i  = (ij/177) % 177 + 2;
  66   int j  = ij % 177 + 2;
  67   int k  = (kl/169) % 178 + 1;
  68   int l  =  kl % 169;
  69
  70   // Initialize random number array.
  71   for (int ii = 0; ii < 97; ++ii) {
  72     double s = 0.;
  73     double t = 0.5;
  74     for (int jj = 0; jj < 48; ++jj) {
  75       int m = (( (i*j)%179 )*k) % 179;
  76       i = j;
  77       j = k;
  78       k = m;
  79       l = (53*l+1) % 169;
  80       if ( (l*m) % 64 >= 32) s += t;
  81       t *= 0.5;
  82     }
  83     u[ii] = s;
  84   }
  85
  86   // Initialize other variables.
  87   double twom24 = 1.;
  88   for (int i24 = 0; i24 < 24; ++i24) twom24 *= 0.5;
  89   c   = 362436. * twom24;
  90   cd  = 7654321. * twom24;
  91   cm  = 16777213. * twom24;
  92   i97 = 96;
  93   j97 = 32;
  94
  95   // Finished.
  96   initRndm = true;
  97
  98 }
  99
 100 //*********
 101
 102 // Generate next random number uniformly between 0 and 1.
 103
 104 double Rndm::flat() {
 105
 106   // Use external random number generator if such has been linked.
 107   if(useExternalRndm) return rndmPtr->flat();
 108
 109   // Ensure that already initialized.
 110   if (!initRndm) init(defaultSeed);
 111
 112   // Find next random number and update saved state.
 113   double uni;
 114   do {
 115     uni = u[i97] - u[j97];
 116     if (uni < 0.) uni += 1.;
 117     u[i97] = uni;
 118     if (--i97 < 0) i97 = 96;
 119     if (--j97 < 0) j97 = 96;
 120     c -= cd;
 121     if (c < 0.) c += cm;
 122     uni -= c;
 123     if(uni < 0.) uni += 1.;
 124    } while (uni <= 0. || uni >= 1.);
 125   return uni;
 126
 127 }
 128
 129 //*********
 130
 131 // Generate random numbers according to exp(-x^2/2).
 132
 133 double Rndm::gauss() {
 134
 135   // Generate pair of Gaussian random numbers.
 136   if (!saveGauss) {
 137     saveGauss  = true;
 138     double r   = sqrt(-2. * log(flat()));
 139     double phi = 2. * M_PI * flat();
 140     save       = r * sin(phi);
 141     return       r * cos(phi);
 142
 143   // Use saved element of pair.
 144   } else {
 145     saveGauss  = false;
 146     return       save;
 147   }
 148
 149 }
 150
 151 //*********
 152
 153 // Pick one option among  vector of (positive) probabilities.
 154
 155 int Rndm::pick(const vector<double>& prob) {
 156
 157   double work = 0.;
 158   for (int i = 0; i < int(prob.size()); ++i) {work += prob[i];}
 159   work *= flat();
 160   int index = -1;
 161   do { work -= prob[++index]; } while (work > 0 && index < int(prob.size()));
 162   return index;
 163
 164 }
 165
 166 //**************************************************************************
 167
 168 // Vec4 class.
 169 // This class implements four-vectors, in energy-momentum space.
 170 // (But could also be used to hold space-time four-vectors.)
 171
 172 //*********
 173
 174 // Constants: could be changed here if desired, but normally should not.
 175 // These are of technical nature, as described for each.
 176
 177 // Small number to avoid division by zero.
 178 const double Vec4::TINY = 1e-20;
 179
 180 //*********
 181
 182 // Rotation (simple).
 183
 184 void Vec4::rot(double thetaIn, double phiIn) {
 185
 186   double cthe = cos(thetaIn);
 187   double sthe = sin(thetaIn);
 188   double cphi = cos(phiIn);
 189   double sphi = sin(phiIn);
 190   double tmpx =  cthe * cphi * xx -    sphi * yy + sthe * cphi * zz;
 191   double tmpy =  cthe * sphi * xx +    cphi * yy + sthe * sphi * zz;
 192   double tmpz = -sthe *        xx +                cthe *        zz;
 193   xx          = tmpx;
 194   yy          = tmpy;
 195   zz          = tmpz;
 196
 197 }
 198
 199 //*********
 200
 201 // Azimuthal rotation phi around an arbitrary axis (nz, ny, nz).
 202
 203 void Vec4::rotaxis(double phiIn, double nx, double ny, double nz) {
 204
 205   double norm = 1./sqrt(nx*nx + ny*ny + nz*nz);
 206   nx         *= norm;
 207   ny         *= norm;
 208   nz         *= norm;
 209   double cphi = cos(phiIn);
 210   double sphi = sin(phiIn);
 211   double comb = (nx * xx + ny * yy + nz * zz) * (1. - cphi);
 212   double tmpx = cphi * xx + comb * nx + sphi * (ny * zz - nz * yy);
 213   double tmpy = cphi * yy + comb * ny + sphi * (nz * xx - nx * zz);
 214   double tmpz = cphi * zz + comb * nz + sphi * (nx * yy - ny * xx);
 215   xx          = tmpx;
 216   yy          = tmpy;
 217   zz          = tmpz;
 218
 219 }
 220
 221 //*********
 222
 223 // Azimuthal rotation phi around an arbitrary (3-vector component of) axis.
 224
 225 void Vec4::rotaxis(double phiIn, const Vec4& n) {
 226
 227   double nx   = n.xx;
 228   double ny   = n.yy;
 229   double nz   = n.zz;
 230   double norm = 1./sqrt(nx*nx + ny*ny + nz*nz);
 231   nx         *= norm;
 232   ny          *=norm;
 233   nz          *=norm;
 234   double cphi = cos(phiIn);
 235   double sphi = sin(phiIn);
 236   double comb = (nx * xx + ny * yy + nz * zz) * (1. - cphi);
 237   double tmpx = cphi * xx + comb * nx + sphi * (ny * zz - nz * yy);
 238   double tmpy = cphi * yy + comb * ny + sphi * (nz * xx - nx * zz);
 239   double tmpz = cphi * zz + comb * nz + sphi * (nx * yy - ny * xx);
 240   xx          = tmpx;
 241   yy          = tmpy;
 242   zz          = tmpz;
 243
 244 }
 245
 246 //*********
 247
 248 // Boost (simple).
 249
 250 void Vec4::bst(double betaX, double betaY, double betaZ) {
 251
 252   double beta2 = betaX*betaX + betaY*betaY + betaZ*betaZ;
 253   double gamma = 1. / sqrt(1. - beta2);
 254   double prod1 = betaX * xx + betaY * yy + betaZ * zz;
 255   double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
 256   xx += prod2 * betaX;
 257   yy += prod2 * betaY;
 258   zz += prod2 * betaZ;
 259   tt = gamma * (tt + prod1);
 260
 261 }
 262
 263 //*********
 264
 265 // Boost (simple, given gamma).
 266
 267 void Vec4::bst(double betaX, double betaY, double betaZ, double gamma) {
 268
 269   double prod1 = betaX * xx + betaY * yy + betaZ * zz;
 270   double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
 271   xx += prod2 * betaX;
 272   yy += prod2 * betaY;
 273   zz += prod2 * betaZ;
 274   tt = gamma * (tt + prod1);
 275
 276 }
 277
 278 //*********
 279
 280 // Boost given by a Vec4 p.
 281
 282 void Vec4::bst(const Vec4& pIn) {
 283
 284   double betaX = pIn.xx / pIn.tt;
 285   double betaY = pIn.yy / pIn.tt;
 286   double betaZ = pIn.zz / pIn.tt;
 287   double beta2 = betaX*betaX + betaY*betaY + betaZ*betaZ;
 288   double gamma = 1. / sqrt(1. - beta2);
 289   double prod1 = betaX * xx + betaY * yy + betaZ * zz;
 290   double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
 291   xx          += prod2 * betaX;
 292   yy          += prod2 * betaY;
 293   zz          += prod2 * betaZ;
 294   tt           = gamma * (tt + prod1);
 295
 296 }
 297
 298 //*********
 299
 300 // Boost given by a Vec4 p and double m.
 301
 302 void Vec4::bst(const Vec4& pIn, double mIn) {
 303
 304   double betaX = pIn.xx / pIn.tt;
 305   double betaY = pIn.yy / pIn.tt;
 306   double betaZ = pIn.zz / pIn.tt;
 307   double gamma = pIn.tt / mIn;
 308   double prod1 = betaX * xx + betaY * yy + betaZ * zz;
 309   double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
 310   xx          += prod2 * betaX;
 311   yy          += prod2 * betaY;
 312   zz          += prod2 * betaZ;
 313   tt           = gamma * (tt + prod1);
 314
 315 }
 316
 317 //*********
 318
 319 // Boost given by a Vec4 p; boost in opposite direction.
 320
 321 void Vec4::bstback(const Vec4& pIn) {
 322
 323   double betaX = -pIn.xx / pIn.tt;
 324   double betaY = -pIn.yy / pIn.tt;
 325   double betaZ = -pIn.zz / pIn.tt;
 326   double beta2 = betaX*betaX + betaY*betaY + betaZ*betaZ;
 327   double gamma = 1. / sqrt(1. - beta2);
 328   double prod1 = betaX * xx + betaY * yy + betaZ * zz;
 329   double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
 330   xx          += prod2 * betaX;
 331   yy          += prod2 * betaY;
 332   zz          += prod2 * betaZ;
 333   tt           = gamma * (tt + prod1);
 334
 335 }
 336
 337 //*********
 338
 339 // Boost given by a Vec4 p and double m; boost in opposite direction.
 340
 341 void Vec4::bstback(const Vec4& pIn, double mIn) {
 342
 343   double betaX = -pIn.xx / pIn.tt;
 344   double betaY = -pIn.yy / pIn.tt;
 345   double betaZ = -pIn.zz / pIn.tt;
 346   double gamma = pIn.tt / mIn;
 347   double prod1 = betaX * xx + betaY * yy + betaZ * zz;
 348   double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
 349   xx          += prod2 * betaX;
 350   yy          += prod2 * betaY;
 351   zz          += prod2 * betaZ;
 352   tt           = gamma * (tt + prod1);
 353
 354 }
 355
 356 //*********
 357
 358 // Arbitrary combination of rotations and boosts defined by 4 * 4 matrix.
 359
 360 void Vec4::rotbst(const RotBstMatrix& M) {
 361
 362   double x = xx; double y = yy; double z = zz; double t = tt;
 363   tt = M.M[0][0] * t + M.M[0][1] * x + M.M[0][2] * y +  M.M[0][3] * z;
 364   xx = M.M[1][0] * t + M.M[1][1] * x + M.M[1][2] * y +  M.M[1][3] * z;
 365   yy = M.M[2][0] * t + M.M[2][1] * x + M.M[2][2] * y +  M.M[2][3] * z;
 366   zz = M.M[3][0] * t + M.M[3][1] * x + M.M[3][2] * y +  M.M[3][3] * z;
 367
 368 }
 369
 370 //*********
 371
 372 // The invariant mass of two four-vectors.
 373
 374 double m(const Vec4& v1, const Vec4& v2) {
 375   double m2 = pow2(v1.tt + v2.tt) - pow2(v1.xx + v2.xx)
 376      - pow2(v1.yy + v2.yy) - pow2(v1.zz + v2.zz);
 377   return (m2 > 0.) ? sqrt(m2) : 0.;
 378 }
 379
 380 //*********
 381
 382 // The squared invariant mass of two four-vectors.
 383
 384 double m2(const Vec4& v1, const Vec4& v2) {
 385   double m2 = pow2(v1.tt + v2.tt) - pow2(v1.xx + v2.xx)
 386      - pow2(v1.yy + v2.yy) - pow2(v1.zz + v2.zz);
 387   return m2;
 388 }
 389
 390 //*********
 391
 392 // The scalar product of two three-vectors.
 393
 394 double dot3(const Vec4& v1, const Vec4& v2) {
 395   return v1.xx*v2.xx + v1.yy*v2.yy + v1.zz*v2.zz;
 396 }
 397
 398 //*********
 399
 400 // The cross product of two three-vectors.
 401
 402 Vec4 cross3(const Vec4& v1, const Vec4& v2) {
 403   Vec4 v;
 404   v.xx = v1.yy * v2.zz - v1.zz * v2.yy;
 405   v.yy = v1.zz * v2.xx - v1.xx * v2.zz;
 406   v.zz = v1.xx * v2.yy - v1.yy * v2.xx; return v;
 407 }
 408
 409 //*********
 410
 411 // Opening angle between two three-vectors.
 412
 413 double theta(const Vec4& v1, const Vec4& v2) {
 414   double cthe = (v1.xx * v2.xx + v1.yy * v2.yy + v1.zz * v2.zz)
 415     / sqrt( (v1.xx*v1.xx + v1.yy*v1.yy + v1.zz*v1.zz)
 416     * (v2.xx*v2.xx + v2.yy*v2.yy + v2.zz*v2.zz) );
 417   cthe = max(-1., min(1., cthe));
 418   return acos(cthe);
 419 }
 420
 421 //*********
 422
 423 // Cosine of the opening angle between two three-vectors.
 424
 425 double costheta(const Vec4& v1, const Vec4& v2) {
 426   double cthe = (v1.xx * v2.xx + v1.yy * v2.yy + v1.zz * v2.zz)
 427     / sqrt( (v1.xx*v1.xx + v1.yy*v1.yy + v1.zz*v1.zz)
 428     * (v2.xx*v2.xx + v2.yy*v2.yy + v2.zz*v2.zz) );
 429   cthe = max(-1., min(1., cthe));
 430   return cthe;
 431 }
 432
 433 //*********
 434
 435 // Azimuthal angle between two three-vectors.
 436
 437 double phi(const Vec4& v1, const Vec4& v2) {
 438   double cphi = (v1.xx * v2.xx + v1.yy * v2.yy) / sqrt( max( Vec4::TINY,
 439     (v1.xx*v1.xx + v1.yy*v1.yy) * (v2.xx*v2.xx + v2.yy*v2.yy) ));
 440   cphi = max(-1., min(1., cphi));
 441   return acos(cphi);
 442 }
 443
 444 //*********
 445
 446 // Cosine of the azimuthal angle between two three-vectors.
 447
 448 double cosphi(const Vec4& v1, const Vec4& v2) {
 449   double cphi = (v1.xx * v2.xx + v1.yy * v2.yy) / sqrt( max( Vec4::TINY,
 450     (v1.xx*v1.xx + v1.yy*v1.yy) * (v2.xx*v2.xx + v2.yy*v2.yy) ));
 451   cphi = max(-1., min(1., cphi));
 452   return cphi;
 453 }
 454
 455 //*********
 456
 457 // Azimuthal angle between two three-vectors around a third.
 458
 459 double phi(const Vec4& v1, const Vec4& v2, const Vec4& n) {
 460   double nx = n.xx; double ny = n.yy; double nz = n.zz;
 461   double norm = 1. / sqrt(nx*nx + ny*ny + nz*nz);
 462   nx *= norm; ny *=norm; nz *=norm;
 463   double v1s = v1.xx * v1.xx + v1.yy * v1.yy + v1.zz * v1.zz;
 464   double v2s = v2.xx * v2.xx + v2.yy * v2.yy + v2.zz * v2.zz;
 465   double v1v2 = v1.xx * v2.xx + v1.yy * v2.yy + v1.zz * v2.zz;
 466   double v1n = v1.xx * nx + v1.yy * ny + v1.zz * nz;
 467   double v2n = v2.xx * nx + v2.yy * ny + v2.zz * nz;
 468   double cphi = (v1v2 - v1n * v2n) / sqrt( max( Vec4::TINY,
 469     (v1s - v1n*v1n) * (v2s - v2n*v2n) ));
 470   cphi = max(-1., min(1., cphi));
 471   return acos(cphi);
 472 }
 473
 474 //*********
 475
 476 // Cosine of the azimuthal angle between two three-vectors around a third.
 477
 478 double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n) {
 479   double nx = n.xx; double ny = n.yy; double nz = n.zz;
 480   double norm = 1. / sqrt(nx*nx + ny*ny + nz*nz);
 481   nx *= norm; ny *=norm; nz *=norm;
 482   double v1s = v1.xx * v1.xx + v1.yy * v1.yy + v1.zz * v1.zz;
 483   double v2s = v2.xx * v2.xx + v2.yy * v2.yy + v2.zz * v2.zz;
 484   double v1v2 = v1.xx * v2.xx + v1.yy * v2.yy + v1.zz * v2.zz;
 485   double v1n = v1.xx * nx + v1.yy * ny + v1.zz * nz;
 486   double v2n = v2.xx * nx + v2.yy * ny + v2.zz * nz;
 487   double cphi = (v1v2 - v1n * v2n) / sqrt( max( Vec4::TINY,
 488     (v1s - v1n*v1n) * (v2s - v2n*v2n) ));
 489   cphi = max(-1., min(1., cphi));
 490   return cphi;
 491 }
 492
 493 //*********
 494
 495 // Print a four-vector: also operator overloading with friend.
 496
 497 ostream& operator<<(ostream& os, const Vec4& v) {
 498   os << fixed << setprecision(3) << setw(10) << v.xx << setw(10) << v.yy
 499      << setw(10) << v.zz << setw(10) << v.tt << "\n";
 500   return os;
 501 }
 502
 503 //**************************************************************************
 504
 505 // RotBstMatrix class.
 506 // This class implements 4 * 4 matrices that encode an arbitrary combination
 507 // of rotations and boosts, that can be applied to Vec4 four-vectors.
 508
 509 //*********
 510
 511 // Constants: could be changed here if desired, but normally should not.
 512 // These are of technical nature, as described for each.
 513
 514 // Small number to avoid division by zero.
 515 const double RotBstMatrix::TINY = 1e-20;
 516
 517 //*********
 518
 519 // Rotate by polar angle theta and azimuthal angle phi.
 520
 521 void RotBstMatrix::rot(double theta, double phi) {
 522
 523   // Set up rotation matrix.
 524   double cthe = cos(theta); double sthe = sin(theta);
 525   double cphi = cos(phi); double sphi = sin(phi);
 526   double Mrot[4][4] = {
 527     {1.,           0.,         0.,          0.},
 528     {0.,  cthe * cphi,     - sphi, sthe * cphi},
 529     {0.,  cthe * sphi,       cphi, sthe * sphi},
 530     {0., -sthe,                0., cthe       } };
 531
 532   // Rotate current matrix accordingly.
 533   double Mtmp[4][4];
 534   for (int i = 0; i < 4; ++i) {
 535     for (int j = 0; j < 4; ++j) {
 536       Mtmp[i][j] = M[i][j];
 537     }
 538   }
 539   for (int i = 0; i < 4; ++i) {
 540     for (int j = 0; j < 4; ++j) {
 541       M[i][j] = Mrot[i][0] * Mtmp[0][j] + Mrot[i][1] * Mtmp[1][j]
 542         + Mrot[i][2] * Mtmp[2][j] + Mrot[i][3] * Mtmp[3][j];
 543     }
 544   }
 545
 546 }
 547
 548 //*********
 549
 550 // Rotate so that vector originally along z axis becomes parallel with p.
 551
 552 void RotBstMatrix::rot(const Vec4& p) {
 553
 554   double theta = p.theta();
 555   double phi = p.phi();
 556   rot(0., -phi);
 557   rot(theta, phi);
 558
 559 }
 560
 561 //*********
 562
 563 // Boost with velocity vector (betaX, betaY, betaZ).
 564
 565 void RotBstMatrix::bst(double betaX, double betaY, double betaZ) {
 566
 567   // Set up boost matrix.
 568   double gm = 1. / sqrt( max( TINY, 1. - betaX*betaX - betaY*betaY
 569     - betaZ*betaZ ) );
 570   double gf = gm*gm / (1. + gm);
 571   double Mbst[4][4] = {
 572     { gm,           gm*betaX,           gm*betaY,          gm*betaZ },
 573     { gm*betaX, 1. + gf*betaX*betaX, gf*betaX*betaY, gf*betaX*betaZ },
 574     { gm*betaY, gf*betaY*betaX, 1. + gf*betaY*betaY, gf*betaY*betaZ },
 575     { gm*betaZ, gf*betaZ*betaX, gf*betaZ*betaY, 1. + gf*betaZ*betaZ } };
 576
 577   // Boost current matrix correspondingly.
 578   double Mtmp[4][4];
 579   for (int i = 0; i < 4; ++i) {
 580     for (int j = 0; j < 4; ++j) {
 581       Mtmp[i][j] = M[i][j];
 582     }
 583   }
 584   for (int i = 0; i < 4; ++i) {
 585     for (int j = 0; j < 4; ++j) {
 586       M[i][j] = Mbst[i][0] * Mtmp[0][j] + Mbst[i][1] * Mtmp[1][j]
 587         + Mbst[i][2] * Mtmp[2][j] + Mbst[i][3] * Mtmp[3][j];
 588     }
 589   }
 590
 591 }
 592
 593 //*********
 594
 595 // Boost so that vector originally at rest obtains same velocity as p.
 596
 597 void RotBstMatrix::bst(const Vec4& p) {
 598   double betaX = p.px() / p.e();
 599   double betaY = p.py() / p.e();
 600   double betaZ = p.pz() / p.e();
 601   bst(betaX, betaY, betaZ);
 602 }
 603
 604 //*********
 605
 606 // Boost so vector originally with same velocity as p is brought to rest.
 607
 608 void RotBstMatrix::bstback(const Vec4& p) {
 609   double betaX = -p.px() / p.e();
 610   double betaY = -p.py() / p.e();
 611   double betaZ = -p.pz() / p.e();
 612   bst(betaX, betaY, betaZ);
 613 }
 614
 615 //*********
 616
 617 // Boost that transforms p1 to p2, where p1^2 = p2^2 is assumed.
 618
 619 void RotBstMatrix::bst(const Vec4& p1, const Vec4& p2) {
 620   double eSum = p1.e() + p2.e();
 621   double betaX = (p2.px() - p1.px()) / eSum;
 622   double betaY = (p2.py() - p1.py()) / eSum;
 623   double betaZ = (p2.pz() - p1.pz()) / eSum;
 624   double fac = 2. / (1. + betaX*betaX + betaY*betaY + betaZ*betaZ);
 625   betaX *= fac; betaY *= fac; betaZ *= fac;
 626   bst(betaX, betaY, betaZ);
 627 }
 628
 629 //*********
 630
 631 // Boost and rotation that transforms from p1 and p2
 632 // to their rest frame with p1 along +z axis.
 633
 634 void RotBstMatrix::toCMframe(const Vec4& p1, const Vec4& p2) {
 635   Vec4 pSum = p1 + p2;
 636   Vec4 dir = p1;
 637   dir.bstback(pSum);
 638   double theta = dir.theta();
 639   double phi = dir.phi();
 640   bstback(pSum);
 641   rot(0., -phi);
 642   rot(-theta, phi);
 643 }
 644
 645 //*********
 646
 647 // Rotation and boost that transforms from rest frame of p1 and p2
 648 // with p1 along +z axis to actual frame of p1 and p2. (Inverse of above.)
 649
 650 void RotBstMatrix::fromCMframe(const Vec4& p1, const Vec4& p2) {
 651   Vec4 pSum = p1 + p2;
 652   Vec4 dir = p1;
 653   dir.bstback(pSum);
 654   double theta = dir.theta();
 655   double phi = dir.phi();
 656   rot(0., -phi);
 657   rot(theta, phi);
 658   bst(pSum);
 659 }
 660
 661 //*********
 662
 663 // Combine existing rotation/boost matrix with another one.
 664
 665 void RotBstMatrix::rotbst(const RotBstMatrix& Mrb) {
 666   double Mtmp[4][4];
 667   for (int i = 0; i < 4; ++i) {
 668     for (int j = 0; j < 4; ++j) {
 669       Mtmp[i][j] = M[i][j];
 670     }
 671   }
 672   for (int i = 0; i < 4; ++i) {
 673     for (int j = 0; j < 4; ++j) {
 674       M[i][j] = Mrb.M[i][0] * Mtmp[0][j] + Mrb.M[i][1] * Mtmp[1][j]
 675         + Mrb.M[i][2] * Mtmp[2][j] + Mrb.M[i][3] * Mtmp[3][j];
 676     }
 677   }
 678 }
 679
 680 //*********
 681
 682 // Invert the rotation and boost.
 683
 684 void RotBstMatrix::invert() {
 685   double Mtmp[4][4];
 686   for (int i = 0; i < 4; ++i) {
 687     for (int j = 0; j < 4; ++j) {
 688       Mtmp[i][j] = M[i][j];
 689     }
 690   }
 691   for (int i = 0; i < 4; ++i) {
 692     for (int j = 0; j < 4; ++j) {
 693       M[i][j] = ( (i == 0 && j > 0) || (i > 0 && j == 0) )
 694         ? - Mtmp[j][i] : Mtmp[j][i];
 695     }
 696   }
 697 }
 698
 699 //*********
 700
 701 // Reset to diagonal matrix.
 702
 703 void RotBstMatrix::reset() {
 704   for (int i = 0; i < 4; ++i) {
 705     for (int j = 0; j < 4; ++j) {
 706       M[i][j] = (i==j) ? 1. : 0.;
 707     }
 708   }
 709 }
 710
 711 //*********
 712
 713 // Crude estimate deviation from unit matrix.
 714
 715 double RotBstMatrix::deviation() const {
 716   double devSum = 0.;
 717   for (int i = 0; i < 4; ++i) {
 718     for (int j = 0; j < 4; ++j) {
 719       devSum += (i==j) ? abs(M[i][j] - 1.) : abs(M[i][j]);
 720     }
 721   }
 722   return devSum;
 723 }
 724
 725 //*********
 726
 727 // Print a rotation and boost matrix: operator overloading with friend.
 728
 729 ostream& operator<<(ostream& os, const RotBstMatrix& M) {
 730   os << fixed << setprecision(5) << "    Rotation/boost matrix: \n";
 731   for (int i = 0; i <4; ++i) {
 732     os << setw(10) << M.M[i][0] << setw(10) << M.M[i][1]
 733        << setw(10) << M.M[i][2] << setw(10) << M.M[i][3] << "\n";
 734   }
 735   return os;
 736 }
 737
 738 //**************************************************************************
 739
 740 // Hist class.
 741 // This class handles a single histogram at a time
 742 // (or a vector of histograms).
 743
 744 //*********
 745
 746 // Constants: could be changed here if desired, but normally should not.
 747 // These are of technical nature, as described for each.
 748
 749 // Maximum number of bins in a histogram.
 750 const int    Hist::NBINMAX   = 100;
 751
 752 // Maximum number of lines a histogram can use at output.
 753 const int    Hist::NLINES    = 30;
 754
 755 // Tolerance in deviation of xMin and xMax between two histograms.
 756 const double Hist::TOLERANCE = 0.001;
 757
 758 // Small number to avoid division by zero.
 759 const double Hist::TINY      = 1e-20;
 760
 761 // When minbin/maxbin < SMALLFRAC the y scale goes down to zero.
 762 const double Hist::SMALLFRAC = 0.1;
 763
 764 // Constants for printout: fixed steps on y scale; filling characters.
 765 const double DYAC[] = {0.04, 0.05, 0.06, 0.08, 0.10,
 766   0.12, 0.15, 0.20, 0.25, 0.30};
 767 const char NUMBER[] = {'0', '1', '2', '3', '4', '5',
 768   '6', '7', '8', '9', 'X' };
 769
 770 //*********
 771
 772 // Book a histogram.
 773
 774 void Hist::book(string titleIn, int nBinIn, double xMinIn,
 775   double xMaxIn) {
 776
 777   title = titleIn;
 778   nBin = nBinIn;
 779   if (nBinIn < 1) nBin = 1;
 780   if (nBinIn > NBINMAX) nBin = NBINMAX;
 781   xMin = xMinIn;
 782   xMax = xMaxIn;
 783   dx = (xMax - xMin)/nBin;
 784   res.resize(nBin);
 785   null();
 786
 787 }
 788
 789 //*********
 790
 791 // Reset bin contents.
 792
 793 void Hist::null() {
 794
 795   nFill = 0;
 796   under = 0.;
 797   inside = 0.;
 798   over = 0.;
 799   for (int ix = 0; ix < nBin; ++ix) {res[ix] = 0.;}
 800
 801 }
 802
 803 //*********
 804
 805 // Fill bin with weight.
 806
 807 void Hist::fill(double x, double w) {
 808
 809   ++nFill;
 810   int iBin = int(floor((x - xMin)/dx));
 811   if (iBin < 0) {under += w; }
 812   else if (iBin >= nBin) {over += w; }
 813   else {inside += w; res[iBin] += w; }
 814
 815 }
 816
 817 //*********
 818
 819 // Get content of specific bin.
 820 // Special values are bin 0 for underflow and bin nBin+1 for overflow.
 821 // All other bins outside proper histogram range return 0.
 822
 823 double Hist::getBinContent(int iBin) {
 824
 825   if (iBin > 0 && iBin <= nBin) return res[iBin - 1];
 826   else if (iBin == 0) return under;
 827   else if (iBin == nBin + 1) return over;
 828   else return 0.;
 829
 830 }
 831
 832 //*********
 833
 834 // Print histogram contents as a table (e.g. for Gnuplot).
 835
 836 void Hist::table(ostream& os) const {
 837
 838   // Print histogram vector bin by bin, with mean x as first column.
 839   os << scientific << setprecision(4);
 840   for (int ix = 0; ix < nBin; ++ix) {
 841     os << setw(12) << xMin + (ix + 0.5) * dx
 842        << setw(12) << res[ix] << "\n";
 843
 844   }
 845 }
 846
 847 //*********
 848
 849 // Check whether another histogram has same size and limits.
 850
 851 bool Hist::sameSize(const Hist& h) const {
 852
 853   if (nBin == h.nBin && abs(xMin - h.xMin) < TOLERANCE * dx &&
 854     abs(xMax - h.xMax) < TOLERANCE * dx) {return true;}
 855   else {return false;}
 856
 857 }
 858
 859 //*********
 860
 861 // Take 10-logarithm or natural logarithm of contents bin by bin.
 862
 863 void Hist::takeLog(bool tenLog) {
 864
 865   // Find smallest positive bin content, and put min a bit below.
 866   double yMin = 1e20;
 867   for (int ix = 0; ix < nBin; ++ix)
 868     if (res[ix] > 1e-20 && res[ix] < yMin ) yMin = res[ix];
 869   yMin *= 0.8;
 870
 871   // Take 10-logarithm bin by bin, but ensure positivity.
 872   if (tenLog) {
 873     for (int ix = 0; ix < nBin; ++ix) res[ix] = log10( max( yMin, res[ix]) );
 874     under  =  log10( max( yMin, under) );
 875     inside =  log10( max( yMin, inside) );
 876     over   =  log10( max( yMin, over) );
 877
 878   // Take natural logarithm bin by bin, but ensure positivity.
 879   } else {
 880     for (int ix = 0; ix < nBin; ++ix) res[ix] = log( max( yMin, res[ix]) );
 881     under  =  log( max( yMin, under) );
 882     inside =  log( max( yMin, inside) );
 883     over   =  log( max( yMin, over) );
 884   }
 885
 886 }
 887
 888 //*********
 889
 890 // Add histogram to existing one.
 891
 892 Hist& Hist::operator+=(const Hist& h) {
 893   if (!sameSize(h)) return *this;
 894   nFill += h.nFill;
 895   under += h.under;
 896   inside += h.inside;
 897   over += h.over;
 898   for (int ix = 0; ix < nBin; ++ix) {res[ix] += h.res[ix];}
 899   return *this;
 900 }
 901
 902 //*********
 903
 904 // Subtract histogram from existing one.
 905
 906 Hist& Hist::operator-=(const Hist& h) {
 907   if (!sameSize(h)) return *this;
 908   nFill += h.nFill;
 909   under -= h.under;
 910   inside -= h.inside;
 911   over -= h.over;
 912   for (int ix = 0; ix < nBin; ++ix) {res[ix] -= h.res[ix];}
 913   return *this;
 914 }
 915
 916 //*********
 917
 918 // Multiply existing histogram by another one.
 919
 920 Hist& Hist::operator*=(const Hist& h) {
 921   if (!sameSize(h)) return *this;
 922   nFill += h.nFill;
 923   under *= h.under;
 924   inside *= h.inside;
 925   over *= h.over;
 926   for (int ix = 0; ix < nBin; ++ix) {res[ix] *= h.res[ix];}
 927   return *this;
 928 }
 929
 930 //*********
 931
 932 // Divide existing histogram by another one.
 933
 934 Hist& Hist::operator/=(const Hist& h) {
 935   if (!sameSize(h)) return *this;
 936   nFill += h.nFill;
 937   under = (abs(h.under) < Hist::TINY) ? 0. : under/h.under;
 938   inside = (abs(h.inside) < Hist::TINY) ? 0. : inside/h.inside;
 939   over = (abs(h.over) < Hist::TINY) ? 0. : over/h.over;
 940   for (int ix = 0; ix < nBin; ++ix) {
 941     res[ix] = (abs(h.res[ix]) < Hist::TINY) ? 0. : res[ix]/h.res[ix];
 942   }
 943   return *this;
 944 }
 945
 946 //*********
 947
 948 // Add constant offset to histogram.
 949
 950 Hist& Hist::operator+=(double f) {
 951   under += f;
 952   inside += nBin * f;
 953   over -= f;
 954   for (int ix = 0; ix < nBin; ++ix) {res[ix] += f;}
 955   return *this;
 956 }
 957
 958 //*********
 959
 960 // Subtract constant offset from histogram.
 961
 962 Hist& Hist::operator-=(double f) {
 963   under -= f;
 964   inside -= nBin * f;
 965   over -= f;
 966   for (int ix = 0; ix < nBin; ++ix) {res[ix] -= f;}
 967   return *this;
 968 }
 969
 970 //*********
 971
 972 // Multiply histogram by constant
 973
 974 Hist& Hist::operator*=(double f) {
 975   under *= f;
 976   inside *= f;
 977   over *= f;
 978   for (int ix = 0; ix < nBin; ++ix) {res[ix] *= f;}
 979   return *this;
 980 }
 981
 982 //*********
 983
 984 // Divide histogram by constant
 985
 986 Hist& Hist::operator/=(double f) {
 987   under /= f;
 988   inside /= f;
 989   over /= f;
 990   for (int ix = 0; ix < nBin; ++ix) {res[ix] /= f;}
 991   return *this;
 992 }
 993
 994 //*********
 995
 996 // Implementation of operator overloading with friends.
 997
 998 Hist operator+(double f, const Hist& h1)
 999   {Hist h = h1; return h += f;}
1000
1001 Hist operator+(const Hist& h1, double f)
1002   {Hist h = h1; return h += f;}
1003
1004 Hist operator+(const Hist& h1, const Hist& h2)
1005   {Hist h = h1; return h += h2;}
1006
1007 Hist operator-(double f, const Hist& h1)
1008   {Hist h = h1;
1009   h.under = f - h1.under;
1010   h.inside = h1.nBin * f - h1.inside;
1011   h.over = f - h1.over;
1012   for (int ix = 0; ix < h1.nBin; ++ix) {h.res[ix] = f - h1.res[ix];}
1013   return h;}
1014
1015 Hist operator-(const Hist& h1, double f)
1016   {Hist h = h1; return h -= f;}
1017
1018 Hist operator-(const Hist& h1, const Hist& h2)
1019   {Hist h = h1; return h -= h2;}
1020
1021 Hist operator*(double f, const Hist& h1)
1022   {Hist h = h1; return h *= f;}
1023
1024 Hist operator*(const Hist& h1, double f)
1025   {Hist h = h1; return h *= f;}
1026
1027 Hist operator*(const Hist& h1, const Hist& h2)
1028   {Hist h = h1; return h *= h2;}
1029
1030 Hist operator/(double f, const Hist& h1) {Hist h = h1;
1031   h.under = (abs(h1.under) < Hist::TINY) ? 0. :  f/h1.under;
1032   h.inside = (abs(h1.inside) < Hist::TINY) ? 0. :  f/h1.inside;
1033   h.over = (abs(h1.over) < Hist::TINY) ? 0. :  f/h1.over;
1034   for (int ix = 0; ix < h1.nBin; ++ix) {
1035     h.res[ix] = (abs(h1.res[ix]) < Hist::TINY) ? 0. : f/h1.res[ix];
1036   }
1037   return h;
1038 }
1039
1040 Hist operator/(const Hist& h1, double f)
1041   {Hist h = h1; return h /= f;}
1042
1043 Hist operator/(const Hist& h1, const Hist& h2)
1044   {Hist h = h1; return h /= h2;}
1045
1046 //*********
1047
1048 // Print a histogram: also operator overloading with friend.
1049
1050 ostream& operator<<(ostream& os, const Hist& h) {
1051
1052   // Do not print empty histograms.
1053   if (h.nFill <= 0) return os;
1054
1055   // Write time and title.
1056   time_t t = time(0);
1057   char date[18];
1058   strftime(date,18,"%Y-%m-%d %H:%M",localtime(&t));
1059   os << "\n\n  " << date << "       " << h.title << "\n\n";
1060
1061   // Find minimum and maximum bin content
1062   double yMin = h.res[0];
1063   double yMax = h.res[0];
1064   for (int i = 1; i < h.nBin; ++i) {
1065     if (h.res[i] < yMin) yMin = h.res[i];
1066     if (h.res[i] > yMax) yMax = h.res[i];
1067   }
1068
1069   // Determine scale and step size for y axis.
1070   if (yMax - yMin > Hist::NLINES * DYAC[0] * 1e-9) {
1071     if (yMin > 0. && yMin < Hist::SMALLFRAC * yMax) yMin = 0.;
1072     if (yMax < 0. && yMax > Hist::SMALLFRAC * yMin) yMax = 0.;
1073     int iPowY = int(floor( log10(yMax - yMin) ));
1074     if (yMax - yMin < Hist::NLINES * DYAC[0] * pow(10.,iPowY))
1075       iPowY = iPowY - 1;
1076     if (yMax - yMin > Hist::NLINES * DYAC[9] * pow(10.,iPowY))
1077       iPowY = iPowY + 1;
1078     double nLinePow = Hist::NLINES * pow(10.,iPowY);
1079     double delY = DYAC[0];
1080     for (int idel = 0; idel < 9; ++idel) {
1081       if (yMax - yMin >= nLinePow * DYAC[idel]) delY = DYAC[idel+1];
1082     }
1083     double dy = delY * pow(10.,iPowY);
1084
1085     // Convert bin contents to integer form; fractional fill in top row.
1086     vector<int> row(h.nBin);
1087     vector<int> frac(h.nBin);
1088     for (int ix = 0; ix < h.nBin ; ++ix) {
1089       double cta = abs(h.res[ix]) / dy;
1090       row[ix] = int(cta + 0.95);
1091       if(h.res[ix] < 0.) row[ix] = - row[ix];
1092       frac[ix] = int(10. * (cta + 1.05 - floor(cta + 0.95)));
1093     }
1094     int rowMin = int(abs(yMin)/dy + 0.95);
1095     if ( yMin < 0) rowMin = - rowMin;
1096     int rowMax = int(abs(yMax)/dy + 0.95);
1097     if ( yMax < 0) rowMax = - rowMax;
1098
1099     // Print histogram row by row.
1100     os << fixed << setprecision(2);
1101     for (int iRow = rowMax; iRow >= rowMin; iRow--) {
1102       if (iRow != 0) {
1103         os << "  " << setw(10) << iRow*delY << "*10^"
1104            << setw(2) << iPowY << "  ";
1105         for (int ix = 0; ix < h.nBin ; ++ix) {
1106           if (iRow == row[ix]) {os << NUMBER[frac[ix]];}
1107           else if (iRow * (row[ix] - iRow) > 0) {os << NUMBER[10];}
1108           else {os << " ";}
1109         } os << "\n";
1110       }
1111     } os << "\n";
1112
1113     // Print sign and value of bin contents
1114     double maxim = log10(max(yMax, -yMin));
1115     int iPowBin = int(floor(maxim + 0.0001));
1116     os << "          Contents  ";
1117     for (int ix = 0; ix < h.nBin ; ++ix) {
1118       if (h.res[ix] < - pow(10., iPowBin-4)) {os << "-";}
1119       else {os << " ";}
1120       row[ix] = int(abs(h.res[ix]) * pow(10.,3-iPowBin) + 0.5);
1121     } os << "\n";
1122     for (int iRow = 3; iRow >= 0; iRow--) {
1123       os << "            *10^" << setw(2) << iPowBin+iRow-3 << "  ";
1124       int mask = int( pow(10., iRow) + 0.5);
1125       for (int ix = 0; ix < h.nBin ; ++ix) {
1126         os << NUMBER[(row[ix] / mask) % 10];
1127       } os << "\n";
1128     } os << "\n";
1129
1130     // Print sign and value of lower bin edge.
1131     maxim = log10(max(-h.xMin, h.xMax - h.dx));
1132     int iPowExp = int(floor(maxim + 0.0001));
1133     os << "          Low edge  ";
1134     for (int ix = 0; ix < h.nBin ; ++ix) {
1135       if (h.xMin + ix * h.dx < - pow(10., iPowExp-3)) {os << "-";}
1136       else {os << " ";}
1137       row[ix] = int(abs(h.xMin + ix * h.dx) * pow(10.,2-iPowExp) + 0.5);
1138     } os << "\n";
1139     for (int iRow = 2; iRow >= 0; iRow--) {
1140       os << "            *10^" << setw(2) << iPowExp+iRow-2 << "  ";
1141       int mask = int( pow(10., iRow) + 0.5);
1142       for (int ix = 0; ix < h.nBin ; ++ix) {
1143         os << NUMBER[(row[ix] / mask) % 10];
1144       } os << "\n";
1145     } os << "\n";
1146   }
1147
1148   // Calculate and print statistics.
1149   double cSum = 0.;
1150   double cxSum = 0.;
1151   double cxxSum = 0.;
1152   for (int ix = 0; ix < h.nBin ; ++ix) {
1153     double cta = abs(h.res[ix]);
1154     double x = h.xMin + (ix + 0.5) * h.dx;
1155     cSum = cSum + cta;
1156     cxSum = cxSum + cta * x;
1157     cxxSum = cxxSum + cta * x * x;
1158   }
1159   double xmean = cxSum / max(cSum, Hist::TINY);
1160   double rms = sqrtpos( cxxSum / max(cSum, Hist::TINY) - xmean*xmean );
1161   os << scientific << setprecision(4)
1162      << "   Entries  =" << setw(12) << h.nFill
1163      << "    Mean =" << setw(12) << xmean
1164      << "    Underflow =" << setw(12) << h.under
1165      << "    Low edge  =" << setw(12) << h.xMin << "\n"
1166      << "   All chan =" << setw(12) << h.inside
1167      << "    Rms  =" << setw(12) << rms
1168      << "    Overflow  =" << setw(12) << h.over
1169      << "    High edge =" << setw(12) << h.xMax << endl;
1170   return os;
1171 }
1172
1173 //**************************************************************************
1174
1175 } // end namespace Pythia8