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