| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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
34 G4AngularDistribution::G4AngularDistribution(G
35 : sym(symmetrize)
36 {
37 // The following are parameters of the model
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
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)
61
62 cPion_3 = -(cm6gp/3.);
63 cPion_2 = -(cm6gp * mPion2/dPion1);
64 cPion_1 = -(cm6gp * mPion2 * (2. * cmPio
65 cPion_m = -(cm6gp * cmPion2 * mPion2 / d
66 cPion_L = -(cm6gp * 2. * cmPion2 * mPion
67 cPion_0 = -(cPion_3 + cPion_2 + cPion_1
68
69 // Definition of constants for sigma-Term
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*gSigma
86
87
88 cSigma_3 = -(cm2gs * dSigma1Sq / 3.);
89 cSigma_2 = -(cm2gs * cmSigma2 * dSigma1
90 cSigma_1 = -(cm2gs * cmSigma4 * (2. * dS
91 cSigma_m = -(cm2gs * cmSigma6 * dSigma2S
92 cSigma_L = -(cm2gs * cmSigma6 * dSigma2
93 cSigma_0 = -(cSigma_3 + cSigma_2 + cSigm
94
95 // Definition of constants for omega-Ter
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 * gOme
111
112 cOmega_3 = cm2go / 3.;
113 cOmega_2 = -(cm2go * cmOmega2 / dOmega3)
114 cOmega_1 = cm2go * cmOmega4 / dOmega3Sq
115 cOmega_m = cm2go * cmOmega6 / (dOmega3S
116 cOmega_L = -(cm2go * cmOmega6 * 4. / (dO
117
118 // Definition of constants for mix-Term
119
120 G4double fac1Tmp = (gSigma * gOmega * cm
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
131 cMix_s1 = fac1 / (cmSigma2 * dMix1Sq
132 cMix_Omega = fac1 / (dOmega3Sq * dMix3Sq
133 cMix_sm = fac1 / (dSigma3Sq * dMix2Sq
134 fac2 = (-fac1) / (dMix1*dMix1Sq * dOmega
135 fac3 = (-fac1) / (dMix1*dMix1Sq * dSigma
136
137 cMix_oLc = fac2 * (3. * cmOmega2*cmOmega
138 - 2. * cmOmega4 * mOmega2
139 + cmOmega2 * mOmega2 * mSigma2
140 - 4. * cmOmega4 * m42
141 + 3. * cmOmega2 * mOmega2 * m42
142 + 3. * cmOmega2 * mSigma2 * m42
143 - 2. * mOmega2 * mSigma2 * m42);
144
145 cMix_oLs = fac2 * (8. * cmOmega4
146 - 6. * cmOmega2 * mOmega2 + 2.
147 - 6. * cmOmega2 * mSigma2 + 2.
148 + 4. * mOmega2 * mSigma2);
149
150 cMix_sLc = fac3 * (cmOmega2 * cmSigma4
151 + 2. * cmSigma4 * mOmega2 + 2.
152 - cmOmega2 * mOmega2 * mSigma2 - cmS
153 - 2. * cmOmega2 * cmSigma2 * m42 + 4.
154 + cmOmega2 * mOmega2 * m42 - 3.
155 + cmOmega2 * mSigma2 * m42 - 3.
156 + 2. * mOmega2 * mSigma2 * m42);
157
158 cMix_sLs = fac3 * (4. * cmOmega2 * cmSig
159 - 2. * cmOmega2 * mOmega2
160 - 2. * cmOmega2 * mSigma2
161 - 4. * mOmega2 * mSigma2);
162 }
163
164
165 G4AngularDistribution::~
166 G4AngularDistribution()
167 { }
168
169
170 G4double G4AngularDistribution::CosTheta(G4dou
171 {
172 G4double random = G4UniformRand();
173 G4double dCosTheta = 2.;
174 G4double cosTheta = -1.;
175
176 // For jmax=12 the accuracy is better than
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,
185 }
186
187 // Randomize in final interval in order to
188 cosTheta += G4UniformRand() * dCosTheta;
189
190
191 if (cosTheta > 1. || cosTheta < -1.)
192 throw G4HadronicException(__FILE__, __LIN
193
194 return cosTheta;
195 }
196
197
198 G4double G4AngularDistribution::DifferentialCr
199 G4double cosTheta) const
200 {
201 // local calculus is in GeV, ie. normalize inp
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 <<
206 // scaling from masses other than p,p.
207 G4double S = sIn - (m_1+m_2) * (m_1+m_2) + m
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. * cmOmeg
215 G4double bOmega_2 = cOmega_2 * ( 2. * cmOmeg
216 G4double bOmega_1 = cOmega_1 * (-4. * cmOmeg
217 - 2. * mOmega2*mOmega2
218 - 2. * (cmOmega2 + 2 * mOmega2) * twoS
219 - 3. * brak1);
220 G4double bOmega_m = cOmega_m * (-2. * mOmega
221 G4double bOmega_L = cOmega_L * (sOmega1 * mO
222 G4double bOmega_0 = -(bOmega_3 + bOmega_2 +
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
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 / cmSigm
235 G4double t2_Sigma = 1. + tMax / mSigma2;
236 G4double t1_Omega = 1. / (1. + tMax / cmOmeg
237 G4double t2_Omega = 1. + tMax / mOmega2;
238
239 G4double norm = Cross(t1_Pion, t1_Sigma, t1_
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 / cmPio
260 G4double t4_Pion = 1. + to / mPion2;
261 G4double t3_Sigma = 1. / (1. + to / cmSi
262 G4double t4_Sigma = 1. + to / mSigma2;
263 G4double t3_Omega = 1. / (1. + to / cmOm
264 G4double t4_Omega = 1. + to / mOmega2;
265
266 dSigma = ( Cross(t1_Pion, t1_Sigma, t1_O
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_Ome
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
295 G4double G4AngularDistribution::Cross(G4double
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) *
317 // G4cout << "cross1 "<< cross<<G4endl;
318 // Sigma
319 cross += ((cSigma_3 * tpSigma + cSigma_2)
320 // G4cout << "cross2 "<< cross<<G4endl;
321 // Omega
322 cross += ((bOmega_3 * tpOmega + bOmega_2)
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 }
348