- Update to pythia8140
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8140 / src / BeamShape.cxx
1 // BeamShape.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2010 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 BeamShape class. 
7
8 #include "BeamShape.h"
9
10 namespace Pythia8 {
11  
12 //==========================================================================
13
14 // The BeamShape class.
15
16 //--------------------------------------------------------------------------
17
18 // Initialize beam parameters.
19
20   void BeamShape::init( Settings& settings, Rndm* rndmPtrIn) { 
21
22   // Save pointer.
23   rndmPtr             = rndmPtrIn;
24
25   // Main flags.
26   allowMomentumSpread = settings.flag("Beams:allowMomentumSpread");
27   allowVertexSpread   = settings.flag("Beams:allowVertexSpread");
28
29   // Parameters for beam A momentum spread.
30   sigmaPxA            = settings.parm("Beams:sigmaPxA");
31   sigmaPyA            = settings.parm("Beams:sigmaPyA");
32   sigmaPzA            = settings.parm("Beams:sigmaPzA");
33   maxDevA             = settings.parm("Beams:maxDevA");
34
35   // Parameters for beam B momentum spread.
36   sigmaPxB            = settings.parm("Beams:sigmaPxB");
37   sigmaPyB            = settings.parm("Beams:sigmaPyB");
38   sigmaPzB            = settings.parm("Beams:sigmaPzB");
39   maxDevB             = settings.parm("Beams:maxDevB");
40  
41   // Parameters for beam vertex spread.
42   sigmaVertexX        = settings.parm("Beams:sigmaVertexX");
43   sigmaVertexY        = settings.parm("Beams:sigmaVertexY");
44   sigmaVertexZ        = settings.parm("Beams:sigmaVertexZ");
45   maxDevVertex        = settings.parm("Beams:maxDevVertex");
46   sigmaTime           = settings.parm("Beams:sigmaTime");
47   maxDevTime          = settings.parm("Beams:maxDevTime");
48  
49   // Parameters for beam vertex offset.
50   offsetX             = settings.parm("Beams:offsetVertexX"); 
51   offsetY             = settings.parm("Beams:offsetVertexY");
52   offsetZ             = settings.parm("Beams:offsetVertexZ");
53   offsetT             = settings.parm("Beams:offsetTime");
54
55 }
56
57 //--------------------------------------------------------------------------
58
59 // Set the two beam momentum deviations and the beam vertex.
60
61 void BeamShape::pick() { 
62
63   // Reset all values.
64   deltaPxA = deltaPyA = deltaPzA = deltaPxB = deltaPyB = deltaPzB
65     = vertexX = vertexY = vertexZ = vertexT = 0.;
66
67   // Set beam A momentum deviation by a three-dimensional Gaussian.
68   if (allowMomentumSpread) {
69     double totalDev, gauss;
70     do {
71       totalDev = 0.;
72       if (sigmaPxA > 0.) {
73         gauss     = rndmPtr->gauss();
74         deltaPxA  = sigmaPxA * gauss;
75         totalDev += gauss * gauss; 
76       }
77       if (sigmaPyA > 0.) {
78         gauss     = rndmPtr->gauss();
79         deltaPyA  = sigmaPyA * gauss;
80         totalDev += gauss * gauss; 
81       }
82       if (sigmaPzA > 0.) {
83         gauss     = rndmPtr->gauss();
84         deltaPzA  = sigmaPzA * gauss;
85         totalDev += gauss * gauss; 
86       }
87     } while (totalDev > maxDevA * maxDevA); 
88
89     // Set beam B momentum deviation by a three-dimensional Gaussian.
90     do {
91       totalDev = 0.;
92       if (sigmaPxB > 0.) {
93         gauss     = rndmPtr->gauss();
94         deltaPxB  = sigmaPxB * gauss;
95         totalDev += gauss * gauss; 
96       }
97       if (sigmaPyB > 0.) {
98         gauss     = rndmPtr->gauss();
99         deltaPyB  = sigmaPyB * gauss;
100         totalDev += gauss * gauss; 
101       }
102       if (sigmaPzB > 0.) {
103         gauss     = rndmPtr->gauss();
104         deltaPzB  = sigmaPzB * gauss;
105         totalDev += gauss * gauss; 
106       }
107     } while (totalDev > maxDevB * maxDevB); 
108   }  
109
110   // Set beam vertex location by a three-dimensional Gaussian.
111   if (allowVertexSpread) {
112     double totalDev, gauss;
113     do {
114       totalDev = 0.;
115       if (sigmaVertexX > 0.) {
116         gauss     = rndmPtr->gauss();
117         vertexX   = sigmaVertexX * gauss;
118         totalDev += gauss * gauss; 
119       }
120       if (sigmaVertexY > 0.) {
121         gauss     = rndmPtr->gauss();
122         vertexY   = sigmaVertexY * gauss;
123         totalDev += gauss * gauss; 
124       }
125       if (sigmaVertexZ > 0.) {
126         gauss     = rndmPtr->gauss();
127         vertexZ   = sigmaVertexZ * gauss;
128         totalDev += gauss * gauss; 
129       }
130     } while (totalDev > maxDevVertex * maxDevVertex);
131
132     // Set beam collision time by a Gaussian.
133     if (sigmaTime > 0.) {
134       do gauss    = rndmPtr->gauss();
135       while (abs(gauss) > maxDevTime);
136       vertexT     = sigmaTime * gauss;
137     }
138
139     // Add offset to beam vertex.
140     vertexX      += offsetX;
141     vertexY      += offsetY;
142     vertexZ      += offsetZ;
143     vertexT      += offsetT;
144   }  
145
146 }
147  
148 //==========================================================================
149
150 } // end namespace Pythia8
151