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