using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / include / Basics.h
CommitLineData
5ad4eb21 1// Basics.h is a part of the PYTHIA event generator.
2// Copyright (C) 2008 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// Header file for basic, often-used helper classes.
7// RndmEngine: base class for external random number generators.
8// Rndm: random number generator (static member functions).
9// Vec4: simple four-vectors.
10// RotBstMatrix: matrices encoding rotations and boosts of Vec4 objects.
11// Hist: simple one-dimensional histograms.
12
13#ifndef Pythia8_Basics_H
14#define Pythia8_Basics_H
15
16#include "PythiaStdlib.h"
17
18namespace Pythia8 {
19
20//**************************************************************************
21
22// RndmEngine is the base class for external random number generators.
23// There is only one pure virtual method, that should do the generation.
24
25class RndmEngine {
26
27public:
28
29 // A pure virtual method, wherein the derived class method
30 // generates a random number uniformly distributed between 1 and 1.
31 virtual double flat() = 0;
32
33protected:
34
35 // Destructor.
36 virtual ~RndmEngine() {}
37
38};
39
40//**************************************************************************
41
42// Rndm class.
43// This class handles random number generation according to the
44// Marsaglia-Zaman-Tsang algorithm.
45
46class Rndm {
47
48public:
49
50 // Constructors.
51 Rndm() {}
52 Rndm(int seedIn) { init(seedIn);}
53
54 // Possibility to pass in pointer for external random number generation.
55 static bool rndmEnginePtr( RndmEngine* rndmPtrIn);
56
57 // Initialize, normally at construction or in first call.
58 static void init(int seedIn = 0) ;
59
60 // Generate next random number uniformly between 0 and 1.
61 static double flat() ;
62
63 // Generate random numbers according to exp(-x).
64 static double exp() { return -log(flat()) ;}
65
66 // Generate random numbers according to x * exp(-x).
67 static double xexp() { return -log(flat() * flat()) ;}
68
69 // Generate random numbers according to exp(-x^2/2).
70 static double gauss() ;
71
72 // Pick one option among vector of (positive) probabilities.
73 static int pick(const vector<double>& prob) ;
74
75private:
76
77 // State of the random number generator.
78 static bool initRndm, saveGauss;
79 static int i97, j97, defaultSeed;
80 static double u[97], c, cd, cm, save;
81
82 // Pointer for external random number generation.
83 static bool useExternalRndm;
84 static RndmEngine* rndmPtr;
85
86};
87
88//**************************************************************************
89
90// Forward reference to RotBstMatrix class.
91class RotBstMatrix;
92
93//**************************************************************************
94
95// Vec4 class.
96// This class implements four-vectors, in energy-momentum space.
97// (But can equally well be used to hold space-time four-vectors.)
98
99class Vec4 {
100
101public:
102
103 // Constructors.
104 Vec4(double xIn = 0., double yIn = 0., double zIn = 0., double tIn = 0.)
105 : xx(xIn), yy(yIn), zz(zIn), tt(tIn) { }
106 Vec4(const Vec4& v) : xx(v.xx), yy(v.yy), zz(v.zz), tt(v.tt) { }
107 Vec4& operator=(const Vec4& v) { if (this != &v) { xx = v.xx; yy = v.yy;
108 zz = v.zz; tt = v.tt; } return *this; }
109 Vec4& operator=(double value) { xx = value; yy = value; zz = value;
110 tt = value; return *this; }
111
112 // Member functions for input.
113 void reset() {xx = 0.; yy = 0.; zz = 0.; tt = 0.;}
114 void p(double xIn, double yIn, double zIn, double tIn)
115 {xx = xIn; yy = yIn; zz = zIn; tt = tIn;}
116 void p(Vec4 pIn) {xx = pIn.xx; yy = pIn.yy; zz = pIn.zz; tt = pIn.tt;}
117 void px(double xIn) {xx = xIn;}
118 void py(double yIn) {yy = yIn;}
119 void pz(double zIn) {zz = zIn;}
120 void e(double tIn) {tt = tIn;}
121
122 // Member functions for output.
123 double px() const {return xx;}
124 double py() const {return yy;}
125 double pz() const {return zz;}
126 double e() const {return tt;}
127 double mCalc() const {double temp = tt*tt - xx*xx - yy*yy - zz*zz;
128 return (temp >= 0) ? sqrt(temp) : -sqrt(-temp);}
129 double m2Calc() const {return tt*tt - xx*xx - yy*yy - zz*zz;}
130 double pT() const {return sqrt(xx*xx + yy*yy);}
131 double pT2() const {return xx*xx + yy*yy;}
132 double pAbs() const {return sqrt(xx*xx + yy*yy + zz*zz);}
133 double pAbs2() const {return xx*xx + yy*yy + zz*zz;}
134 double theta() const {return atan2(sqrt(xx*xx + yy*yy), zz);}
135 double phi() const {return atan2(yy,xx);}
136 double thetaXZ() const {return atan2(xx,zz);}
137 double pPlus() const {return tt + zz;}
138 double pMinus() const {return tt - zz;}
139
140 // Member functions that perform operations.
141 void rescale3(double fac) {xx *= fac; yy *= fac; zz *= fac;}
142 void rescale4(double fac) {xx *= fac; yy *= fac; zz *= fac; tt *= fac;}
143 void flip3() {xx = -xx; yy = -yy; zz = -zz;}
144 void flip4() {xx = -xx; yy = -yy; zz = -zz; tt = -tt;}
145 void rot(double thetaIn, double phiIn);
146 void rotaxis(double phiIn, double nx, double ny, double nz);
147 void rotaxis(double phiIn, const Vec4& n);
148 void bst(double betaX, double betaY, double betaZ);
149 void bst(double betaX, double betaY, double betaZ, double gamma);
150 void bst(const Vec4& pIn);
151 void bst(const Vec4& pIn, double mIn);
152 void bstback(const Vec4& pIn);
153 void bstback(const Vec4& pIn, double mIn);
154 void rotbst(const RotBstMatrix& M);
155
156 // Operator overloading with member functions
157 Vec4& operator-() {xx = -xx; yy = -yy; zz = -zz; tt = -tt; return *this;}
158 Vec4& operator+=(const Vec4& v) {xx += v.xx; yy += v.yy; zz += v.zz;
159 tt += v.tt; return *this;}
160 Vec4& operator-=(const Vec4& v) {xx -= v.xx; yy -= v.yy; zz -= v.zz;
161 tt -= v.tt; return *this;}
162 Vec4& operator*=(double f) {xx *= f; yy *= f; zz *= f;
163 tt *= f; return *this;}
164 Vec4& operator/=(double f) {xx /= f; yy /= f; zz /= f;
165 tt /= f; return *this;}
166
167 // Operator overloading with friends
168 friend Vec4 operator+(const Vec4& v1, const Vec4& v2);
169 friend Vec4 operator-(const Vec4& v1, const Vec4& v2);
170 friend Vec4 operator*(double f, const Vec4& v1);
171 friend Vec4 operator*(const Vec4& v1, double f);
172 friend Vec4 operator/(const Vec4& v1, double f);
173 friend double operator*(const Vec4& v1, const Vec4& v2);
174
175 // Invariant mass of a pair and its square.
176 friend double m(const Vec4& v1, const Vec4& v2);
177 friend double m2(const Vec4& v1, const Vec4& v2);
178
179 // Scalar and cross product of 3-vector parts.
180 friend double dot3(const Vec4& v1, const Vec4& v2);
181 friend Vec4 cross3(const Vec4& v1, const Vec4& v2);
182
183 // theta is polar angle between v1 and v2.
184 friend double theta(const Vec4& v1, const Vec4& v2);
185 friend double costheta(const Vec4& v1, const Vec4& v2);
186
187 // phi is azimuthal angle between v1 and v2 around z axis.
188 friend double phi(const Vec4& v1, const Vec4& v2);
189 friend double cosphi(const Vec4& v1, const Vec4& v2);
190
191 // phi is azimuthal angle between v1 and v2 around n axis.
192 friend double phi(const Vec4& v1, const Vec4& v2, const Vec4& n);
193 friend double cosphi(const Vec4& v1, const Vec4& v2, const Vec4& n);
194
195 // Print a four-vector
196 friend ostream& operator<<(ostream&, const Vec4& v) ;
197
198private:
199
200 // Constants: could only be changed in the code itself.
201 static const double TINY;
202
203 // The four-vector data members.
204 double xx, yy, zz, tt;
205
206};
207
208// Implementation of operator overloading with friends.
209
210inline Vec4 operator+(const Vec4& v1, const Vec4& v2)
211 {Vec4 v = v1 ; return v += v2;}
212
213inline Vec4 operator-(const Vec4& v1, const Vec4& v2)
214 {Vec4 v = v1 ; return v -= v2;}
215
216inline Vec4 operator*(double f, const Vec4& v1)
217 {Vec4 v = v1; return v *= f;}
218
219inline Vec4 operator*(const Vec4& v1, double f)
220 {Vec4 v = v1; return v *= f;}
221
222inline Vec4 operator/(const Vec4& v1, double f)
223 {Vec4 v = v1; return v /= f;}
224
225inline double operator*(const Vec4& v1, const Vec4& v2)
226 {return v1.tt*v2.tt - v1.xx*v2.xx - v1.yy*v2.yy - v1.zz*v2.zz;}
227
228//**************************************************************************
229
230// RotBstMatrix class.
231// This class implements 4 * 4 matrices that encode an arbitrary combination
232// of rotations and boosts, that can be applied to Vec4 four-vectors.
233
234class RotBstMatrix {
235
236public:
237
238 // Constructors.
239 RotBstMatrix() {for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j)
240 { M[i][j] = (i==j) ? 1. : 0.; } } }
241 RotBstMatrix(const RotBstMatrix& Min) {
242 for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
243 M[i][j] = Min.M[i][j]; } } }
244 RotBstMatrix& operator=(const RotBstMatrix& Min) {if (this != &Min) {
245 for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) {
246 M[i][j] = Min.M[i][j]; } } } return *this; }
247
248 // Member functions.
249 void rot(double = 0., double = 0.);
250 void rot(const Vec4& p);
251 void bst(double = 0., double = 0., double = 0.);
252 void bst(const Vec4&);
253 void bstback(const Vec4&);
254 void bst(const Vec4&, const Vec4&);
255 void toCMframe(const Vec4&, const Vec4&);
256 void fromCMframe(const Vec4&, const Vec4&);
257 void rotbst(const RotBstMatrix&);
258 void invert();
259 void reset();
260
261 // Crude estimate deviation from unit matrix.
262 double deviation() const;
263
264 // Print a transformation matrix.
265 friend ostream& operator<<(ostream&, const RotBstMatrix&) ;
266
267 // Private members to be accessible from Vec4.
268 friend class Vec4;
269
270private:
271
272 // Constants: could only be changed in the code itself.
273 static const double TINY;
274
275 // The rotation-and-boost matrix data members.
276 double M[4][4];
277
278};
279
280//**************************************************************************
281
282// Hist class.
283// This class handles a single histogram at a time.
284
285class Hist{
286
287public:
288
289 // Constructors, including copy constructors.
290 Hist() {;}
291 Hist(string titleIn, int nBinIn = 100, double xMinIn = 0.,
292 double xMaxIn = 1.) {
293 book(titleIn, nBinIn, xMinIn, xMaxIn);}
294 Hist(const Hist& h)
295 : title(h.title), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
296 xMax(h.xMax), dx(h.dx), under(h.under), inside(h.inside),
297 over(h.over), res(h.res) { }
298 Hist(string titleIn, const Hist& h)
299 : title(titleIn), nBin(h.nBin), nFill(h.nFill), xMin(h.xMin),
300 xMax(h.xMax), dx(h.dx), under(h.under), inside(h.inside),
301 over(h.over), res(h.res) { }
302 Hist& operator=(const Hist& h) { if(this != &h) {
303 nBin = h.nBin; nFill = h.nFill; xMin = h.xMin; xMax = h.xMax;
304 dx = h.dx; under = h.under; inside = h.inside; over = h.over;
305 res = h.res; } return *this; }
306
307 // Book a histogram.
308 void book(string titleIn = " ", int nBinIn = 100, double xMinIn = 0.,
309 double xMaxIn = 1.) ;
310
311 // Set title of a histogram.
312 void name(string titleIn = " ") {title = titleIn; }
313
314 // Reset bin contents.
315 void null() ;
316
317 // Fill bin with weight.
318 void fill(double x, double w = 1.) ;
319
320 // Return content of specific bin: -1 gives underflow and nBin overflow.
321 double getBinContent(int iBin) ;
322
323 // Return number of entries
324 double getEntries() {return nFill; }
325
326 // Print histogram contents as a table (e.g. for Gnuplot).
327 void table(ostream& os = cout) const ;
328
329 // Check whether another histogram has same size and limits.
330 bool sameSize(const Hist& h) const ;
331
332 // Take logarithm (base 10 or e) of bin contents.
333 void takeLog(bool tenLog = true) ;
334
335 // Operator overloading with member functions
336 Hist& operator+=(const Hist& h) ;
337 Hist& operator-=(const Hist& h) ;
338 Hist& operator*=(const Hist& h) ;
339 Hist& operator/=(const Hist& h) ;
340 Hist& operator+=(double f) ;
341 Hist& operator-=(double f) ;
342 Hist& operator*=(double f) ;
343 Hist& operator/=(double f) ;
344
345 // Operator overloading with friends
346 friend Hist operator+(double f, const Hist& h1);
347 friend Hist operator+(const Hist& h1, double f);
348 friend Hist operator+(const Hist& h1, const Hist& h2);
349 friend Hist operator-(double f, const Hist& h1);
350 friend Hist operator-(const Hist& h1, double f);
351 friend Hist operator-(const Hist& h1, const Hist& h2);
352 friend Hist operator*(double f, const Hist& h1);
353 friend Hist operator*(const Hist& h1, double f);
354 friend Hist operator*(const Hist& h1, const Hist& h2);
355 friend Hist operator/(double f, const Hist& h1);
356 friend Hist operator/(const Hist& h1, double f);
357 friend Hist operator/(const Hist& h1, const Hist& h2);
358
359 // Print a histogram with overloaded << operator.
360 friend ostream& operator<<(ostream& os, const Hist& h) ;
361
362private:
363
364 // Constants: could only be changed in the code itself.
365 static const int NBINMAX, NLINES;
366 static const double TOLERANCE, TINY, SMALLFRAC, DYAC[];
367 static const char NUMBER[];
368
369 // Properties and contents of a histogram.
370 string title;
371 int nBin, nFill;
372 double xMin, xMax, dx, under, inside, over;
373 vector<double> res;
374
375};
376
377//**************************************************************************
378
379} // end namespace Pythia8
380
381#endif // end Pythia8_Basics_H