StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
BeamShape.cc
1 // BeamShape.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2012 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