Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4AngularDistribution.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// hpw: done, but low quality at present.
27
28#include "globals.hh"
29#include "G4Log.hh"
30#include "G4SystemOfUnits.hh"
31#include "G4AngularDistribution.hh"
32#include "Randomize.hh"
33
34G4AngularDistribution::G4AngularDistribution(G4bool symmetrize)
35 : sym(symmetrize)
36{
37 // The following are parameters of the model - not to be confused with the PDG values!
38
39 mSigma = 0.55;
40 cmSigma = 1.20;
41 gSigma = 9.4;
42
43 mOmega = 0.783;
44 cmOmega = 0.808;
45 gOmega = 10.95;
46
47 mPion = 0.138;
48 cmPion = 0.51;
49 gPion = 7.27;
50
51 mNucleon = 0.938;
52
53 // Definition of constants for pion-Term (no s-dependence)
54
55 m42 = 4. * mNucleon * mNucleon;
56 mPion2 = mPion * mPion;
57 cmPion2 = cmPion * cmPion;
58 dPion1 = cmPion2-mPion2;
59 dPion2 = dPion1 * dPion1;
60 cm6gp = 1.5 * (cmPion2*cmPion2*cmPion2) * (gPion*gPion*gPion*gPion) * m42 * m42 / dPion2;
61
62 cPion_3 = -(cm6gp/3.);
63 cPion_2 = -(cm6gp * mPion2/dPion1);
64 cPion_1 = -(cm6gp * mPion2 * (2. * cmPion2 + mPion2) / dPion2);
65 cPion_m = -(cm6gp * cmPion2 * mPion2 / dPion2);
66 cPion_L = -(cm6gp * 2. * cmPion2 * mPion2 * (cmPion2 + mPion2) / dPion2 / dPion1);
67 cPion_0 = -(cPion_3 + cPion_2 + cPion_1 + cPion_m);
68
69 // Definition of constants for sigma-Term (no s-dependence)
70
71 G4double gSigmaSq = gSigma * gSigma;
72
73 mSigma2 = mSigma * mSigma;
74 cmSigma2 = cmSigma * cmSigma;
75 cmSigma4 = cmSigma2 * cmSigma2;
76 cmSigma6 = cmSigma2 * cmSigma4;
77 dSigma1 = m42 - cmSigma2;
78 dSigma2 = m42 - mSigma2;
79 dSigma3 = cmSigma2 - mSigma2;
80
81 G4double dSigma1Sq = dSigma1 * dSigma1;
82 G4double dSigma2Sq = dSigma2 * dSigma2;
83 G4double dSigma3Sq = dSigma3 * dSigma3;
84
85 cm2gs = 0.5 * cmSigma2 * gSigmaSq*gSigmaSq / dSigma3Sq;
86
87
88 cSigma_3 = -(cm2gs * dSigma1Sq / 3.);
89 cSigma_2 = -(cm2gs * cmSigma2 * dSigma1 * dSigma2 / dSigma3);
90 cSigma_1 = -(cm2gs * cmSigma4 * (2. * dSigma1 + dSigma2) * dSigma2 / dSigma3Sq);
91 cSigma_m = -(cm2gs * cmSigma6 * dSigma2Sq / mSigma2 / dSigma3Sq);
92 cSigma_L = -(cm2gs * cmSigma6 * dSigma2 * (dSigma1 + dSigma2) * 2. / (dSigma3 * dSigma3Sq));
93 cSigma_0 = -(cSigma_3 + cSigma_2 + cSigma_1 + cSigma_m);
94
95 // Definition of constants for omega-Term
96
97 G4double gOmegaSq = gOmega * gOmega;
98
99 mOmega2 = mOmega * mOmega;
100 cmOmega2 = cmOmega * cmOmega;
101 cmOmega4 = cmOmega2 * cmOmega2;
102 cmOmega6 = cmOmega2 * cmOmega4;
103 dOmega1 = m42 - cmOmega2;
104 dOmega2 = m42 - mOmega2;
105 dOmega3 = cmOmega2 - mOmega2;
106 sOmega1 = cmOmega2 + mOmega2;
107
108 G4double dOmega3Sq = dOmega3 * dOmega3;
109
110 cm2go = 0.5 * cmOmega2 * gOmegaSq * gOmegaSq / dOmega3Sq;
111
112 cOmega_3 = cm2go / 3.;
113 cOmega_2 = -(cm2go * cmOmega2 / dOmega3);
114 cOmega_1 = cm2go * cmOmega4 / dOmega3Sq;
115 cOmega_m = cm2go * cmOmega6 / (dOmega3Sq * mOmega2);
116 cOmega_L = -(cm2go * cmOmega6 * 4. / (dOmega3 * dOmega3Sq));
117
118 // Definition of constants for mix-Term
119
120 G4double fac1Tmp = (gSigma * gOmega * cmSigma2 * cmOmega2);
121 fac1 = -(fac1Tmp * fac1Tmp * m42);
122 dMix1 = cmOmega2 - cmSigma2;
123 dMix2 = cmOmega2 - mSigma2;
124 dMix3 = cmSigma2 - mOmega2;
125
126 G4double dMix1Sq = dMix1 * dMix1;
127 G4double dMix2Sq = dMix2 * dMix2;
128 G4double dMix3Sq = dMix3 * dMix3;
129
130 cMix_o1 = fac1 / (cmOmega2 * dMix1Sq * dMix2 * dOmega3);
131 cMix_s1 = fac1 / (cmSigma2 * dMix1Sq * dMix3 * dSigma3);
132 cMix_Omega = fac1 / (dOmega3Sq * dMix3Sq * (mOmega2 - mSigma2));
133 cMix_sm = fac1 / (dSigma3Sq * dMix2Sq * (mSigma2 - mOmega2));
134 fac2 = (-fac1) / (dMix1*dMix1Sq * dOmega3Sq * dMix2Sq);
135 fac3 = (-fac1) / (dMix1*dMix1Sq * dSigma3Sq * dMix3Sq);
136
137 cMix_oLc = fac2 * (3. * cmOmega2*cmOmega4 - cmOmega4 * cmSigma2
138 - 2. * cmOmega4 * mOmega2 - 2. * cmOmega4 * mSigma2
139 + cmOmega2 * mOmega2 * mSigma2 + cmSigma2 * mOmega2 * mSigma2
140 - 4. * cmOmega4 * m42 + 2. * cmOmega2 * cmSigma2 * m42
141 + 3. * cmOmega2 * mOmega2 * m42 - cmSigma2 * mOmega2 * m42
142 + 3. * cmOmega2 * mSigma2 * m42 - cmSigma2 * mSigma2 * m42
143 - 2. * mOmega2 * mSigma2 * m42);
144
145 cMix_oLs = fac2 * (8. * cmOmega4 - 4. * cmOmega2 * cmSigma2
146 - 6. * cmOmega2 * mOmega2 + 2. * cmSigma2 * mOmega2
147 - 6. * cmOmega2 * mSigma2 + 2. * cmSigma2 * mSigma2
148 + 4. * mOmega2 * mSigma2);
149
150 cMix_sLc = fac3 * (cmOmega2 * cmSigma4 - 3. * cmSigma6
151 + 2. * cmSigma4 * mOmega2 + 2. * cmSigma4 * mSigma2
152 - cmOmega2 * mOmega2 * mSigma2 - cmSigma2 * mOmega2 * mSigma2
153 - 2. * cmOmega2 * cmSigma2 * m42 + 4. * cmSigma4 * m42
154 + cmOmega2 * mOmega2 * m42 - 3. * cmSigma2 * mOmega2 * m42
155 + cmOmega2 * mSigma2 * m42 - 3. * cmSigma2 * mSigma2 * m42
156 + 2. * mOmega2 * mSigma2 * m42);
157
158 cMix_sLs = fac3 * (4. * cmOmega2 * cmSigma2 - 8. * cmSigma4
159 - 2. * cmOmega2 * mOmega2 + 6. * cmSigma2 * mOmega2
160 - 2. * cmOmega2 * mSigma2 + 6. * cmSigma2 * mSigma2
161 - 4. * mOmega2 * mSigma2);
162}
163
164
165G4AngularDistribution::~
166G4AngularDistribution()
167{ }
168
169
170G4double G4AngularDistribution::CosTheta(G4double S, G4double m_1, G4double m_2) const
171{
172 G4double random = G4UniformRand();
173 G4double dCosTheta = 2.;
174 G4double cosTheta = -1.;
175
176 // For jmax=12 the accuracy is better than 0.1 degree
177 G4int jMax = 12;
178
179 for (G4int j = 1; j <= jMax; ++j)
180 {
181 // Accuracy is 2^-jmax
182 dCosTheta *= 0.5;
183 G4double cosTh = cosTheta + dCosTheta;
184 if(DifferentialCrossSection(S, m_1, m_2, cosTh) <= random) cosTheta = cosTh;
185 }
186
187 // Randomize in final interval in order to avoid discrete angles
188 cosTheta += G4UniformRand() * dCosTheta;
189
190
191 if (cosTheta > 1. || cosTheta < -1.)
192 throw G4HadronicException(__FILE__, __LINE__, "G4AngularDistribution::CosTheta - std::cos(theta) outside allowed range");
193
194 return cosTheta;
195}
196
197
198G4double G4AngularDistribution::DifferentialCrossSection(G4double sIn, G4double m_1, G4double m_2,
199 G4double cosTheta) const
200{
201// local calculus is in GeV, ie. normalize input
202 sIn = sIn/sqr(GeV)+m42/2.;
203 m_1 = m_1/GeV;
204 m_2 = m_2/GeV;
205// G4cout << "Here we go"<<sIn << " "<<m1 << " " << m2 <<" " m42<< G4endl;
206// scaling from masses other than p,p.
207 G4double S = sIn - (m_1+m_2) * (m_1+m_2) + m42;
208 G4double tMax = S - m42;
209 G4double tp = 0.5 * (cosTheta + 1.) * tMax;
210 G4double twoS = 2. * S;
211
212 // Define s-dependent stuff for omega-Term
213 G4double brak1 = (twoS-m42) * (twoS-m42);
214 G4double bOmega_3 = cOmega_3 * (-2. * cmOmega4 - 2. * cmOmega2 * twoS - brak1);
215 G4double bOmega_2 = cOmega_2 * ( 2. * cmOmega2 * mOmega2 + sOmega1 * twoS + brak1);
216 G4double bOmega_1 = cOmega_1 * (-4. * cmOmega2 * mOmega2
217 - 2. * mOmega2*mOmega2
218 - 2. * (cmOmega2 + 2 * mOmega2) * twoS
219 - 3. * brak1);
220 G4double bOmega_m = cOmega_m * (-2. * mOmega2*mOmega2 - 2. * mOmega2 * twoS - brak1);
221 G4double bOmega_L = cOmega_L * (sOmega1 * mOmega2 + (cmOmega2 + 3. * mOmega2) * S + brak1);
222 G4double bOmega_0 = -(bOmega_3 + bOmega_2 + bOmega_1 + bOmega_m);
223
224 // Define s-dependent stuff for mix-Term
225 G4double bMix_o1 = cMix_o1 * (dOmega1 - twoS);
226 G4double bMix_s1 = cMix_s1 * (dSigma1 - twoS);
227 G4double bMix_Omega = cMix_Omega * (dOmega2 - twoS);
228 G4double bMix_sm = cMix_sm * (dSigma2 - twoS);
229 G4double bMix_oL = cMix_oLc + cMix_oLs * S;
230 G4double bMix_sL = cMix_sLc + cMix_sLs * S;
231
232 G4double t1_Pion = 1. / (1. + tMax / cmPion2);
233 G4double t2_Pion = 1. + tMax / mPion2;
234 G4double t1_Sigma = 1. / (1. + tMax / cmSigma2);
235 G4double t2_Sigma = 1. + tMax / mSigma2;
236 G4double t1_Omega = 1. / (1. + tMax / cmOmega2);
237 G4double t2_Omega = 1. + tMax / mOmega2;
238
239 G4double norm = Cross(t1_Pion, t1_Sigma, t1_Omega,
240 t2_Pion, t2_Sigma, t2_Omega,
241 bMix_o1, bMix_s1, bMix_Omega,
242 bMix_sm, bMix_oL, bMix_sL,
243 bOmega_0, bOmega_1, bOmega_2,
244 bOmega_3, bOmega_m, bOmega_L);
245
246 t1_Pion = 1. / (1. + tp / cmPion2);
247 t2_Pion = 1. + tp / mPion2;
248 t1_Sigma = 1. / (1. + tp / cmSigma2);
249 t2_Sigma = 1. + tp / mSigma2;
250 t1_Omega = 1. / (1. + tp / cmOmega2);
251 t2_Omega = 1. + tp / mOmega2;
252
253 G4double dSigma;
254 if (sym)
255 {
256 G4double to;
257 norm = 2. * norm;
258 to = tMax - tp;
259 G4double t3_Pion = 1. / (1. + to / cmPion2);
260 G4double t4_Pion = 1. + to / mPion2;
261 G4double t3_Sigma = 1. / (1. + to / cmSigma2);
262 G4double t4_Sigma = 1. + to / mSigma2;
263 G4double t3_Omega = 1. / (1. + to / cmOmega2);
264 G4double t4_Omega = 1. + to / mOmega2;
265
266 dSigma = ( Cross(t1_Pion, t1_Sigma, t1_Omega,
267 t2_Pion,t2_Sigma, t2_Omega,
268 bMix_o1, bMix_s1, bMix_Omega,
269 bMix_sm, bMix_oL, bMix_sL,
270 bOmega_0, bOmega_1, bOmega_2,
271 bOmega_3, bOmega_m, bOmega_L) -
272 Cross(t3_Pion,t3_Sigma, t3_Omega,
273 t4_Pion, t4_Sigma, t4_Omega,
274 bMix_o1, bMix_s1, bMix_Omega,
275 bMix_sm, bMix_oL, bMix_sL,
276 bOmega_0, bOmega_1, bOmega_2,
277 bOmega_3, bOmega_m, bOmega_L) )
278 / norm + 0.5;
279 }
280 else
281 {
282 dSigma = Cross(t1_Pion, t1_Sigma, t1_Omega,
283 t2_Pion, t2_Sigma, t2_Omega,
284 bMix_o1, bMix_s1, bMix_Omega,
285 bMix_sm, bMix_oL, bMix_sL,
286 bOmega_0, bOmega_1, bOmega_2,
287 bOmega_3, bOmega_m, bOmega_L)
288 / norm;
289 }
290
291 return dSigma;
292}
293
294
295G4double G4AngularDistribution::Cross(G4double tpPion,
296 G4double tpSigma,
297 G4double tpOmega,
298 G4double tmPion,
299 G4double tmSigma,
300 G4double tmOmega,
301 G4double bMix_o1,
302 G4double bMix_s1,
303 G4double bMix_Omega,
304 G4double bMix_sm,
305 G4double bMix_oL,
306 G4double bMix_sL,
307 G4double bOmega_0,
308 G4double bOmega_1,
309 G4double bOmega_2,
310 G4double bOmega_3,
311 G4double bOmega_m,
312 G4double bOmega_L) const
313{
314 G4double cross = 0;
315 // Pion
316 cross += ((cPion_3 * tpPion + cPion_2) * tpPion + cPion_1) * tpPion + cPion_m/tmPion + cPion_0 + cPion_L * G4Log(tpPion*tmPion);
317// G4cout << "cross1 "<< cross<<G4endl;
318 // Sigma
319 cross += ((cSigma_3 * tpSigma + cSigma_2) * tpSigma + cSigma_1) * tpSigma + cSigma_m/tmSigma + cSigma_0 + cSigma_L * G4Log(tpSigma*tmSigma);
320// G4cout << "cross2 "<< cross<<G4endl;
321 // Omega
322 cross += ((bOmega_3 * tpOmega + bOmega_2) * tpOmega + bOmega_1) * tpOmega + bOmega_m/tmOmega + bOmega_0 + bOmega_L * G4Log(tpOmega*tmOmega)
323 // Mix
324 + bMix_o1 * (tpOmega - 1.)
325 + bMix_s1 * (tpSigma - 1.)
326 + bMix_Omega * G4Log(tmOmega)
327 + bMix_sm * G4Log(tmSigma)
328 + bMix_oL * G4Log(tpOmega)
329 + bMix_sL * G4Log(tpSigma);
330/* G4cout << "cross3 "<< cross<<" "
331 <<bMix_o1<<" "
332 <<bMix_s1<<" "
333 <<bMix_Omega<<" "
334 <<bMix_sm<<" "
335 <<bMix_oL<<" "
336 <<bMix_sL<<" "
337 <<tpOmega<<" "
338 <<tpSigma<<" "
339 <<tmOmega<<" "
340 <<tmSigma<<" "
341 <<tpOmega<<" "
342 <<tpSigma
343 <<G4endl;
344*/
345 return cross;
346
347}
S
G4double S(G4double temp)
Definition: G4DNAElectronHoleRecombination.cc:76
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4AngularDistribution::G4AngularDistribution
G4AngularDistribution(G4bool symmetrize)
Definition: G4AngularDistribution.cc:34
G4AngularDistribution::CosTheta
virtual G4double CosTheta(G4double s, G4double m1, G4double m2) const
Definition: G4AngularDistribution.cc:170
G4AngularDistribution::Cross
G4double Cross(G4double tpPion, G4double tpSigma, G4double tpOmega, G4double tmPion, G4double tmSigma, G4double tmOmega, G4double bMix_o1, G4double bMix_s1, G4double bMix_Omega, G4double bMix_sm, G4double bMix_oL, G4double bMix_sL, G4double bOmega_0, G4double bOmega_1, G4double bOmega_2, G4double bOmega_3, G4double bOmega_m, G4double bOmega_L) const
Definition: G4AngularDistribution.cc:295
G4AngularDistribution::DifferentialCrossSection
G4double DifferentialCrossSection(G4double sIn, G4double m1, G4double m2, G4double cosTheta) const
Definition: G4AngularDistribution.cc:198
sqr
T sqr(const T &x)
Definition: templates.hh:128