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