StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Basics.cc
1 // Basics.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 "Pythia8/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 // Distance in cylindrical (y, phi) coordinates.
538 
539 double RRapPhi(const Vec4& v1, const Vec4& v2) {
540  double dRap = abs(v1.rap() - v2.rap());
541  double dPhi = abs(v1.phi() - v2.phi());
542  if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
543  return sqrt(dRap*dRap + dPhi*dPhi);
544 }
545 
546 //--------------------------------------------------------------------------
547 
548 // Distance in cylindrical (eta, phi) coordinates.
549 
550 double REtaPhi(const Vec4& v1, const Vec4& v2) {
551  double dEta = abs(v1.eta() - v2.eta());
552  double dPhi = abs(v1.phi() - v2.phi());
553  if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
554  return sqrt(dEta*dEta + dPhi*dPhi);
555 }
556 
557 //--------------------------------------------------------------------------
558 
559 // Print a four-vector: also operator overloading with friend.
560 
561 ostream& operator<<(ostream& os, const Vec4& v) {
562  os << fixed << setprecision(3) << " " << setw(9) << v.xx << " "
563  << setw(9) << v.yy << " " << setw(9) << v.zz << " " << setw(9)
564  << v.tt << " (" << setw(9) << v.mCalc() << ")\n";
565  return os;
566 }
567 
568 //==========================================================================
569 
570 // RotBstMatrix class.
571 // This class implements 4 * 4 matrices that encode an arbitrary combination
572 // of rotations and boosts, that can be applied to Vec4 four-vectors.
573 
574 //--------------------------------------------------------------------------
575 
576 // Constants: could be changed here if desired, but normally should not.
577 // These are of technical nature, as described for each.
578 
579 // Small number to avoid division by zero.
580 const double RotBstMatrix::TINY = 1e-20;
581 
582 //--------------------------------------------------------------------------
583 
584 // Rotate by polar angle theta and azimuthal angle phi.
585 
586 void RotBstMatrix::rot(double theta, double phi) {
587 
588  // Set up rotation matrix.
589  double cthe = cos(theta);
590  double sthe = sin(theta);
591  double cphi = cos(phi);
592  double sphi = sin(phi);
593  double Mrot[4][4] = {
594  {1., 0., 0., 0.},
595  {0., cthe * cphi, - sphi, sthe * cphi},
596  {0., cthe * sphi, cphi, sthe * sphi},
597  {0., -sthe, 0., cthe } };
598 
599  // Rotate current matrix accordingly.
600  double Mtmp[4][4];
601  for (int i = 0; i < 4; ++i)
602  for (int j = 0; j < 4; ++j)
603  Mtmp[i][j] = M[i][j];
604  for (int i = 0; i < 4; ++i)
605  for (int j = 0; j < 4; ++j)
606  M[i][j] = Mrot[i][0] * Mtmp[0][j] + Mrot[i][1] * Mtmp[1][j]
607  + Mrot[i][2] * Mtmp[2][j] + Mrot[i][3] * Mtmp[3][j];
608 
609 }
610 
611 //--------------------------------------------------------------------------
612 
613 // Rotate so that vector originally along z axis becomes parallel with p.
614 
615 void RotBstMatrix::rot(const Vec4& p) {
616 
617  double theta = p.theta();
618  double phi = p.phi();
619  rot(0., -phi);
620  rot(theta, phi);
621 
622 }
623 
624 //--------------------------------------------------------------------------
625 
626 // Boost with velocity vector (betaX, betaY, betaZ).
627 
628 void RotBstMatrix::bst(double betaX, double betaY, double betaZ) {
629 
630  // Set up boost matrix.
631  double gm = 1. / sqrt( max( TINY, 1. - betaX*betaX - betaY*betaY
632  - betaZ*betaZ ) );
633  double gf = gm*gm / (1. + gm);
634  double Mbst[4][4] = {
635  { gm, gm*betaX, gm*betaY, gm*betaZ },
636  { gm*betaX, 1. + gf*betaX*betaX, gf*betaX*betaY, gf*betaX*betaZ },
637  { gm*betaY, gf*betaY*betaX, 1. + gf*betaY*betaY, gf*betaY*betaZ },
638  { gm*betaZ, gf*betaZ*betaX, gf*betaZ*betaY, 1. + gf*betaZ*betaZ } };
639 
640  // Boost current matrix correspondingly.
641  double Mtmp[4][4];
642  for (int i = 0; i < 4; ++i)
643  for (int j = 0; j < 4; ++j)
644  Mtmp[i][j] = M[i][j];
645  for (int i = 0; i < 4; ++i)
646  for (int j = 0; j < 4; ++j)
647  M[i][j] = Mbst[i][0] * Mtmp[0][j] + Mbst[i][1] * Mtmp[1][j]
648  + Mbst[i][2] * Mtmp[2][j] + Mbst[i][3] * Mtmp[3][j];
649 
650 }
651 
652 //--------------------------------------------------------------------------
653 
654 // Boost so that vector originally at rest obtains same velocity as p.
655 
656 void RotBstMatrix::bst(const Vec4& p) {
657  double betaX = p.px() / p.e();
658  double betaY = p.py() / p.e();
659  double betaZ = p.pz() / p.e();
660  bst(betaX, betaY, betaZ);
661 }
662 
663 //--------------------------------------------------------------------------
664 
665 // Boost so vector originally with same velocity as p is brought to rest.
666 
667 void RotBstMatrix::bstback(const Vec4& p) {
668  double betaX = -p.px() / p.e();
669  double betaY = -p.py() / p.e();
670  double betaZ = -p.pz() / p.e();
671  bst(betaX, betaY, betaZ);
672 }
673 
674 //--------------------------------------------------------------------------
675 
676 // Boost that transforms p1 to p2, where p1^2 = p2^2 is assumed.
677 
678 void RotBstMatrix::bst(const Vec4& p1, const Vec4& p2) {
679  double eSum = p1.e() + p2.e();
680  double betaX = (p2.px() - p1.px()) / eSum;
681  double betaY = (p2.py() - p1.py()) / eSum;
682  double betaZ = (p2.pz() - p1.pz()) / eSum;
683  double fac = 2. / (1. + betaX*betaX + betaY*betaY + betaZ*betaZ);
684  betaX *= fac; betaY *= fac; betaZ *= fac;
685  bst(betaX, betaY, betaZ);
686 }
687 
688 //--------------------------------------------------------------------------
689 
690 // Boost and rotation that transforms from p1 and p2
691 // to their rest frame with p1 along +z axis.
692 
693 void RotBstMatrix::toCMframe(const Vec4& p1, const Vec4& p2) {
694  Vec4 pSum = p1 + p2;
695  Vec4 dir = p1;
696  dir.bstback(pSum);
697  double theta = dir.theta();
698  double phi = dir.phi();
699  bstback(pSum);
700  rot(0., -phi);
701  rot(-theta, phi);
702 }
703 
704 //--------------------------------------------------------------------------
705 
706 // Rotation and boost that transforms from rest frame of p1 and p2
707 // with p1 along +z axis to actual frame of p1 and p2. (Inverse of above.)
708 
709 void RotBstMatrix::fromCMframe(const Vec4& p1, const Vec4& p2) {
710  Vec4 pSum = p1 + p2;
711  Vec4 dir = p1;
712  dir.bstback(pSum);
713  double theta = dir.theta();
714  double phi = dir.phi();
715  rot(0., -phi);
716  rot(theta, phi);
717  bst(pSum);
718 }
719 
720 //--------------------------------------------------------------------------
721 
722 // Combine existing rotation/boost matrix with another one.
723 
724 void RotBstMatrix::rotbst(const RotBstMatrix& Mrb) {
725  double Mtmp[4][4];
726  for (int i = 0; i < 4; ++i)
727  for (int j = 0; j < 4; ++j)
728  Mtmp[i][j] = M[i][j];
729  for (int i = 0; i < 4; ++i)
730  for (int j = 0; j < 4; ++j)
731  M[i][j] = Mrb.M[i][0] * Mtmp[0][j] + Mrb.M[i][1] * Mtmp[1][j]
732  + Mrb.M[i][2] * Mtmp[2][j] + Mrb.M[i][3] * Mtmp[3][j];
733 }
734 
735 //--------------------------------------------------------------------------
736 
737 // Invert the rotation and boost.
738 
739 void RotBstMatrix::invert() {
740  double Mtmp[4][4];
741  for (int i = 0; i < 4; ++i)
742  for (int j = 0; j < 4; ++j)
743  Mtmp[i][j] = M[i][j];
744  for (int i = 0; i < 4; ++i)
745  for (int j = 0; j < 4; ++j)
746  M[i][j] = ( (i == 0 && j > 0) || (i > 0 && j == 0) )
747  ? - Mtmp[j][i] : Mtmp[j][i];
748 }
749 
750 //--------------------------------------------------------------------------
751 
752 // Reset to diagonal matrix.
753 
754 void RotBstMatrix::reset() {
755  for (int i = 0; i < 4; ++i)
756  for (int j = 0; j < 4; ++j)
757  M[i][j] = (i==j) ? 1. : 0.;
758 }
759 
760 //--------------------------------------------------------------------------
761 
762 // Crude estimate deviation from unit matrix.
763 
764 double RotBstMatrix::deviation() const {
765  double devSum = 0.;
766  for (int i = 0; i < 4; ++i)
767  for (int j = 0; j < 4; ++j)
768  devSum += (i==j) ? abs(M[i][j] - 1.) : abs(M[i][j]);
769  return devSum;
770 }
771 
772 //--------------------------------------------------------------------------
773 
774 // Print a rotation and boost matrix: operator overloading with friend.
775 
776 ostream& operator<<(ostream& os, const RotBstMatrix& M) {
777  os << fixed << setprecision(5) << " Rotation/boost matrix: \n";
778  for (int i = 0; i <4; ++i)
779  os << setw(10) << M.M[i][0] << setw(10) << M.M[i][1]
780  << setw(10) << M.M[i][2] << setw(10) << M.M[i][3] << "\n";
781  return os;
782 }
783 
784 //==========================================================================
785 
786 // Hist class.
787 // This class handles a single histogram at a time
788 // (or a vector of histograms).
789 
790 //--------------------------------------------------------------------------
791 
792 // Constants: could be changed here if desired, but normally should not.
793 // These are of technical nature, as described for each.
794 
795 // Maximum number of bins in a histogram.
796 const int Hist::NBINMAX = 1000;
797 
798 // Maximum number of columns that can be printed for a histogram.
799 const int Hist::NCOLMAX = 100;
800 
801 // Maximum number of lines a histogram can use at output.
802 const int Hist::NLINES = 30;
803 
804 // Tolerance in deviation of xMin and xMax between two histograms.
805 const double Hist::TOLERANCE = 0.001;
806 
807 // Small and large numbers to avoid division by zero and overflow.
808 const double Hist::TINY = 1e-20;
809 const double Hist::LARGE = 1e20;
810 
811 // When minbin/maxbin < SMALLFRAC the y scale goes down to zero.
812 const double Hist::SMALLFRAC = 0.1;
813 
814 // Constants for printout: fixed steps on y scale; filling characters.
815 const double DYAC[] = {0.04, 0.05, 0.06, 0.08, 0.10,
816  0.12, 0.15, 0.20, 0.25, 0.30};
817 const char NUMBER[] = {'0', '1', '2', '3', '4', '5',
818  '6', '7', '8', '9', 'X' };
819 
820 //--------------------------------------------------------------------------
821 
822 // Book a histogram.
823 
824 void Hist::book(string titleIn, int nBinIn, double xMinIn,
825  double xMaxIn) {
826 
827  title = titleIn;
828  nBin = nBinIn;
829  if (nBinIn < 1) nBin = 1;
830  if (nBinIn > NBINMAX) nBin = NBINMAX;
831  xMin = xMinIn;
832  xMax = xMaxIn;
833  dx = (xMax - xMin)/nBin;
834  res.resize(nBin);
835  null();
836 
837 }
838 
839 //--------------------------------------------------------------------------
840 
841 // Reset bin contents.
842 
843 void Hist::null() {
844 
845  nFill = 0;
846  under = 0.;
847  inside = 0.;
848  over = 0.;
849  for (int ix = 0; ix < nBin; ++ix) res[ix] = 0.;
850 
851 }
852 
853 //--------------------------------------------------------------------------
854 
855 // Fill bin with weight.
856 
857 void Hist::fill(double x, double w) {
858 
859  ++nFill;
860  int iBin = int(floor((x - xMin)/dx));
861  if (iBin < 0) under += w;
862  else if (iBin >= nBin) over += w;
863  else {inside += w; res[iBin] += w; }
864 
865 }
866 
867 //--------------------------------------------------------------------------
868 
869 // Print a histogram: also operator overloading with friend.
870 
871 ostream& operator<<(ostream& os, const Hist& h) {
872 
873  // Do not print empty histograms.
874  if (h.nFill <= 0) return os;
875 
876  // Write time and title.
877  time_t t = time(0);
878  char date[18];
879  strftime(date,18,"%Y-%m-%d %H:%M",localtime(&t));
880  os << "\n\n " << date << " " << h.title << "\n\n";
881 
882  // Group bins, where required, to make printout have fewer columns.
883  // Avoid overflow.
884  int nGroup = 1 + (h.nBin - 1) / Hist::NCOLMAX;
885  int nCol = 1 + (h.nBin - 1) / nGroup;
886  vector<double> resCol(nCol);
887  for (int iCol = 0; iCol < nCol; ++iCol) {
888  resCol[iCol] = 0.;
889  for (int ix = nGroup * iCol; ix < min( h.nBin, nGroup * (iCol + 1)); ++ix)
890  resCol[iCol] += h.res[ix];
891  resCol[iCol] = max( -Hist::LARGE, min( Hist::LARGE, resCol[iCol] ) );
892  }
893 
894  // Find minimum and maximum bin content.
895  double yMin = resCol[0];
896  double yMax = resCol[0];
897  for (int iCol = 1; iCol < nCol; ++iCol) {
898  if (resCol[iCol] < yMin) yMin = resCol[iCol];
899  if (resCol[iCol] > yMax) yMax = resCol[iCol];
900  }
901 
902  // Determine scale and step size for y axis.
903  if (yMax - yMin > Hist::NLINES * DYAC[0] * 1e-9) {
904  if (yMin > 0. && yMin < Hist::SMALLFRAC * yMax) yMin = 0.;
905  if (yMax < 0. && yMax > Hist::SMALLFRAC * yMin) yMax = 0.;
906  int iPowY = int(floor( log10(yMax - yMin) ));
907  if (yMax - yMin < Hist::NLINES * DYAC[0] * pow(10.,iPowY))
908  iPowY = iPowY - 1;
909  if (yMax - yMin > Hist::NLINES * DYAC[9] * pow(10.,iPowY))
910  iPowY = iPowY + 1;
911  double nLinePow = Hist::NLINES * pow(10.,iPowY);
912  double delY = DYAC[0];
913  for (int idel = 0; idel < 9; ++idel)
914  if (yMax - yMin >= nLinePow * DYAC[idel]) delY = DYAC[idel+1];
915  double dy = delY * pow(10.,iPowY);
916 
917  // Convert bin contents to integer form; fractional fill in top row.
918  vector<int> row(nCol);
919  vector<int> frac(nCol);
920  for (int iCol = 0; iCol < nCol ; ++iCol) {
921  double cta = abs(resCol[iCol]) / dy;
922  row[iCol] = int(cta + 0.95);
923  if(resCol[iCol] < 0.) row[iCol] = - row[iCol];
924  frac[iCol] = int(10. * (cta + 1.05 - floor(cta + 0.95)));
925  }
926  int rowMin = int(abs(yMin)/dy + 0.95);
927  if ( yMin < 0) rowMin = - rowMin;
928  int rowMax = int(abs(yMax)/dy + 0.95);
929  if ( yMax < 0) rowMax = - rowMax;
930 
931  // Print histogram row by row.
932  os << fixed << setprecision(2);
933  for (int iRow = rowMax; iRow >= rowMin; iRow--) if (iRow != 0) {
934  os << " " << setw(10) << iRow*delY << "*10^"
935  << setw(2) << iPowY << " ";
936  for (int iCol = 0; iCol < nCol ; ++iCol) {
937  if (iRow == row[iCol]) os << NUMBER[frac[iCol]];
938  else if (iRow * (row[iCol] - iRow) > 0) os << NUMBER[10];
939  else os << " ";
940  } os << "\n";
941  } os << "\n";
942 
943  // Print sign and value of bin contents
944  double maxim = log10(max(yMax, -yMin));
945  int iPowBin = int(floor(maxim + 0.0001));
946  os << " Contents ";
947  for (int iCol = 0; iCol < nCol ; ++iCol) {
948  if (resCol[iCol] < - pow(10., iPowBin - 4)) os << "-";
949  else os << " ";
950  row[iCol] = int(abs(resCol[iCol]) * pow(10., 3 - iPowBin) + 0.5);
951  } os << "\n";
952  for (int iRow = 3; iRow >= 0; iRow--) {
953  os << " *10^" << setw(2) << iPowBin + iRow - 3 << " ";
954  int mask = int( pow(10., iRow) + 0.5);
955  for (int iCol = 0; iCol < nCol ; ++iCol) {
956  os << NUMBER[(row[iCol] / mask) % 10];
957  } os << "\n";
958  } os << "\n";
959 
960  // Print sign and value of lower bin edge.
961  maxim = log10( max( -h.xMin, h.xMax - h.dx));
962  int iPowExp = int(floor(maxim + 0.0001));
963  os << " Low edge ";
964  for (int iCol = 0; iCol < nCol ; ++iCol) {
965  if (h.xMin + iCol * nGroup * h.dx < - pow(10., iPowExp - 3)) os << "-";
966  else os << " ";
967  row[iCol] = int(abs(h.xMin + iCol * nGroup * h.dx)
968  * pow(10., 2 - iPowExp) + 0.5);
969  } os << "\n";
970  for (int iRow = 2; iRow >= 0; iRow--) {
971  os << " *10^" << setw(2) << iPowExp + iRow - 2 << " ";
972  int mask = int( pow(10., iRow) + 0.5);
973  for (int iCol = 0; iCol < nCol ; ++iCol)
974  os << NUMBER[(row[iCol] / mask) % 10];
975  os << "\n";
976  } os << "\n";
977 
978  // Print explanation if histogram cannot be shown.
979  } else os << " Histogram not shown since lowest value" << scientific
980  << setprecision(4) << setw(12) << yMin << " and highest value"
981  << setw(12) << yMax << " are too close \n \n";
982 
983  // Calculate and print statistics.
984  double cSum = 0.;
985  double cxSum = 0.;
986  double cxxSum = 0.;
987  for (int ix = 0; ix < h.nBin ; ++ix) {
988  double cta = abs(h.res[ix]);
989  double x = h.xMin + (ix + 0.5) * h.dx;
990  cSum = cSum + cta;
991  cxSum = cxSum + cta * x;
992  cxxSum = cxxSum + cta * x * x;
993  }
994  double xmean = cxSum / max(cSum, Hist::TINY);
995  double rms = sqrtpos( cxxSum / max(cSum, Hist::TINY) - xmean*xmean );
996  os << scientific << setprecision(4)
997  << " Entries =" << setw(12) << h.nFill
998  << " Mean =" << setw(12) << xmean
999  << " Underflow =" << setw(12) << h.under
1000  << " Low edge =" << setw(12) << h.xMin << "\n"
1001  << " All chan =" << setw(12) << h.inside
1002  << " Rms =" << setw(12) << rms
1003  << " Overflow =" << setw(12) << h.over
1004  << " High edge =" << setw(12) << h.xMax << endl;
1005  return os;
1006 }
1007 
1008 //--------------------------------------------------------------------------
1009 
1010 // Print histogram contents as a table (e.g. for Gnuplot).
1011 
1012 void Hist::table(ostream& os, bool printOverUnder, bool xMidBin) const {
1013 
1014  // Print histogram vector bin by bin, with mean x as first column.
1015  os << scientific << setprecision(4);
1016  double xBeg = (xMidBin) ? xMin + 0.5 * dx : xMin;
1017  if (printOverUnder)
1018  os << setw(12) << xBeg - dx << setw(12) << under << "\n";
1019  for (int ix = 0; ix < nBin; ++ix)
1020  os << setw(12) << xBeg + ix * dx << setw(12) << res[ix] << "\n";
1021  if (printOverUnder)
1022  os << setw(12) << xBeg + nBin * dx << setw(12) << over << "\n";
1023 
1024 }
1025 
1026 //--------------------------------------------------------------------------
1027 
1028 // Print a table out of two histograms with same x axis (e.g. for Gnuplot).
1029 
1030 void table(const Hist& h1, const Hist& h2, ostream& os, bool printOverUnder,
1031  bool xMidBin) {
1032 
1033  // Require histogram x axes to agree.
1034  if (h1.nBin != h2.nBin || abs(h1.xMin - h2.xMin) > Hist::TOLERANCE * h1.dx
1035  || abs(h1.xMax - h2.xMax) > Hist::TOLERANCE * h1.dx) return;
1036 
1037  // Print histogram vectors bin by bin, with mean x as first column.
1038  os << scientific << setprecision(4);
1039  double xBeg = (xMidBin) ? h1.xMin + 0.5 * h1.dx : h1.xMin;
1040  if (printOverUnder)
1041  os << setw(12) << xBeg - h1.dx << setw(12) << h1.under
1042  << setw(12) << h2.under << "\n";
1043  for (int ix = 0; ix < h1.nBin; ++ix)
1044  os << setw(12) << xBeg + ix * h1.dx << setw(12) << h1.res[ix]
1045  << setw(12) << h2.res[ix] << "\n";
1046  if (printOverUnder)
1047  os << setw(12) << xBeg + h1.nBin * h1.dx << setw(12) << h1.over
1048  << setw(12) << h2.over << "\n";
1049 
1050 }
1051 
1052 void table(const Hist& h1, const Hist& h2, string fileName,
1053  bool printOverUnder, bool xMidBin) {
1054 
1055  ofstream streamName(fileName.c_str());
1056  table( h1, h2, streamName, printOverUnder, xMidBin);
1057 
1058 }
1059 
1060 //--------------------------------------------------------------------------
1061 
1062 // Get content of specific bin.
1063 // Special values are bin 0 for underflow and bin nBin+1 for overflow.
1064 // All other bins outside proper histogram range return 0.
1065 
1066 double Hist::getBinContent(int iBin) const {
1067 
1068  if (iBin > 0 && iBin <= nBin) return res[iBin - 1];
1069  else if (iBin == 0) return under;
1070  else if (iBin == nBin + 1) return over;
1071  else return 0.;
1072 
1073 }
1074 
1075 //--------------------------------------------------------------------------
1076 
1077 // Check whether another histogram has same size and limits.
1078 
1079 bool Hist::sameSize(const Hist& h) const {
1080 
1081  if (nBin == h.nBin && abs(xMin - h.xMin) < TOLERANCE * dx &&
1082  abs(xMax - h.xMax) < TOLERANCE * dx) return true;
1083  else return false;
1084 
1085 }
1086 
1087 //--------------------------------------------------------------------------
1088 
1089 // Take 10-logarithm or natural logarithm of contents bin by bin.
1090 
1091 void Hist::takeLog(bool tenLog) {
1092 
1093  // Find smallest positive bin content, and put min a bit below.
1094  double yMin = Hist::LARGE;
1095  for (int ix = 0; ix < nBin; ++ix)
1096  if (res[ix] > Hist::TINY && res[ix] < yMin ) yMin = res[ix];
1097  yMin *= 0.8;
1098 
1099  // Take 10-logarithm bin by bin, but ensure positivity.
1100  if (tenLog) {
1101  for (int ix = 0; ix < nBin; ++ix)
1102  res[ix] = log10( max( yMin, res[ix]) );
1103  under = log10( max( yMin, under) );
1104  inside = log10( max( yMin, inside) );
1105  over = log10( max( yMin, over) );
1106 
1107  // Take natural logarithm bin by bin, but ensure positivity.
1108  } else {
1109  for (int ix = 0; ix < nBin; ++ix)
1110  res[ix] = log( max( yMin, res[ix]) );
1111  under = log( max( yMin, under) );
1112  inside = log( max( yMin, inside) );
1113  over = log( max( yMin, over) );
1114  }
1115 
1116 }
1117 
1118 //--------------------------------------------------------------------------
1119 
1120 // Take square root of contents bin by bin; set 0 for negative content.
1121 
1122 void Hist::takeSqrt() {
1123 
1124  for (int ix = 0; ix < nBin; ++ix) res[ix] = sqrtpos(res[ix]);
1125  under = sqrtpos(under);
1126  inside = sqrtpos(inside);
1127  over = sqrtpos(over);
1128 
1129 }
1130 
1131 //--------------------------------------------------------------------------
1132 
1133 // Add histogram to existing one.
1134 
1135 Hist& Hist::operator+=(const Hist& h) {
1136  if (!sameSize(h)) return *this;
1137  nFill += h.nFill;
1138  under += h.under;
1139  inside += h.inside;
1140  over += h.over;
1141  for (int ix = 0; ix < nBin; ++ix) res[ix] += h.res[ix];
1142  return *this;
1143 }
1144 
1145 //--------------------------------------------------------------------------
1146 
1147 // Subtract histogram from existing one.
1148 
1149 Hist& Hist::operator-=(const Hist& h) {
1150  if (!sameSize(h)) return *this;
1151  nFill += h.nFill;
1152  under -= h.under;
1153  inside -= h.inside;
1154  over -= h.over;
1155  for (int ix = 0; ix < nBin; ++ix) res[ix] -= h.res[ix];
1156  return *this;
1157 }
1158 
1159 //--------------------------------------------------------------------------
1160 
1161 // Multiply existing histogram by another one.
1162 
1163 Hist& Hist::operator*=(const Hist& h) {
1164  if (!sameSize(h)) return *this;
1165  nFill += h.nFill;
1166  under *= h.under;
1167  inside *= h.inside;
1168  over *= h.over;
1169  for (int ix = 0; ix < nBin; ++ix) res[ix] *= h.res[ix];
1170  return *this;
1171 }
1172 
1173 //--------------------------------------------------------------------------
1174 
1175 // Divide existing histogram by another one.
1176 
1177 Hist& Hist::operator/=(const Hist& h) {
1178  if (!sameSize(h)) return *this;
1179  nFill += h.nFill;
1180  under = (abs(h.under) < Hist::TINY) ? 0. : under/h.under;
1181  inside = (abs(h.inside) < Hist::TINY) ? 0. : inside/h.inside;
1182  over = (abs(h.over) < Hist::TINY) ? 0. : over/h.over;
1183  for (int ix = 0; ix < nBin; ++ix)
1184  res[ix] = (abs(h.res[ix]) < Hist::TINY) ? 0. : res[ix]/h.res[ix];
1185  return *this;
1186 }
1187 
1188 //--------------------------------------------------------------------------
1189 
1190 // Add constant offset to histogram.
1191 
1192 Hist& Hist::operator+=(double f) {
1193  under += f;
1194  inside += nBin * f;
1195  over += f;
1196  for (int ix = 0; ix < nBin; ++ix) res[ix] += f;
1197  return *this;
1198 }
1199 
1200 //--------------------------------------------------------------------------
1201 
1202 // Subtract constant offset from histogram.
1203 
1204 Hist& Hist::operator-=(double f) {
1205  under -= f;
1206  inside -= nBin * f;
1207  over -= f;
1208  for (int ix = 0; ix < nBin; ++ix) res[ix] -= f;
1209  return *this;
1210 }
1211 
1212 //--------------------------------------------------------------------------
1213 
1214 // Multiply histogram by constant.
1215 
1216 Hist& Hist::operator*=(double f) {
1217  under *= f;
1218  inside *= f;
1219  over *= f;
1220  for (int ix = 0; ix < nBin; ++ix) res[ix] *= f;
1221  return *this;
1222 }
1223 
1224 //--------------------------------------------------------------------------
1225 
1226 // Divide histogram by constant.
1227 
1228 Hist& Hist::operator/=(double f) {
1229  if (abs(f) > Hist::TINY) {
1230  under /= f;
1231  inside /= f;
1232  over /= f;
1233  for (int ix = 0; ix < nBin; ++ix) res[ix] /= f;
1234  // Set empty contents when division by zero.
1235  } else {
1236  under = 0.;
1237  inside = 0.;
1238  over = 0.;
1239  for (int ix = 0; ix < nBin; ++ix) res[ix] = 0.;
1240  }
1241  return *this;
1242 }
1243 
1244 //--------------------------------------------------------------------------
1245 
1246 // Implementation of operator overloading with friends.
1247 
1248 Hist operator+(double f, const Hist& h1) {
1249  Hist h = h1; return h += f;}
1250 
1251 Hist operator+(const Hist& h1, double f) {
1252  Hist h = h1; return h += f;}
1253 
1254 Hist operator+(const Hist& h1, const Hist& h2) {
1255  Hist h = h1; return h += h2;}
1256 
1257 Hist operator-(double f, const Hist& h1) {
1258  Hist h = h1;
1259  h.under = f - h1.under;
1260  h.inside = h1.nBin * f - h1.inside;
1261  h.over = f - h1.over;
1262  for (int ix = 0; ix < h1.nBin; ++ix) h.res[ix] = f - h1.res[ix];
1263  return h;}
1264 
1265 Hist operator-(const Hist& h1, double f) {
1266  Hist h = h1; return h -= f;}
1267 
1268 Hist operator-(const Hist& h1, const Hist& h2) {
1269  Hist h = h1; return h -= h2;}
1270 
1271 Hist operator*(double f, const Hist& h1) {
1272  Hist h = h1; return h *= f;}
1273 
1274 Hist operator*(const Hist& h1, double f) {
1275  Hist h = h1; return h *= f;}
1276 
1277 Hist operator*(const Hist& h1, const Hist& h2) {
1278  Hist h = h1; return h *= h2;}
1279 
1280 Hist operator/(double f, const Hist& h1) {
1281  Hist h = h1;
1282  h.under = (abs(h1.under) < Hist::TINY) ? 0. : f/h1.under;
1283  h.inside = (abs(h1.inside) < Hist::TINY) ? 0. : f/h1.inside;
1284  h.over = (abs(h1.over) < Hist::TINY) ? 0. : f/h1.over;
1285  for (int ix = 0; ix < h1.nBin; ++ix)
1286  h.res[ix] = (abs(h1.res[ix]) < Hist::TINY) ? 0. : f/h1.res[ix];
1287  return h;
1288 }
1289 
1290 Hist operator/(const Hist& h1, double f) {
1291  Hist h = h1; return h /= f;}
1292 
1293 Hist operator/(const Hist& h1, const Hist& h2) {
1294  Hist h = h1; return h /= h2;}
1295 
1296 //==========================================================================
1297 
1298 } // end namespace Pythia8
Generated by   doxygen 1.8.5