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