PYTHIA8/pythia8170/src/Basics.cxx

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