Line data Source code
1 : // BeamShape.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 BeamShape class.
7 :
8 : #include "Pythia8/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 0 : rndmPtr = rndmPtrIn;
24 :
25 : // Main flags.
26 0 : allowMomentumSpread = settings.flag("Beams:allowMomentumSpread");
27 0 : allowVertexSpread = settings.flag("Beams:allowVertexSpread");
28 :
29 : // Parameters for beam A momentum spread.
30 0 : sigmaPxA = settings.parm("Beams:sigmaPxA");
31 0 : sigmaPyA = settings.parm("Beams:sigmaPyA");
32 0 : sigmaPzA = settings.parm("Beams:sigmaPzA");
33 0 : maxDevA = settings.parm("Beams:maxDevA");
34 :
35 : // Parameters for beam B momentum spread.
36 0 : sigmaPxB = settings.parm("Beams:sigmaPxB");
37 0 : sigmaPyB = settings.parm("Beams:sigmaPyB");
38 0 : sigmaPzB = settings.parm("Beams:sigmaPzB");
39 0 : maxDevB = settings.parm("Beams:maxDevB");
40 :
41 : // Parameters for beam vertex spread.
42 0 : sigmaVertexX = settings.parm("Beams:sigmaVertexX");
43 0 : sigmaVertexY = settings.parm("Beams:sigmaVertexY");
44 0 : sigmaVertexZ = settings.parm("Beams:sigmaVertexZ");
45 0 : maxDevVertex = settings.parm("Beams:maxDevVertex");
46 0 : sigmaTime = settings.parm("Beams:sigmaTime");
47 0 : maxDevTime = settings.parm("Beams:maxDevTime");
48 :
49 : // Parameters for beam vertex offset.
50 0 : offsetX = settings.parm("Beams:offsetVertexX");
51 0 : offsetY = settings.parm("Beams:offsetVertexY");
52 0 : offsetZ = settings.parm("Beams:offsetVertexZ");
53 0 : offsetT = settings.parm("Beams:offsetTime");
54 :
55 0 : }
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 0 : deltaPxA = deltaPyA = deltaPzA = deltaPxB = deltaPyB = deltaPzB
65 0 : = vertexX = vertexY = vertexZ = vertexT = 0.;
66 :
67 : // Set beam A momentum deviation by a three-dimensional Gaussian.
68 0 : if (allowMomentumSpread) {
69 : double totalDev, gauss;
70 0 : do {
71 : totalDev = 0.;
72 0 : if (sigmaPxA > 0.) {
73 0 : gauss = rndmPtr->gauss();
74 0 : deltaPxA = sigmaPxA * gauss;
75 0 : totalDev += gauss * gauss;
76 0 : }
77 0 : if (sigmaPyA > 0.) {
78 0 : gauss = rndmPtr->gauss();
79 0 : deltaPyA = sigmaPyA * gauss;
80 0 : totalDev += gauss * gauss;
81 0 : }
82 0 : if (sigmaPzA > 0.) {
83 0 : gauss = rndmPtr->gauss();
84 0 : deltaPzA = sigmaPzA * gauss;
85 0 : totalDev += gauss * gauss;
86 0 : }
87 0 : } while (totalDev > maxDevA * maxDevA);
88 :
89 : // Set beam B momentum deviation by a three-dimensional Gaussian.
90 : do {
91 : totalDev = 0.;
92 0 : if (sigmaPxB > 0.) {
93 0 : gauss = rndmPtr->gauss();
94 0 : deltaPxB = sigmaPxB * gauss;
95 0 : totalDev += gauss * gauss;
96 0 : }
97 0 : if (sigmaPyB > 0.) {
98 0 : gauss = rndmPtr->gauss();
99 0 : deltaPyB = sigmaPyB * gauss;
100 0 : totalDev += gauss * gauss;
101 0 : }
102 0 : if (sigmaPzB > 0.) {
103 0 : gauss = rndmPtr->gauss();
104 0 : deltaPzB = sigmaPzB * gauss;
105 0 : totalDev += gauss * gauss;
106 0 : }
107 0 : } while (totalDev > maxDevB * maxDevB);
108 0 : }
109 :
110 : // Set beam vertex location by a three-dimensional Gaussian.
111 0 : if (allowVertexSpread) {
112 : double totalDev, gauss;
113 0 : do {
114 : totalDev = 0.;
115 0 : if (sigmaVertexX > 0.) {
116 0 : gauss = rndmPtr->gauss();
117 0 : vertexX = sigmaVertexX * gauss;
118 0 : totalDev += gauss * gauss;
119 0 : }
120 0 : if (sigmaVertexY > 0.) {
121 0 : gauss = rndmPtr->gauss();
122 0 : vertexY = sigmaVertexY * gauss;
123 0 : totalDev += gauss * gauss;
124 0 : }
125 0 : if (sigmaVertexZ > 0.) {
126 0 : gauss = rndmPtr->gauss();
127 0 : vertexZ = sigmaVertexZ * gauss;
128 0 : totalDev += gauss * gauss;
129 0 : }
130 0 : } while (totalDev > maxDevVertex * maxDevVertex);
131 :
132 : // Set beam collision time by a Gaussian.
133 0 : if (sigmaTime > 0.) {
134 0 : do gauss = rndmPtr->gauss();
135 0 : while (abs(gauss) > maxDevTime);
136 0 : vertexT = sigmaTime * gauss;
137 0 : }
138 :
139 : // Add offset to beam vertex.
140 0 : vertexX += offsetX;
141 0 : vertexY += offsetY;
142 0 : vertexZ += offsetZ;
143 0 : vertexT += offsetT;
144 0 : }
145 :
146 0 : }
147 :
148 : //==========================================================================
149 :
150 : } // end namespace Pythia8
|
