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) 2018 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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 #include <limits>
14 
15 namespace Pythia8 {
16 
17 //==========================================================================
18 
19 // Rndm class.
20 // This class handles random number generation according to the
21 // Marsaglia-Zaman-Tsang algorithm
22 
23 //--------------------------------------------------------------------------
24 
25 // Constants: could be changed here if desired, but normally should not.
26 // These are of technical nature, as described for each.
27 
28 // The default seed, i.e. the Marsaglia-Zaman random number sequence.
29 const int Rndm::DEFAULTSEED = 19780503;
30 
31 //--------------------------------------------------------------------------
32 
33 // Method to pass in pointer for external random number generation.
34 
35 bool Rndm::rndmEnginePtr( RndmEngine* rndmEngPtrIn) {
36 
37  // Save pointer.
38  if (rndmEngPtrIn == 0) return false;
39  rndmEngPtr = rndmEngPtrIn;
40  useExternalRndm = true;
41 
42  // Done.
43  return true;
44 
45 }
46 
47 //--------------------------------------------------------------------------
48 
49 // Initialize, normally at construction or in first call.
50 
51 void Rndm::init(int seedIn) {
52 
53  // Pick seed in convenient way. Assure it to be non-negative.
54  int seed = seedIn;
55  if (seedIn < 0) seed = DEFAULTSEED;
56  else if (seedIn == 0) seed = int(time(0));
57  if (seed < 0) seed = -seed;
58 
59  // Unpack seed.
60  int ij = (seed/30082) % 31329;
61  int kl = seed % 30082;
62  int i = (ij/177) % 177 + 2;
63  int j = ij % 177 + 2;
64  int k = (kl/169) % 178 + 1;
65  int l = kl % 169;
66 
67  // Initialize random number array.
68  for (int ii = 0; ii < 97; ++ii) {
69  double s = 0.;
70  double t = 0.5;
71  for (int jj = 0; jj < 48; ++jj) {
72  int m = (( (i*j)%179 )*k) % 179;
73  i = j;
74  j = k;
75  k = m;
76  l = (53*l+1) % 169;
77  if ( (l*m) % 64 >= 32) s += t;
78  t *= 0.5;
79  }
80  u[ii] = s;
81  }
82 
83  // Initialize other variables.
84  double twom24 = 1.;
85  for (int i24 = 0; i24 < 24; ++i24) twom24 *= 0.5;
86  c = 362436. * twom24;
87  cd = 7654321. * twom24;
88  cm = 16777213. * twom24;
89  i97 = 96;
90  j97 = 32;
91 
92  // Finished.
93  initRndm = true;
94  seedSave = seed;
95  sequence = 0;
96 
97 }
98 
99 //--------------------------------------------------------------------------
100 
101 // Generate next random number uniformly between 0 and 1.
102 
103 double Rndm::flat() {
104 
105  // Use external random number generator if such has been linked.
106  if (useExternalRndm) return rndmEngPtr->flat();
107 
108  // Ensure that already initialized.
109  if (!initRndm) init(DEFAULTSEED);
110 
111  // Find next random number and update saved state.
112  ++sequence;
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 // Pick one option among vector of (positive) probabilities.
132 
133 int Rndm::pick(const vector<double>& prob) {
134 
135  double work = 0.;
136  for (int i = 0; i < int(prob.size()); ++i) work += prob[i];
137  work *= flat();
138  int index = -1;
139  do work -= prob[++index];
140  while (work > 0. && index < int(prob.size()));
141  return index;
142 
143 }
144 
145 //--------------------------------------------------------------------------
146 
147 // Save current state of the random number generator to a binary file.
148 
149 bool Rndm::dumpState(string fileName) {
150 
151  // Open file as output stream.
152  const char* fn = fileName.c_str();
153  ofstream ofs(fn, ios::binary);
154 
155  if (!ofs.good()) {
156  cout << " Rndm::dumpState: could not open output file" << endl;
157  return false;
158  }
159 
160  // Write the state of the generator on the file.
161  ofs.write((char *) &seedSave, sizeof(int));
162  ofs.write((char *) &sequence, sizeof(long));
163  ofs.write((char *) &i97, sizeof(int));
164  ofs.write((char *) &j97, sizeof(int));
165  ofs.write((char *) &c, sizeof(double));
166  ofs.write((char *) &cd, sizeof(double));
167  ofs.write((char *) &cm, sizeof(double));
168  ofs.write((char *) &u, sizeof(double) * 97);
169 
170  // Write confirmation on cout.
171  cout << " PYTHIA Rndm::dumpState: seed = " << seedSave
172  << ", sequence no = " << sequence << endl;
173  return true;
174 
175 }
176 
177 //--------------------------------------------------------------------------
178 
179 // Read in the state of the random number generator from a binary file.
180 
181 bool Rndm::readState(string fileName) {
182 
183  // Open file as input stream.
184  const char* fn = fileName.c_str();
185  ifstream ifs(fn, ios::binary);
186 
187  if (!ifs.good()) {
188  cout << " Rndm::readState: could not open input file" << endl;
189  return false;
190  }
191 
192  // Read the state of the generator from the file.
193  ifs.read((char *) &seedSave, sizeof(int));
194  ifs.read((char *) &sequence, sizeof(long));
195  ifs.read((char *) &i97, sizeof(int));
196  ifs.read((char *) &j97, sizeof(int));
197  ifs.read((char *) &c, sizeof(double));
198  ifs.read((char *) &cd, sizeof(double));
199  ifs.read((char *) &cm, sizeof(double));
200  ifs.read((char *) &u, sizeof(double) *97);
201 
202  // Write confirmation on cout.
203  cout << " PYTHIA Rndm::readState: seed " << seedSave
204  << ", sequence no = " << sequence << endl;
205  return true;
206 
207 }
208 
209 //==========================================================================
210 
211 // Vec4 class.
212 // This class implements four-vectors, in energy-momentum space.
213 // (But could also be used to hold space-time four-vectors.)
214 
215 //--------------------------------------------------------------------------
216 
217 // Constants: could be changed here if desired, but normally should not.
218 // These are of technical nature, as described for each.
219 
220 // Small number to avoid division by zero.
221 const double Vec4::TINY = 1e-20;
222 
223 //--------------------------------------------------------------------------
224 
225 // Rotation (simple).
226 
227 void Vec4::rot(double thetaIn, double phiIn) {
228 
229  double cthe = cos(thetaIn);
230  double sthe = sin(thetaIn);
231  double cphi = cos(phiIn);
232  double sphi = sin(phiIn);
233  double tmpx = cthe * cphi * xx - sphi * yy + sthe * cphi * zz;
234  double tmpy = cthe * sphi * xx + cphi * yy + sthe * sphi * zz;
235  double tmpz = -sthe * xx + cthe * zz;
236  xx = tmpx;
237  yy = tmpy;
238  zz = tmpz;
239 
240 }
241 
242 //--------------------------------------------------------------------------
243 
244 // Azimuthal rotation phi around an arbitrary axis (nz, ny, nz).
245 
246 void Vec4::rotaxis(double phiIn, double nx, double ny, double nz) {
247 
248  double norm = 1./sqrt(nx*nx + ny*ny + nz*nz);
249  nx *= norm;
250  ny *= norm;
251  nz *= norm;
252  double cphi = cos(phiIn);
253  double sphi = sin(phiIn);
254  double comb = (nx * xx + ny * yy + nz * zz) * (1. - cphi);
255  double tmpx = cphi * xx + comb * nx + sphi * (ny * zz - nz * yy);
256  double tmpy = cphi * yy + comb * ny + sphi * (nz * xx - nx * zz);
257  double tmpz = cphi * zz + comb * nz + sphi * (nx * yy - ny * xx);
258  xx = tmpx;
259  yy = tmpy;
260  zz = tmpz;
261 
262 }
263 
264 //--------------------------------------------------------------------------
265 
266 // Azimuthal rotation phi around an arbitrary (3-vector component of) axis.
267 
268 void Vec4::rotaxis(double phiIn, const Vec4& n) {
269 
270  double nx = n.xx;
271  double ny = n.yy;
272  double nz = n.zz;
273  double norm = 1./sqrt(nx*nx + ny*ny + nz*nz);
274  nx *= norm;
275  ny *=norm;
276  nz *=norm;
277  double cphi = cos(phiIn);
278  double sphi = sin(phiIn);
279  double comb = (nx * xx + ny * yy + nz * zz) * (1. - cphi);
280  double tmpx = cphi * xx + comb * nx + sphi * (ny * zz - nz * yy);
281  double tmpy = cphi * yy + comb * ny + sphi * (nz * xx - nx * zz);
282  double tmpz = cphi * zz + comb * nz + sphi * (nx * yy - ny * xx);
283  xx = tmpx;
284  yy = tmpy;
285  zz = tmpz;
286 
287 }
288 
289 //--------------------------------------------------------------------------
290 
291 // Boost (simple).
292 
293 void Vec4::bst(double betaX, double betaY, double betaZ) {
294 
295  double beta2 = betaX*betaX + betaY*betaY + betaZ*betaZ;
296  if (beta2 >= 1.) return;
297  double gamma = 1. / sqrt(1. - beta2);
298  double prod1 = betaX * xx + betaY * yy + betaZ * zz;
299  double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
300  xx += prod2 * betaX;
301  yy += prod2 * betaY;
302  zz += prod2 * betaZ;
303  tt = gamma * (tt + prod1);
304 
305 }
306 
307 //--------------------------------------------------------------------------
308 
309 // Boost (simple, given gamma).
310 
311 void Vec4::bst(double betaX, double betaY, double betaZ, double gamma) {
312 
313  double prod1 = betaX * xx + betaY * yy + betaZ * zz;
314  double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
315  xx += prod2 * betaX;
316  yy += prod2 * betaY;
317  zz += prod2 * betaZ;
318  tt = gamma * (tt + prod1);
319 
320 }
321 
322 //--------------------------------------------------------------------------
323 
324 // Boost given by a Vec4 p.
325 
326 void Vec4::bst(const Vec4& pIn) {
327 
328  if (abs(pIn.tt) < Vec4::TINY) return;
329  double betaX = pIn.xx / pIn.tt;
330  double betaY = pIn.yy / pIn.tt;
331  double betaZ = pIn.zz / pIn.tt;
332  double beta2 = betaX*betaX + betaY*betaY + betaZ*betaZ;
333  if (beta2 >= 1.) return;
334  double gamma = 1. / sqrt(1. - beta2);
335  double prod1 = betaX * xx + betaY * yy + betaZ * zz;
336  double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
337  xx += prod2 * betaX;
338  yy += prod2 * betaY;
339  zz += prod2 * betaZ;
340  tt = gamma * (tt + prod1);
341 
342 }
343 
344 //--------------------------------------------------------------------------
345 
346 // Boost given by a Vec4 p and double m.
347 
348 void Vec4::bst(const Vec4& pIn, double mIn) {
349 
350  if (abs(pIn.tt) < Vec4::TINY) return;
351  double betaX = pIn.xx / pIn.tt;
352  double betaY = pIn.yy / pIn.tt;
353  double betaZ = pIn.zz / pIn.tt;
354  double gamma = pIn.tt / mIn;
355  double prod1 = betaX * xx + betaY * yy + betaZ * zz;
356  double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
357  xx += prod2 * betaX;
358  yy += prod2 * betaY;
359  zz += prod2 * betaZ;
360  tt = gamma * (tt + prod1);
361 
362 }
363 
364 //--------------------------------------------------------------------------
365 
366 // Boost given by a Vec4 p; boost in opposite direction.
367 
368 void Vec4::bstback(const Vec4& pIn) {
369 
370  if (abs(pIn.tt) < Vec4::TINY) return;
371  double betaX = -pIn.xx / pIn.tt;
372  double betaY = -pIn.yy / pIn.tt;
373  double betaZ = -pIn.zz / pIn.tt;
374  double beta2 = betaX*betaX + betaY*betaY + betaZ*betaZ;
375  if (beta2 >= 1.) return;
376  double gamma = 1. / sqrt(1. - beta2);
377  double prod1 = betaX * xx + betaY * yy + betaZ * zz;
378  double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
379  xx += prod2 * betaX;
380  yy += prod2 * betaY;
381  zz += prod2 * betaZ;
382  tt = gamma * (tt + prod1);
383 
384 }
385 
386 //--------------------------------------------------------------------------
387 
388 // Boost given by a Vec4 p and double m; boost in opposite direction.
389 
390 void Vec4::bstback(const Vec4& pIn, double mIn) {
391 
392  if (abs(pIn.tt) < Vec4::TINY) return;
393  double betaX = -pIn.xx / pIn.tt;
394  double betaY = -pIn.yy / pIn.tt;
395  double betaZ = -pIn.zz / pIn.tt;
396  double gamma = pIn.tt / mIn;
397  double prod1 = betaX * xx + betaY * yy + betaZ * zz;
398  double prod2 = gamma * (gamma * prod1 / (1. + gamma) + tt);
399  xx += prod2 * betaX;
400  yy += prod2 * betaY;
401  zz += prod2 * betaZ;
402  tt = gamma * (tt + prod1);
403 
404 }
405 
406 //--------------------------------------------------------------------------
407 
408 // Arbitrary combination of rotations and boosts defined by 4 * 4 matrix.
409 
410 void Vec4::rotbst(const RotBstMatrix& M) {
411 
412  double x = xx; double y = yy; double z = zz; double t = tt;
413  tt = M.M[0][0] * t + M.M[0][1] * x + M.M[0][2] * y + M.M[0][3] * z;
414  xx = M.M[1][0] * t + M.M[1][1] * x + M.M[1][2] * y + M.M[1][3] * z;
415  yy = M.M[2][0] * t + M.M[2][1] * x + M.M[2][2] * y + M.M[2][3] * z;
416  zz = M.M[3][0] * t + M.M[3][1] * x + M.M[3][2] * y + M.M[3][3] * z;
417 
418 }
419 
420 //--------------------------------------------------------------------------
421 
422 // Print a four-vector: also operator overloading with friend.
423 
424 ostream& operator<<(ostream& os, const Vec4& v) {
425  os << fixed << setprecision(3) << " " << setw(9) << v.xx << " "
426  << setw(9) << v.yy << " " << setw(9) << v.zz << " " << setw(9)
427  << v.tt << " (" << setw(9) << v.mCalc() << ")\n";
428  return os;
429 }
430 
431 //--------------------------------------------------------------------------
432 
433 // The invariant mass of two four-vectors.
434 
435 double m(const Vec4& v1, const Vec4& v2) {
436  double m2 = pow2(v1.tt + v2.tt) - pow2(v1.xx + v2.xx)
437  - pow2(v1.yy + v2.yy) - pow2(v1.zz + v2.zz);
438  return (m2 > 0.) ? sqrt(m2) : 0.;
439 }
440 
441 //--------------------------------------------------------------------------
442 
443 // The squared invariant mass of two four-vectors.
444 
445 double m2(const Vec4& v1, const Vec4& v2) {
446  double m2 = pow2(v1.tt + v2.tt) - pow2(v1.xx + v2.xx)
447  - pow2(v1.yy + v2.yy) - pow2(v1.zz + v2.zz);
448  return m2;
449 }
450 
451 //--------------------------------------------------------------------------
452 
453 // The scalar product of two three-vectors.
454 
455 double dot3(const Vec4& v1, const Vec4& v2) {
456  return v1.xx*v2.xx + v1.yy*v2.yy + v1.zz*v2.zz;
457 }
458 
459 //--------------------------------------------------------------------------
460 
461 // The cross product of two three-vectors.
462 
463 Vec4 cross3(const Vec4& v1, const Vec4& v2) {
464  Vec4 v;
465  v.xx = v1.yy * v2.zz - v1.zz * v2.yy;
466  v.yy = v1.zz * v2.xx - v1.xx * v2.zz;
467  v.zz = v1.xx * v2.yy - v1.yy * v2.xx; return v;
468 }
469 
470 
471 //--------------------------------------------------------------------------
472 
473 // Cross-product of three 4-vectors ( p_i = epsilon_{iabc} p_a p_b p_c)
474 
475 Vec4 cross4(const Vec4& a, const Vec4& b, const Vec4& c) {
476  Vec4 v(0.,0.,0.,0.);
477  v.tt = a.xx*b.yy*c.zz + a.yy*b.zz*c.xx + a.zz*b.xx*c.yy
478  - a.xx*b.zz*c.yy - a.zz*b.yy*c.xx - a.yy*b.xx*c.zz;
479  v.xx = -(- a.tt*b.yy*c.zz - a.yy*b.zz*c.tt - a.zz*b.tt*c.yy
480  + a.tt*b.zz*c.yy + a.zz*b.yy*c.tt + a.yy*b.tt*c.zz);
481  v.yy = -(- a.xx*b.tt*c.zz - a.tt*b.zz*c.xx - a.zz*b.xx*c.tt
482  + a.xx*b.zz*c.tt + a.zz*b.tt*c.xx + a.tt*b.xx*c.zz);
483  v.zz = -(- a.xx*b.yy*c.tt - a.yy*b.tt*c.xx - a.tt*b.xx*c.yy
484  + a.xx*b.tt*c.yy + a.tt*b.yy*c.xx + a.yy*b.xx*c.tt);
485  return v;
486 }
487 
488 //--------------------------------------------------------------------------
489 
490 // Opening angle between two three-vectors.
491 
492 double theta(const Vec4& v1, const Vec4& v2) {
493  double cthe = (v1.xx * v2.xx + v1.yy * v2.yy + v1.zz * v2.zz)
494  / sqrt( (v1.xx*v1.xx + v1.yy*v1.yy + v1.zz*v1.zz)
495  * (v2.xx*v2.xx + v2.yy*v2.yy + v2.zz*v2.zz) );
496  cthe = max(-1., min(1., cthe));
497  return acos(cthe);
498 }
499 
500 //--------------------------------------------------------------------------
501 
502 // Cosine of the opening angle between two three-vectors.
503 
504 double costheta(const Vec4& v1, const Vec4& v2) {
505  double cthe = (v1.xx * v2.xx + v1.yy * v2.yy + v1.zz * v2.zz)
506  / sqrt( (v1.xx*v1.xx + v1.yy*v1.yy + v1.zz*v1.zz)
507  * (v2.xx*v2.xx + v2.yy*v2.yy + v2.zz*v2.zz) );
508  cthe = max(-1., min(1., cthe));
509  return cthe;
510 }
511 
512 //--------------------------------------------------------------------------
513 
514 // Azimuthal angle between two three-vectors.
515 
516 double phi(const Vec4& v1, const Vec4& v2) {
517  double cphi = (v1.xx * v2.xx + v1.yy * v2.yy) / sqrt( max( Vec4::TINY,
518  (v1.xx*v1.xx + v1.yy*v1.yy) * (v2.xx*v2.xx + v2.yy*v2.yy) ));
519  cphi = max(-1., min(1., cphi));
520  return acos(cphi);
521 }
522 
523 //--------------------------------------------------------------------------
524 
525 // Cosine of the azimuthal angle between two three-vectors.
526 
527 double cosphi(const Vec4& v1, const Vec4& v2) {
528  double cphi = (v1.xx * v2.xx + v1.yy * v2.yy) / sqrt( max( Vec4::TINY,
529  (v1.xx*v1.xx + v1.yy*v1.yy) * (v2.xx*v2.xx + v2.yy*v2.yy) ));
530  cphi = max(-1., min(1., cphi));
531  return cphi;
532 }
533 
534 //--------------------------------------------------------------------------
535 
536 // Azimuthal angle between two three-vectors around a third.
537 
538 double phi(const Vec4& v1, const Vec4& v2, const Vec4& n) {
539  double nx = n.xx; double ny = n.yy; double nz = n.zz;
540  double norm = 1. / sqrt(nx*nx + ny*ny + nz*nz);
541  nx *= norm; ny *=norm; nz *=norm;
542  double v1s = v1.xx * v1.xx + v1.yy * v1.yy + v1.zz * v1.zz;
543  double v2s = v2.xx * v2.xx + v2.yy * v2.yy + v2.zz * v2.zz;
544  double v1v2 = v1.xx * v2.xx + v1.yy * v2.yy + v1.zz * v2.zz;
545  double v1n = v1.xx * nx + v1.yy * ny + v1.zz * nz;
546  double v2n = v2.xx * nx + v2.yy * ny + v2.zz * nz;
547  double cphi = (v1v2 - v1n * v2n) / sqrt( max( Vec4::TINY,
548  (v1s - v1n*v1n) * (v2s - v2n*v2n) ));
549  cphi = max(-1., min(1., cphi));
550  return acos(cphi);
551 }
552 
553 //--------------------------------------------------------------------------
554 
555 // Cosine of the azimuthal angle between two three-vectors around a third.
556 
557 double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n) {
558  double nx = n.xx; double ny = n.yy; double nz = n.zz;
559  double norm = 1. / sqrt(nx*nx + ny*ny + nz*nz);
560  nx *= norm; ny *=norm; nz *=norm;
561  double v1s = v1.xx * v1.xx + v1.yy * v1.yy + v1.zz * v1.zz;
562  double v2s = v2.xx * v2.xx + v2.yy * v2.yy + v2.zz * v2.zz;
563  double v1v2 = v1.xx * v2.xx + v1.yy * v2.yy + v1.zz * v2.zz;
564  double v1n = v1.xx * nx + v1.yy * ny + v1.zz * nz;
565  double v2n = v2.xx * nx + v2.yy * ny + v2.zz * nz;
566  double cphi = (v1v2 - v1n * v2n) / sqrt( max( Vec4::TINY,
567  (v1s - v1n*v1n) * (v2s - v2n*v2n) ));
568  cphi = max(-1., min(1., cphi));
569  return cphi;
570 }
571 
572 //--------------------------------------------------------------------------
573 
574 // Distance in cylindrical (y, phi) coordinates.
575 
576 double RRapPhi(const Vec4& v1, const Vec4& v2) {
577  double dRap = abs(v1.rap() - v2.rap());
578  double dPhi = abs(v1.phi() - v2.phi());
579  if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
580  return sqrt(dRap*dRap + dPhi*dPhi);
581 }
582 
583 //--------------------------------------------------------------------------
584 
585 // Distance in cylindrical (eta, phi) coordinates.
586 
587 double REtaPhi(const Vec4& v1, const Vec4& v2) {
588  double dEta = abs(v1.eta() - v2.eta());
589  double dPhi = abs(v1.phi() - v2.phi());
590  if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
591  return sqrt(dEta*dEta + dPhi*dPhi);
592 }
593 
594 //--------------------------------------------------------------------------
595 
596 // Shift four-momenta within pair from old to new masses.
597 // Note that p1Move and p2Move change values during operation.
598 
599 bool pShift( Vec4& p1Move, Vec4& p2Move, double m1New, double m2New) {
600 
601  // Standard kinematics variables.
602  double sH = (p1Move + p2Move).m2Calc();
603  double r1 = p1Move.m2Calc() / sH;
604  double r2 = p2Move.m2Calc() / sH;
605  double r3 = m1New * m1New / sH;
606  double r4 = m2New * m2New / sH;
607  double l12 = sqrtpos(pow2(1. - r1 - r2) - 4. * r1 * r2);
608  double l34 = sqrtpos(pow2(1. - r3 - r4) - 4. * r3 * r4);
609 
610  // Check that shift operation possible.
611  if (sH <= pow2(m1New + m2New) || l12 < Vec4::TINY || l34 < Vec4::TINY)
612  return false;
613 
614  // Calculate needed shift and apply it.
615  double c1 = 0.5 * ( (1. - r1 + r2) * l34 / l12 - (1. - r3 + r4) );
616  double c2 = 0.5 * ( (1. + r1 - r2) * l34 / l12 - (1. + r3 - r4) );
617  Vec4 pSh = c1 * p1Move - c2 * p2Move;
618  p1Move += pSh;
619  p2Move -= pSh;
620  return true;
621 }
622 
623 //--------------------------------------------------------------------------
624 
625 // Create two vectors that are perpendicular to both input vectors.
626 
627 pair<Vec4,Vec4> getTwoPerpendicular(const Vec4& v1, const Vec4& v2) {
628 
629  // One perpendicular vector from three-dimensional cross-product.
630  Vec4 nPerp( cross3(v1,v2) );
631  double TINY = std::numeric_limits<double>::epsilon();
632  if ( abs(nPerp.pAbs()) < TINY) {
633  Vec4 aux;
634  if (v1.px() != 0.) aux.p(v1.yy,v1.px(),v1.pz(),v1.e());
635  else if (v1.py() != 0.) aux.p(v1.px(),v1.pz(),v1.py(),v1.e());
636  else if (v1.pz() != 0.) aux.p(v1.pz(),v1.py(),v1.px(),v1.e());
637  nPerp.p( cross3(v1,aux) );
638  }
639  nPerp /= abs(nPerp.pAbs());
640 
641  // Second perpendicular vector from four-dimensional cross-product.
642  Vec4 lPerp( cross4(v1,v2,nPerp) );
643  lPerp /= sqrt(abs(lPerp.m2Calc()));
644  return make_pair(nPerp,lPerp);
645 }
646 
647 //==========================================================================
648 
649 // RotBstMatrix class.
650 // This class implements 4 * 4 matrices that encode an arbitrary combination
651 // of rotations and boosts, that can be applied to Vec4 four-vectors.
652 
653 //--------------------------------------------------------------------------
654 
655 // Constants: could be changed here if desired, but normally should not.
656 // These are of technical nature, as described for each.
657 
658 // Small number to avoid division by zero.
659 const double RotBstMatrix::TINY = 1e-20;
660 
661 //--------------------------------------------------------------------------
662 
663 // Rotate by polar angle theta and azimuthal angle phi.
664 
665 void RotBstMatrix::rot(double theta, double phi) {
666 
667  // Set up rotation matrix.
668  double cthe = cos(theta);
669  double sthe = sin(theta);
670  double cphi = cos(phi);
671  double sphi = sin(phi);
672  double Mrot[4][4] = {
673  {1., 0., 0., 0.},
674  {0., cthe * cphi, - sphi, sthe * cphi},
675  {0., cthe * sphi, cphi, sthe * sphi},
676  {0., -sthe, 0., cthe } };
677 
678  // Rotate current matrix accordingly.
679  double Mtmp[4][4];
680  for (int i = 0; i < 4; ++i)
681  for (int j = 0; j < 4; ++j)
682  Mtmp[i][j] = M[i][j];
683  for (int i = 0; i < 4; ++i)
684  for (int j = 0; j < 4; ++j)
685  M[i][j] = Mrot[i][0] * Mtmp[0][j] + Mrot[i][1] * Mtmp[1][j]
686  + Mrot[i][2] * Mtmp[2][j] + Mrot[i][3] * Mtmp[3][j];
687 
688 }
689 
690 //--------------------------------------------------------------------------
691 
692 // Rotate so that vector originally along z axis becomes parallel with p.
693 
694 void RotBstMatrix::rot(const Vec4& p) {
695 
696  double theta = p.theta();
697  double phi = p.phi();
698  rot(0., -phi);
699  rot(theta, phi);
700 
701 }
702 
703 //--------------------------------------------------------------------------
704 
705 // Boost with velocity vector (betaX, betaY, betaZ).
706 
707 void RotBstMatrix::bst(double betaX, double betaY, double betaZ) {
708 
709  // Set up boost matrix.
710  double gm = 1. / sqrt( max( TINY, 1. - betaX*betaX - betaY*betaY
711  - betaZ*betaZ ) );
712  double gf = gm*gm / (1. + gm);
713  double Mbst[4][4] = {
714  { gm, gm*betaX, gm*betaY, gm*betaZ },
715  { gm*betaX, 1. + gf*betaX*betaX, gf*betaX*betaY, gf*betaX*betaZ },
716  { gm*betaY, gf*betaY*betaX, 1. + gf*betaY*betaY, gf*betaY*betaZ },
717  { gm*betaZ, gf*betaZ*betaX, gf*betaZ*betaY, 1. + gf*betaZ*betaZ } };
718 
719  // Boost current matrix correspondingly.
720  double Mtmp[4][4];
721  for (int i = 0; i < 4; ++i)
722  for (int j = 0; j < 4; ++j)
723  Mtmp[i][j] = M[i][j];
724  for (int i = 0; i < 4; ++i)
725  for (int j = 0; j < 4; ++j)
726  M[i][j] = Mbst[i][0] * Mtmp[0][j] + Mbst[i][1] * Mtmp[1][j]
727  + Mbst[i][2] * Mtmp[2][j] + Mbst[i][3] * Mtmp[3][j];
728 
729 }
730 
731 //--------------------------------------------------------------------------
732 
733 // Boost so that vector originally at rest obtains same velocity as p.
734 
735 void RotBstMatrix::bst(const Vec4& p) {
736  double betaX = p.px() / p.e();
737  double betaY = p.py() / p.e();
738  double betaZ = p.pz() / p.e();
739  bst(betaX, betaY, betaZ);
740 }
741 
742 //--------------------------------------------------------------------------
743 
744 // Boost so vector originally with same velocity as p is brought to rest.
745 
746 void RotBstMatrix::bstback(const Vec4& p) {
747  double betaX = -p.px() / p.e();
748  double betaY = -p.py() / p.e();
749  double betaZ = -p.pz() / p.e();
750  bst(betaX, betaY, betaZ);
751 }
752 
753 //--------------------------------------------------------------------------
754 
755 // Boost that transforms p1 to p2, where p1^2 = p2^2 is assumed.
756 
757 void RotBstMatrix::bst(const Vec4& p1, const Vec4& p2) {
758  double eSum = p1.e() + p2.e();
759  double betaX = (p2.px() - p1.px()) / eSum;
760  double betaY = (p2.py() - p1.py()) / eSum;
761  double betaZ = (p2.pz() - p1.pz()) / eSum;
762  double fac = 2. / (1. + betaX*betaX + betaY*betaY + betaZ*betaZ);
763  betaX *= fac; betaY *= fac; betaZ *= fac;
764  bst(betaX, betaY, betaZ);
765 }
766 
767 //--------------------------------------------------------------------------
768 
769 // Boost and rotation that transforms from p1 and p2
770 // to their rest frame with p1 along +z axis.
771 
772 void RotBstMatrix::toCMframe(const Vec4& p1, const Vec4& p2) {
773  Vec4 pSum = p1 + p2;
774  Vec4 dir = p1;
775  dir.bstback(pSum);
776  double theta = dir.theta();
777  double phi = dir.phi();
778  bstback(pSum);
779  rot(0., -phi);
780  rot(-theta, phi);
781 }
782 
783 //--------------------------------------------------------------------------
784 
785 // Rotation and boost that transforms from rest frame of p1 and p2
786 // with p1 along +z axis to actual frame of p1 and p2. (Inverse of above.)
787 
788 void RotBstMatrix::fromCMframe(const Vec4& p1, const Vec4& p2) {
789  Vec4 pSum = p1 + p2;
790  Vec4 dir = p1;
791  dir.bstback(pSum);
792  double theta = dir.theta();
793  double phi = dir.phi();
794  rot(0., -phi);
795  rot(theta, phi);
796  bst(pSum);
797 }
798 
799 //--------------------------------------------------------------------------
800 
801 // Combine existing rotation/boost matrix with another one.
802 
803 void RotBstMatrix::rotbst(const RotBstMatrix& Mrb) {
804  double Mtmp[4][4];
805  for (int i = 0; i < 4; ++i)
806  for (int j = 0; j < 4; ++j)
807  Mtmp[i][j] = M[i][j];
808  for (int i = 0; i < 4; ++i)
809  for (int j = 0; j < 4; ++j)
810  M[i][j] = Mrb.M[i][0] * Mtmp[0][j] + Mrb.M[i][1] * Mtmp[1][j]
811  + Mrb.M[i][2] * Mtmp[2][j] + Mrb.M[i][3] * Mtmp[3][j];
812 }
813 
814 //--------------------------------------------------------------------------
815 
816 // Invert the rotation and boost.
817 
818 void RotBstMatrix::invert() {
819  double Mtmp[4][4];
820  for (int i = 0; i < 4; ++i)
821  for (int j = 0; j < 4; ++j)
822  Mtmp[i][j] = M[i][j];
823  for (int i = 0; i < 4; ++i)
824  for (int j = 0; j < 4; ++j)
825  M[i][j] = ( (i == 0 && j > 0) || (i > 0 && j == 0) )
826  ? - Mtmp[j][i] : Mtmp[j][i];
827 }
828 
829 //--------------------------------------------------------------------------
830 
831 // Reset to diagonal matrix.
832 
833 void RotBstMatrix::reset() {
834  for (int i = 0; i < 4; ++i)
835  for (int j = 0; j < 4; ++j)
836  M[i][j] = (i==j) ? 1. : 0.;
837 }
838 
839 //--------------------------------------------------------------------------
840 
841 // Crude estimate deviation from unit matrix.
842 
843 double RotBstMatrix::deviation() const {
844  double devSum = 0.;
845  for (int i = 0; i < 4; ++i)
846  for (int j = 0; j < 4; ++j)
847  devSum += (i==j) ? abs(M[i][j] - 1.) : abs(M[i][j]);
848  return devSum;
849 }
850 
851 //--------------------------------------------------------------------------
852 
853 // Print a rotation and boost matrix: operator overloading with friend.
854 
855 ostream& operator<<(ostream& os, const RotBstMatrix& M) {
856  os << fixed << setprecision(5) << " Rotation/boost matrix: \n";
857  for (int i = 0; i <4; ++i)
858  os << setw(10) << M.M[i][0] << setw(10) << M.M[i][1]
859  << setw(10) << M.M[i][2] << setw(10) << M.M[i][3] << "\n";
860  return os;
861 }
862 
863 //==========================================================================
864 
865 // Hist class.
866 // This class handles a single histogram at a time
867 // (or a vector of histograms).
868 
869 //--------------------------------------------------------------------------
870 
871 // Constants: could be changed here if desired, but normally should not.
872 // These are of technical nature, as described for each.
873 
874 // Maximum number of bins in a histogram.
875 const int Hist::NBINMAX = 10000;
876 
877 // Maximum number of columns that can be printed for a histogram.
878 const int Hist::NCOLMAX = 100;
879 
880 // Maximum number of lines a histogram can use at output.
881 const int Hist::NLINES = 30;
882 
883 // Tolerance in deviation of xMin and xMax between two histograms.
884 const double Hist::TOLERANCE = 0.001;
885 
886 // Small and large numbers to avoid division by zero and overflow.
887 const double Hist::TINY = 1e-20;
888 const double Hist::LARGE = 1e20;
889 
890 // When minbin/maxbin < SMALLFRAC the y scale goes down to zero.
891 const double Hist::SMALLFRAC = 0.1;
892 
893 // Constants for printout: fixed steps on y scale; filling characters.
894 const double DYAC[] = {0.04, 0.05, 0.06, 0.08, 0.10,
895  0.12, 0.15, 0.20, 0.25, 0.30};
896 const char NUMBER[] = {'0', '1', '2', '3', '4', '5',
897  '6', '7', '8', '9', 'X' };
898 
899 //--------------------------------------------------------------------------
900 
901 // Book a histogram.
902 
903 void Hist::book(string titleIn, int nBinIn, double xMinIn,
904  double xMaxIn, bool logXIn) {
905 
906  titleSave = titleIn;
907  nBin = nBinIn;
908  if (nBinIn < 1) nBin = 1;
909  if (nBinIn > NBINMAX) {
910  nBin = NBINMAX;
911  cout << " Warning: number of bins for histogram " << titleIn
912  << " reduced to " << nBin << endl;
913  }
914  linX = !logXIn;
915  xMin = xMinIn;
916  xMax = xMaxIn;
917  if (!linX && xMin < TINY) {
918  xMin = TINY;
919  cout << " Warning: lower x border of histogram " << titleIn
920  << " increased to " << xMin << endl;
921  }
922  if (xMax < xMin + TINY) {
923  xMax = 2. * xMin;
924  cout << " Warning: upper x border of histogram " << titleIn
925  << " increased to " << xMax << endl;
926  }
927  dx = (linX) ? (xMax - xMin) / nBin : log10(xMax / xMin) / nBin;
928  res.resize(nBin);
929  null();
930 
931 }
932 
933 //--------------------------------------------------------------------------
934 
935 // Reset bin contents.
936 
937 void Hist::null() {
938 
939  nFill = 0;
940  under = 0.;
941  inside = 0.;
942  over = 0.;
943  for (int ix = 0; ix < nBin; ++ix) res[ix] = 0.;
944 
945 }
946 
947 //--------------------------------------------------------------------------
948 
949 // Fill bin with weight.
950 
951 void Hist::fill(double x, double w) {
952 
953  ++nFill;
954  if (x < xMin) {under += w; return;}
955  if (x > xMax) {over += w; return;}
956  int iBin = (linX) ? int( floor( (x - xMin) / dx) )
957  : int( floor( log10(x / xMin) / dx) );
958  if (iBin < 0) under += w;
959  else if (iBin >= nBin) over += w;
960  else {inside += w; res[iBin] += w; }
961 
962 }
963 
964 //--------------------------------------------------------------------------
965 
966 // Print a histogram: also operator overloading with friend.
967 
968 ostream& operator<<(ostream& os, const Hist& h) {
969 
970  // Do not print empty histograms.
971  if (h.nFill <= 0) return os;
972 
973  // Write time and title.
974  time_t t = time(0);
975  char date[18];
976  strftime(date,18,"%Y-%m-%d %H:%M",localtime(&t));
977  os << "\n\n " << date << " " << h.titleSave << "\n\n";
978 
979  // Group bins, where required, to make printout have fewer columns.
980  // Avoid overflow.
981  int nGroup = 1 + (h.nBin - 1) / Hist::NCOLMAX;
982  int nCol = 1 + (h.nBin - 1) / nGroup;
983  vector<double> resCol(nCol);
984  for (int iCol = 0; iCol < nCol; ++iCol) {
985  resCol[iCol] = 0.;
986  for (int ix = nGroup * iCol; ix < min( h.nBin, nGroup * (iCol + 1)); ++ix)
987  resCol[iCol] += h.res[ix];
988  resCol[iCol] = max( -Hist::LARGE, min( Hist::LARGE, resCol[iCol] ) );
989  }
990 
991  // Find minimum and maximum bin content.
992  double yMin = resCol[0];
993  double yMax = resCol[0];
994  for (int iCol = 1; iCol < nCol; ++iCol) {
995  if (resCol[iCol] < yMin) yMin = resCol[iCol];
996  if (resCol[iCol] > yMax) yMax = resCol[iCol];
997  }
998 
999  // Determine scale and step size for y axis.
1000  if (yMax - yMin > Hist::NLINES * DYAC[0] * 1e-9) {
1001  if (yMin > 0. && yMin < Hist::SMALLFRAC * yMax) yMin = 0.;
1002  if (yMax < 0. && yMax > Hist::SMALLFRAC * yMin) yMax = 0.;
1003  int iPowY = int(floor( log10(yMax - yMin) ));
1004  if (yMax - yMin < Hist::NLINES * DYAC[0] * pow(10.,iPowY))
1005  iPowY = iPowY - 1;
1006  if (yMax - yMin > Hist::NLINES * DYAC[9] * pow(10.,iPowY))
1007  iPowY = iPowY + 1;
1008  double nLinePow = Hist::NLINES * pow(10.,iPowY);
1009  double delY = DYAC[0];
1010  for (int idel = 0; idel < 9; ++idel)
1011  if (yMax - yMin >= nLinePow * DYAC[idel]) delY = DYAC[idel+1];
1012  double dy = delY * pow(10.,iPowY);
1013 
1014  // Convert bin contents to integer form; fractional fill in top row.
1015  vector<int> row(nCol);
1016  vector<int> frac(nCol);
1017  for (int iCol = 0; iCol < nCol ; ++iCol) {
1018  double cta = abs(resCol[iCol]) / dy;
1019  row[iCol] = int(cta + 0.95);
1020  if(resCol[iCol] < 0.) row[iCol] = - row[iCol];
1021  frac[iCol] = int(10. * (cta + 1.05 - floor(cta + 0.95)));
1022  }
1023  int rowMin = int(abs(yMin)/dy + 0.95);
1024  if ( yMin < 0) rowMin = - rowMin;
1025  int rowMax = int(abs(yMax)/dy + 0.95);
1026  if ( yMax < 0) rowMax = - rowMax;
1027 
1028  // Print histogram row by row.
1029  os << fixed << setprecision(2);
1030  for (int iRow = rowMax; iRow >= rowMin; iRow--) if (iRow != 0) {
1031  os << " " << setw(10) << iRow*delY << "*10^"
1032  << setw(2) << iPowY << " ";
1033  for (int iCol = 0; iCol < nCol ; ++iCol) {
1034  if (iRow == row[iCol]) os << NUMBER[frac[iCol]];
1035  else if (iRow * (row[iCol] - iRow) > 0) os << NUMBER[10];
1036  else os << " ";
1037  } os << "\n";
1038  } os << "\n";
1039 
1040  // Print sign and value of bin contents
1041  double maxim = log10(max(yMax, -yMin));
1042  int iPowBin = int(floor(maxim + 0.0001));
1043  os << " Contents ";
1044  for (int iCol = 0; iCol < nCol ; ++iCol) {
1045  if (resCol[iCol] < - pow(10., iPowBin - 4)) os << "-";
1046  else os << " ";
1047  row[iCol] = int(abs(resCol[iCol]) * pow(10., 3 - iPowBin) + 0.5);
1048  } os << "\n";
1049  for (int iRow = 3; iRow >= 0; iRow--) {
1050  os << " *10^" << setw(2) << iPowBin + iRow - 3 << " ";
1051  int mask = int( pow(10., iRow) + 0.5);
1052  for (int iCol = 0; iCol < nCol ; ++iCol) {
1053  os << NUMBER[(row[iCol] / mask) % 10];
1054  } os << "\n";
1055  } os << "\n";
1056 
1057  // Print sign and value of lower bin edge.
1058  maxim = (h.linX) ? log10( max( -h.xMin, h.xMax - h.dx))
1059  : log10( max( -log10(h.xMin), log10(h.xMax) - h.dx ) );
1060  int iPowExp = int(floor(maxim + 0.0001));
1061  os << " Low edge ";
1062  for (int iCol = 0; iCol < nCol ; ++iCol) {
1063  double edgeNow = (h.linX) ? h.xMin + iCol * nGroup * h.dx
1064  : log10(h.xMin) + iCol * nGroup * h.dx;
1065  os << ( ( edgeNow < - pow(10., iPowExp - 3) ) ? "-" : " " );
1066  row[iCol] = int( abs(edgeNow) * pow(10., 2 - iPowExp) + 0.5 );
1067  } os << "\n";
1068  for (int iRow = 2; iRow >= 0; iRow--) {
1069  os << " *10^" << setw(2) << iPowExp + iRow - 2 << " ";
1070  int mask = int( pow(10., iRow) + 0.5);
1071  for (int iCol = 0; iCol < nCol ; ++iCol)
1072  os << NUMBER[(row[iCol] / mask) % 10];
1073  os << "\n";
1074  } os << "\n";
1075 
1076  // Print explanation if histogram cannot be shown.
1077  } else os << " Histogram not shown since lowest value" << scientific
1078  << setprecision(4) << setw(12) << yMin << " and highest value"
1079  << setw(12) << yMax << " are too close \n \n";
1080 
1081  // Calculate and print statistics.
1082  double cSum = 0.;
1083  double cxSum = 0.;
1084  double cxxSum = 0.;
1085  for (int ix = 0; ix < h.nBin ; ++ix) {
1086  double cta = abs(h.res[ix]);
1087  double x = (h.linX) ? h.xMin + (ix + 0.5) * h.dx
1088  : h.xMin * pow( 10., (ix + 0.5) * h.dx);
1089  cSum = cSum + cta;
1090  cxSum = cxSum + cta * x;
1091  cxxSum = cxxSum + cta * x * x;
1092  }
1093  double xmean = cxSum / max(cSum, Hist::TINY);
1094  double rms = sqrtpos( cxxSum / max(cSum, Hist::TINY) - xmean*xmean );
1095  os << scientific << setprecision(4)
1096  << " Entries =" << setw(12) << h.nFill
1097  << " Mean =" << setw(12) << xmean
1098  << " Underflow =" << setw(12) << h.under
1099  << " Low edge =" << setw(12) << h.xMin << "\n"
1100  << " All chan =" << setw(12) << h.inside
1101  << " Rms =" << setw(12) << rms
1102  << " Overflow =" << setw(12) << h.over
1103  << " High edge =" << setw(12) << h.xMax << endl;
1104  return os;
1105 }
1106 
1107 //--------------------------------------------------------------------------
1108 
1109 // Print histogram contents as a table (e.g. for Gnuplot).
1110 
1111 void Hist::table(ostream& os, bool printOverUnder, bool xMidBin) const {
1112 
1113  // Print histogram vector bin by bin, with mean x as first column.
1114  os << scientific << setprecision(4);
1115  double xBeg = (xMidBin) ? xMin + 0.5 * dx : xMin;
1116  if (!linX && xMidBin) xBeg = xMin * pow( 10., 0.5 * dx);
1117  if (printOverUnder)
1118  os << setw(12) << (linX ? xBeg - dx : xBeg * pow(10., -dx))
1119  << setw(12) << under << "\n";
1120  for (int ix = 0; ix < nBin; ++ix)
1121  os << setw(12) << (linX ? xBeg + ix * dx : xBeg * pow(10., ix * dx))
1122  << setw(12) << res[ix] << "\n";
1123  if (printOverUnder)
1124  os << setw(12) << (linX ? xBeg + nBin * dx : xBeg * pow(10., nBin * dx))
1125  << setw(12) << over << "\n";
1126 
1127 }
1128 
1129 //--------------------------------------------------------------------------
1130 
1131 // Print histogram contents as a table, in Rivet's *.dat style.
1132 
1133 void Hist::rivetTable(ostream& os, bool printError) const {
1134 
1135  // Print histogram vector bin by bin, with x range in first two columns
1136  // and +- error in last two (assuming that contents is number of events).
1137  os << scientific << setprecision(4);
1138  double xBeg = xMin;
1139  double xEnd = (linX) ? xMin + dx : xMin * pow(10., dx);
1140  for (int ix = 0; ix < nBin; ++ix) {
1141  double err = (printError) ? sqrtpos(res[ix]) : 0.0;
1142  os << setw(12) << (linX ? xBeg + ix * dx : xBeg * pow(10., ix * dx))
1143  << setw(12) << (linX ? xEnd + ix * dx : xEnd * pow(10., ix * dx))
1144  << setw(12) << res[ix] << setw(12) << err << setw(12) << err << "\n";
1145  }
1146 
1147 }
1148 
1149 //--------------------------------------------------------------------------
1150 
1151 // Print histogram contents as a table, as appropriate for Pyplot.
1152 
1153 void Hist::pyplotTable(ostream& os, bool isHist) const {
1154 
1155  // Set precision.
1156  os << scientific << setprecision(4);
1157 
1158  // For plotting as a histogram one needs bin edges as last column.
1159  double xBeg = (linX) ? xMin + 0.5 * dx : xMin * pow( 10., 0.5 * dx);
1160  double xNow, xEdge;
1161  for (int ix = 0; ix < nBin; ++ix) {
1162  xNow = (linX) ? xBeg + ix * dx : xBeg * pow(10., ix * dx);
1163  xEdge = (linX) ? xMin + ix * dx : xMin * pow(10., ix * dx);
1164  os << setw(12) << xNow << setw(12) << res[ix];
1165  if (isHist) os << setw(12) << xEdge << "\n";
1166  else os << "\n";
1167  }
1168 
1169  // And also an extra no-weights line to give final upper bin edge.
1170  if (isHist) {
1171  double xEnd = (linX) ? xMax - 0.5 * dx : xMax * pow( 10., -0.5 * dx);
1172  os << setw(12) << xEnd << setw(12) << 0. << setw(12) << xMax << "\n";
1173  }
1174 
1175 }
1176 
1177 //--------------------------------------------------------------------------
1178 
1179 // Print a table out of two histograms with same x axis (e.g. for Gnuplot).
1180 
1181 void table(const Hist& h1, const Hist& h2, ostream& os, bool printOverUnder,
1182  bool xMidBin) {
1183 
1184  // Require histogram x axes to agree.
1185  int nBin = h1.nBin;
1186  double dx = h1.dx;
1187  if (nBin != h2.nBin || abs(h1.xMin - h2.xMin) > Hist::TOLERANCE * dx
1188  || abs(h1.xMax - h2.xMax) > Hist::TOLERANCE * dx || h1.linX != h2.linX)
1189  return;
1190 
1191  // Print histogram vectors bin by bin, with mean x as first column.
1192  os << scientific << setprecision(4);
1193  double xBeg = (xMidBin) ? h1.xMin + 0.5 * dx : h1.xMin;
1194  if (!h1.linX && xMidBin) xBeg = h1.xMin * pow(10., 0.5 * dx);
1195  if (printOverUnder)
1196  os << setw(12) << (h1.linX ? xBeg - dx : xBeg * pow(10., -dx))
1197  << setw(12) << h1.under << setw(12) << h2.under << "\n";
1198  for (int ix = 0; ix < nBin; ++ix)
1199  os << setw(12) << (h1.linX ? xBeg + ix * dx : xBeg * pow(10., ix * dx))
1200  << setw(12) << h1.res[ix] << setw(12) << h2.res[ix] << "\n";
1201  if (printOverUnder)
1202  os << setw(12) << (h1.linX ? xBeg + nBin * dx : xBeg * pow(10., nBin * dx))
1203  << setw(12) << h1.over << setw(12) << h2.over << "\n";
1204 
1205 }
1206 
1207 void table(const Hist& h1, const Hist& h2, string fileName,
1208  bool printOverUnder, bool xMidBin) {
1209 
1210  ofstream streamName(fileName.c_str());
1211  table( h1, h2, streamName, printOverUnder, xMidBin);
1212 
1213 }
1214 
1215 //--------------------------------------------------------------------------
1216 
1217 // Get content of specific bin.
1218 // Special values are bin 0 for underflow and bin nBin+1 for overflow.
1219 // All other bins outside proper histogram range return 0.
1220 
1221 double Hist::getBinContent(int iBin) const {
1222 
1223  if (iBin > 0 && iBin <= nBin) return res[iBin - 1];
1224  else if (iBin == 0) return under;
1225  else if (iBin == nBin + 1) return over;
1226  else return 0.;
1227 
1228 }
1229 
1230 //--------------------------------------------------------------------------
1231 
1232 // Check whether another histogram has same size and limits.
1233 
1234 bool Hist::sameSize(const Hist& h) const {
1235 
1236  if (nBin == h.nBin && abs(xMin - h.xMin) < TOLERANCE * dx &&
1237  abs(xMax - h.xMax) < TOLERANCE * dx) return true;
1238  else return false;
1239 
1240 }
1241 
1242 //--------------------------------------------------------------------------
1243 
1244 // Take 10-logarithm or natural logarithm of contents bin by bin.
1245 
1246 void Hist::takeLog(bool tenLog) {
1247 
1248  // Find smallest positive bin content, and put min a bit below.
1249  double yMin = Hist::LARGE;
1250  for (int ix = 0; ix < nBin; ++ix)
1251  if (res[ix] > Hist::TINY && res[ix] < yMin ) yMin = res[ix];
1252  yMin *= 0.8;
1253 
1254  // Take 10-logarithm bin by bin, but ensure positivity.
1255  if (tenLog) {
1256  for (int ix = 0; ix < nBin; ++ix)
1257  res[ix] = log10( max( yMin, res[ix]) );
1258  under = log10( max( yMin, under) );
1259  inside = log10( max( yMin, inside) );
1260  over = log10( max( yMin, over) );
1261 
1262  // Take natural logarithm bin by bin, but ensure positivity.
1263  } else {
1264  for (int ix = 0; ix < nBin; ++ix)
1265  res[ix] = log( max( yMin, res[ix]) );
1266  under = log( max( yMin, under) );
1267  inside = log( max( yMin, inside) );
1268  over = log( max( yMin, over) );
1269  }
1270 
1271 }
1272 
1273 //--------------------------------------------------------------------------
1274 
1275 // Take square root of contents bin by bin; set 0 for negative content.
1276 
1277 void Hist::takeSqrt() {
1278 
1279  for (int ix = 0; ix < nBin; ++ix) res[ix] = sqrtpos(res[ix]);
1280  under = sqrtpos(under);
1281  inside = sqrtpos(inside);
1282  over = sqrtpos(over);
1283 
1284 }
1285 
1286 //--------------------------------------------------------------------------
1287 
1288 // Find smallest nonzero absolute value of bin contents.
1289 
1290 double Hist::smallestAbsValue() const {
1291 
1292  double smallest = 1e20; double yAbs;
1293  for (int ix = 0; ix < nBin; ++ix) { yAbs = abs(res[ix]);
1294  if (yAbs > 1e-20 && yAbs < smallest) smallest = yAbs; }
1295  return smallest;
1296 
1297 }
1298 
1299 //--------------------------------------------------------------------------
1300 
1301 // Add histogram to existing one.
1302 
1303 Hist& Hist::operator+=(const Hist& h) {
1304  if (!sameSize(h)) return *this;
1305  nFill += h.nFill;
1306  under += h.under;
1307  inside += h.inside;
1308  over += h.over;
1309  for (int ix = 0; ix < nBin; ++ix) res[ix] += h.res[ix];
1310  return *this;
1311 }
1312 
1313 //--------------------------------------------------------------------------
1314 
1315 // Subtract histogram from existing one.
1316 
1317 Hist& Hist::operator-=(const Hist& h) {
1318  if (!sameSize(h)) return *this;
1319  nFill += h.nFill;
1320  under -= h.under;
1321  inside -= h.inside;
1322  over -= h.over;
1323  for (int ix = 0; ix < nBin; ++ix) res[ix] -= h.res[ix];
1324  return *this;
1325 }
1326 
1327 //--------------------------------------------------------------------------
1328 
1329 // Multiply existing histogram by another one.
1330 
1331 Hist& Hist::operator*=(const Hist& h) {
1332  if (!sameSize(h)) return *this;
1333  nFill += h.nFill;
1334  under *= h.under;
1335  inside *= h.inside;
1336  over *= h.over;
1337  for (int ix = 0; ix < nBin; ++ix) res[ix] *= h.res[ix];
1338  return *this;
1339 }
1340 
1341 //--------------------------------------------------------------------------
1342 
1343 // Divide existing histogram by another one.
1344 
1345 Hist& Hist::operator/=(const Hist& h) {
1346  if (!sameSize(h)) return *this;
1347  nFill += h.nFill;
1348  under = (abs(h.under) < Hist::TINY) ? 0. : under/h.under;
1349  inside = (abs(h.inside) < Hist::TINY) ? 0. : inside/h.inside;
1350  over = (abs(h.over) < Hist::TINY) ? 0. : over/h.over;
1351  for (int ix = 0; ix < nBin; ++ix)
1352  res[ix] = (abs(h.res[ix]) < Hist::TINY) ? 0. : res[ix]/h.res[ix];
1353  return *this;
1354 }
1355 
1356 //--------------------------------------------------------------------------
1357 
1358 // Add constant offset to histogram.
1359 
1360 Hist& Hist::operator+=(double f) {
1361  under += f;
1362  inside += nBin * f;
1363  over += f;
1364  for (int ix = 0; ix < nBin; ++ix) res[ix] += f;
1365  return *this;
1366 }
1367 
1368 //--------------------------------------------------------------------------
1369 
1370 // Subtract constant offset from histogram.
1371 
1372 Hist& Hist::operator-=(double f) {
1373  under -= f;
1374  inside -= nBin * f;
1375  over -= f;
1376  for (int ix = 0; ix < nBin; ++ix) res[ix] -= f;
1377  return *this;
1378 }
1379 
1380 //--------------------------------------------------------------------------
1381 
1382 // Multiply histogram by constant.
1383 
1384 Hist& Hist::operator*=(double f) {
1385  under *= f;
1386  inside *= f;
1387  over *= f;
1388  for (int ix = 0; ix < nBin; ++ix) res[ix] *= f;
1389  return *this;
1390 }
1391 
1392 //--------------------------------------------------------------------------
1393 
1394 // Divide histogram by constant.
1395 
1396 Hist& Hist::operator/=(double f) {
1397  if (abs(f) > Hist::TINY) {
1398  under /= f;
1399  inside /= f;
1400  over /= f;
1401  for (int ix = 0; ix < nBin; ++ix) res[ix] /= f;
1402  // Set empty contents when division by zero.
1403  } else {
1404  under = 0.;
1405  inside = 0.;
1406  over = 0.;
1407  for (int ix = 0; ix < nBin; ++ix) res[ix] = 0.;
1408  }
1409  return *this;
1410 }
1411 
1412 Hist Hist::operator+(double f) const {
1413  Hist h = *this; return h += f;}
1414 
1415 Hist Hist::operator+(const Hist& h2) const {
1416  Hist h = *this; return h += h2;}
1417 
1418 Hist Hist::operator-(double f) const {
1419  Hist h = *this; return h -= f;}
1420 
1421 Hist Hist::operator-(const Hist& h2) const {
1422  Hist h = *this; return h -= h2;}
1423 
1424 Hist Hist::operator*(double f) const {
1425  Hist h = *this; return h *= f;}
1426 
1427 Hist Hist::operator*(const Hist& h2) const {
1428  Hist h = *this; return h *= h2;}
1429 
1430 Hist Hist::operator/(double f) const {
1431  Hist h = *this; return h /= f;}
1432 
1433 Hist Hist::operator/(const Hist& h2) const {
1434  Hist h = *this; return h /= h2;}
1435 
1436 //--------------------------------------------------------------------------
1437 
1438 // Implementation of operator overloading with friends.
1439 
1440 Hist operator+(double f, const Hist& h1) {
1441  Hist h = h1; return h += f;}
1442 
1443 Hist operator-(double f, const Hist& h1) {
1444  Hist h = h1;
1445  h.under = f - h1.under;
1446  h.inside = h1.nBin * f - h1.inside;
1447  h.over = f - h1.over;
1448  for (int ix = 0; ix < h1.nBin; ++ix) h.res[ix] = f - h1.res[ix];
1449  return h;}
1450 
1451 Hist operator*(double f, const Hist& h1) {
1452  Hist h = h1; return h *= f;}
1453 
1454 Hist operator/(double f, const Hist& h1) {
1455  Hist h = h1;
1456  h.under = (abs(h1.under) < Hist::TINY) ? 0. : f/h1.under;
1457  h.inside = (abs(h1.inside) < Hist::TINY) ? 0. : f/h1.inside;
1458  h.over = (abs(h1.over) < Hist::TINY) ? 0. : f/h1.over;
1459  for (int ix = 0; ix < h1.nBin; ++ix)
1460  h.res[ix] = (abs(h1.res[ix]) < Hist::TINY) ? 0. : f/h1.res[ix];
1461  return h;
1462 }
1463 
1464 //==========================================================================
1465 
1466 // HistPlot class.
1467 // Writes a Python program that can generate PDF plots from Hist histograms.
1468 
1469 //--------------------------------------------------------------------------
1470 
1471 // Generate the Python code for plotting a frame.
1472 
1473 void HistPlot::plot( bool logY) {
1474 
1475  // Start new file or add to existing one.
1476  if (frameName != "") {
1477  if (nPDF > 0) toPython << "pp.close()" << endl;
1478  ++nPDF;
1479  fileName = frameName;
1480  toPython << "pp = PdfPages('" << fileName << ".pdf')" << endl;
1481  nFrame = 1;
1482  nTable = 0;
1483  } else {
1484  ++nFrame;
1485  }
1486  toPython << "tmp" << nFrame << " = plt.figure(" << nFrame << ")" << endl;
1487 
1488  // Loop through the vector of histograms.
1489  double yAbsMin = 1e20;
1490  for (int iHist = 0; iHist < int(histos.size()); ++iHist) {
1491 
1492  // Histogram information for plotting.
1493  string legendNow = (legends[iHist] != "void") ? legends[iHist]
1494  : histos[iHist]->getTitle();
1495  stringstream nBin;
1496  nBin << histos[iHist]->getBinNumber();
1497  double yAbsNow = histos[iHist]->smallestAbsValue();
1498  if (yAbsNow < yAbsMin) yAbsMin = yAbsNow;
1499 
1500  // Split plotting style and potential colour information.
1501  string styleNow = (styles[iHist] == "") ? "h" : styles[iHist];
1502  string style1 = styleNow;
1503  string style2 = "";
1504  if (styleNow.find(",") != string::npos) {
1505  int iComma = styleNow.find(",");
1506  style1 = (iComma > 0) ? styleNow.substr( 0, iComma) : "h";
1507  if (iComma + 1 < int(styleNow.length()))
1508  style2 = styleNow.substr( iComma + 1);
1509  }
1510 
1511  // Write histogram itself to a data file as two columns of (x,y) values.
1512  stringstream encode;
1513  encode << fileName << "-" << nTable + iHist << ".dat";
1514  histos[iHist]->pyplotTable( encode.str(), (style1 == "h") );
1515 
1516  // Write code to plot histogram.
1517  toPython << "plot = open('" << encode.str() << "')" << endl
1518  << "plot = [line.split() for line in plot]" << endl
1519  << "valx = [float(x[0]) for x in plot]" << endl
1520  << "valy = [float(x[1]) for x in plot]" << endl;
1521  if (style1 == "h") toPython << "vale = [float(x[2]) for x in plot]"
1522  << endl << "plt.hist( valx, vale, weights = valy,"
1523  << " histtype='step',";
1524  else toPython << "plt.plot( valx, valy, '" << style1 << "',";
1525  if (style2 != "") toPython << " color='" << style2 << "',";
1526  toPython << " label=r'" << legendNow << "')" << endl;
1527  }
1528 
1529  // Write title, axes and create plot.
1530  if (!histos[0]->getLinX()) toPython << "plt.xscale('log')" << endl;
1531  if (logY) toPython << "plt.yscale('symlog', linthreshy=" << scientific
1532  << setprecision(2) << yAbsMin << ")" << endl;
1533  else toPython << "plt.ticklabel_format(axis='y', style='sci', "
1534  << "scilimits=(-2,3))" << endl;
1535  toPython << "plt.legend(frameon=False,loc='best')" << endl
1536  << "plt.title(r'" << title << "')" << endl
1537  << "plt.xlabel(r'" << xLabel << "')" << endl
1538  << "plt.ylabel(r\'" << yLabel << "')" << endl
1539  << "pp.savefig(tmp" << nFrame << ",bbox_inches='tight')"
1540  << endl << "plt.clf()" << endl;
1541 
1542  // Update table counter. Done.
1543  nTable += histos.size();
1544 }
1545 
1546 //==========================================================================
1547 
1548 } // end namespace Pythia8
Generated by   doxygen 1.8.5