source: HiSusy/trunk/Pythia8/pythia8170/src/Basics.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

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