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 "G4AngularDistribution.hh" |
---|
30 | #include "Randomize.hh" |
---|
31 | |
---|
32 | G4AngularDistribution::G4AngularDistribution(G4bool symmetrize) |
---|
33 | : sym(symmetrize) |
---|
34 | { |
---|
35 | // The following are parameters of the model - not to be confused with the PDG values! |
---|
36 | |
---|
37 | mSigma = 0.55; |
---|
38 | cmSigma = 1.20; |
---|
39 | gSigma = 9.4; |
---|
40 | |
---|
41 | mOmega = 0.783; |
---|
42 | cmOmega = 0.808; |
---|
43 | gOmega = 10.95; |
---|
44 | |
---|
45 | mPion = 0.138; |
---|
46 | cmPion = 0.51; |
---|
47 | gPion = 7.27; |
---|
48 | |
---|
49 | mNucleon = 0.938; |
---|
50 | |
---|
51 | // Definition of constants for pion-Term (no s-dependence) |
---|
52 | |
---|
53 | m42 = 4. * mNucleon * mNucleon; |
---|
54 | mPion2 = mPion * mPion; |
---|
55 | cmPion2 = cmPion * cmPion; |
---|
56 | dPion1 = cmPion2-mPion2; |
---|
57 | dPion2 = dPion1 * dPion1; |
---|
58 | cm6gp = 1.5 * (cmPion2*cmPion2*cmPion2) * (gPion*gPion*gPion*gPion) * m42 * m42 / dPion2; |
---|
59 | |
---|
60 | cPion_3 = -(cm6gp/3.); |
---|
61 | cPion_2 = -(cm6gp * mPion2/dPion1); |
---|
62 | cPion_1 = -(cm6gp * mPion2 * (2. * cmPion2 + mPion2) / dPion2); |
---|
63 | cPion_m = -(cm6gp * cmPion2 * mPion2 / dPion2); |
---|
64 | cPion_L = -(cm6gp * 2. * cmPion2 * mPion2 * (cmPion2 + mPion2) / dPion2 / dPion1); |
---|
65 | cPion_0 = -(cPion_3 + cPion_2 + cPion_1 + cPion_m); |
---|
66 | |
---|
67 | // Definition of constants for sigma-Term (no s-dependence) |
---|
68 | |
---|
69 | G4double gSigmaSq = gSigma * gSigma; |
---|
70 | |
---|
71 | mSigma2 = mSigma * mSigma; |
---|
72 | cmSigma2 = cmSigma * cmSigma; |
---|
73 | cmSigma4 = cmSigma2 * cmSigma2; |
---|
74 | cmSigma6 = cmSigma2 * cmSigma4; |
---|
75 | dSigma1 = m42 - cmSigma2; |
---|
76 | dSigma2 = m42 - mSigma2; |
---|
77 | dSigma3 = cmSigma2 - mSigma2; |
---|
78 | |
---|
79 | G4double dSigma1Sq = dSigma1 * dSigma1; |
---|
80 | G4double dSigma2Sq = dSigma2 * dSigma2; |
---|
81 | G4double dSigma3Sq = dSigma3 * dSigma3; |
---|
82 | |
---|
83 | cm2gs = 0.5 * cmSigma2 * gSigmaSq*gSigmaSq / dSigma3Sq; |
---|
84 | |
---|
85 | |
---|
86 | cSigma_3 = -(cm2gs * dSigma1Sq / 3.); |
---|
87 | cSigma_2 = -(cm2gs * cmSigma2 * dSigma1 * dSigma2 / dSigma3); |
---|
88 | cSigma_1 = -(cm2gs * cmSigma4 * (2. * dSigma1 + dSigma2) * dSigma2 / dSigma3Sq); |
---|
89 | cSigma_m = -(cm2gs * cmSigma6 * dSigma2Sq / mSigma2 / dSigma3Sq); |
---|
90 | cSigma_L = -(cm2gs * cmSigma6 * dSigma2 * (dSigma1 + dSigma2) * 2. / (dSigma3 * dSigma3Sq)); |
---|
91 | cSigma_0 = -(cSigma_3 + cSigma_2 + cSigma_1 + cSigma_m); |
---|
92 | |
---|
93 | // Definition of constants for omega-Term |
---|
94 | |
---|
95 | G4double gOmegaSq = gOmega * gOmega; |
---|
96 | |
---|
97 | mOmega2 = mOmega * mOmega; |
---|
98 | cmOmega2 = cmOmega * cmOmega; |
---|
99 | cmOmega4 = cmOmega2 * cmOmega2; |
---|
100 | cmOmega6 = cmOmega2 * cmOmega4; |
---|
101 | dOmega1 = m42 - cmOmega2; |
---|
102 | dOmega2 = m42 - mOmega2; |
---|
103 | dOmega3 = cmOmega2 - mOmega2; |
---|
104 | sOmega1 = cmOmega2 + mOmega2; |
---|
105 | |
---|
106 | G4double dOmega3Sq = dOmega3 * dOmega3; |
---|
107 | |
---|
108 | cm2go = 0.5 * cmOmega2 * gOmegaSq * gOmegaSq / dOmega3Sq; |
---|
109 | |
---|
110 | cOmega_3 = cm2go / 3.; |
---|
111 | cOmega_2 = -(cm2go * cmOmega2 / dOmega3); |
---|
112 | cOmega_1 = cm2go * cmOmega4 / dOmega3Sq; |
---|
113 | cOmega_m = cm2go * cmOmega6 / (dOmega3Sq * mOmega2); |
---|
114 | cOmega_L = -(cm2go * cmOmega6 * 4. / (dOmega3 * dOmega3Sq)); |
---|
115 | |
---|
116 | // Definition of constants for mix-Term |
---|
117 | |
---|
118 | G4double fac1Tmp = (gSigma * gOmega * cmSigma2 * cmOmega2); |
---|
119 | fac1 = -(fac1Tmp * fac1Tmp * m42); |
---|
120 | dMix1 = cmOmega2 - cmSigma2; |
---|
121 | dMix2 = cmOmega2 - mSigma2; |
---|
122 | dMix3 = cmSigma2 - mOmega2; |
---|
123 | |
---|
124 | G4double dMix1Sq = dMix1 * dMix1; |
---|
125 | G4double dMix2Sq = dMix2 * dMix2; |
---|
126 | G4double dMix3Sq = dMix3 * dMix3; |
---|
127 | |
---|
128 | cMix_o1 = fac1 / (cmOmega2 * dMix1Sq * dMix2 * dOmega3); |
---|
129 | cMix_s1 = fac1 / (cmSigma2 * dMix1Sq * dMix3 * dSigma3); |
---|
130 | cMix_Omega = fac1 / (dOmega3Sq * dMix3Sq * (mOmega2 - mSigma2)); |
---|
131 | cMix_sm = fac1 / (dSigma3Sq * dMix2Sq * (mSigma2 - mOmega2)); |
---|
132 | fac2 = (-fac1) / (dMix1*dMix1Sq * dOmega3Sq * dMix2Sq); |
---|
133 | fac3 = (-fac1) / (dMix1*dMix1Sq * dSigma3Sq * dMix3Sq); |
---|
134 | |
---|
135 | cMix_oLc = fac2 * (3. * cmOmega2*cmOmega4 - cmOmega4 * cmSigma2 |
---|
136 | - 2. * cmOmega4 * mOmega2 - 2. * cmOmega4 * mSigma2 |
---|
137 | + cmOmega2 * mOmega2 * mSigma2 + cmSigma2 * mOmega2 * mSigma2 |
---|
138 | - 4. * cmOmega4 * m42 + 2. * cmOmega2 * cmSigma2 * m42 |
---|
139 | + 3. * cmOmega2 * mOmega2 * m42 - cmSigma2 * mOmega2 * m42 |
---|
140 | + 3. * cmOmega2 * mSigma2 * m42 - cmSigma2 * mSigma2 * m42 |
---|
141 | - 2. * mOmega2 * mSigma2 * m42); |
---|
142 | |
---|
143 | cMix_oLs = fac2 * (8. * cmOmega4 - 4. * cmOmega2 * cmSigma2 |
---|
144 | - 6. * cmOmega2 * mOmega2 + 2. * cmSigma2 * mOmega2 |
---|
145 | - 6. * cmOmega2 * mSigma2 + 2. * cmSigma2 * mSigma2 |
---|
146 | + 4. * mOmega2 * mSigma2); |
---|
147 | |
---|
148 | cMix_sLc = fac3 * (cmOmega2 * cmSigma4 - 3. * cmSigma6 |
---|
149 | + 2. * cmSigma4 * mOmega2 + 2. * cmSigma4 * mSigma2 |
---|
150 | - cmOmega2 * mOmega2 * mSigma2 - cmSigma2 * mOmega2 * mSigma2 |
---|
151 | - 2. * cmOmega2 * cmSigma2 * m42 + 4. * cmSigma4 * m42 |
---|
152 | + cmOmega2 * mOmega2 * m42 - 3. * cmSigma2 * mOmega2 * m42 |
---|
153 | + cmOmega2 * mSigma2 * m42 - 3. * cmSigma2 * mSigma2 * m42 |
---|
154 | + 2. * mOmega2 * mSigma2 * m42); |
---|
155 | |
---|
156 | cMix_sLs = fac3 * (4. * cmOmega2 * cmSigma2 - 8. * cmSigma4 |
---|
157 | - 2. * cmOmega2 * mOmega2 + 6. * cmSigma2 * mOmega2 |
---|
158 | - 2. * cmOmega2 * mSigma2 + 6. * cmSigma2 * mSigma2 |
---|
159 | - 4. * mOmega2 * mSigma2); |
---|
160 | } |
---|
161 | |
---|
162 | |
---|
163 | G4AngularDistribution::~ |
---|
164 | G4AngularDistribution() |
---|
165 | { } |
---|
166 | |
---|
167 | |
---|
168 | G4double G4AngularDistribution::CosTheta(G4double s, G4double m1, G4double m2) const |
---|
169 | { |
---|
170 | G4double random = G4UniformRand(); |
---|
171 | G4double dCosTheta = 2.; |
---|
172 | G4double cosTheta = -1.; |
---|
173 | |
---|
174 | // For jmax=12 the accuracy is better than 0.1 degree |
---|
175 | G4int jMax = 12; |
---|
176 | |
---|
177 | for (G4int j = 1; j <= jMax; ++j) |
---|
178 | { |
---|
179 | // Accuracy is 2^-jmax |
---|
180 | dCosTheta *= 0.5; |
---|
181 | G4double cosTh = cosTheta + dCosTheta; |
---|
182 | if(DifferentialCrossSection(s, m1, m2, cosTh) <= random) cosTheta = cosTh; |
---|
183 | } |
---|
184 | |
---|
185 | // Randomize in final interval in order to avoid discrete angles |
---|
186 | cosTheta += G4UniformRand() * dCosTheta; |
---|
187 | |
---|
188 | |
---|
189 | if (cosTheta > 1. || cosTheta < -1.) |
---|
190 | throw G4HadronicException(__FILE__, __LINE__, "G4AngularDistribution::CosTheta - std::cos(theta) outside allowed range"); |
---|
191 | |
---|
192 | return cosTheta; |
---|
193 | } |
---|
194 | |
---|
195 | |
---|
196 | G4double G4AngularDistribution::DifferentialCrossSection(G4double sIn, G4double m1, G4double m2, |
---|
197 | G4double cosTheta) const |
---|
198 | { |
---|
199 | // local calculus is in GeV, ie. normalize input |
---|
200 | sIn = sIn/sqr(GeV)+m42/2.; |
---|
201 | m1 = m1/GeV; |
---|
202 | m2 = m2/GeV; |
---|
203 | // G4cout << "Here we go"<<sIn << " "<<m1 << " " << m2 <<" " m42<< G4endl; |
---|
204 | // scaling from masses other than p,p. |
---|
205 | G4double s = sIn - (m1+m2) * (m1+m2) + m42; |
---|
206 | G4double tMax = s - m42; |
---|
207 | G4double tp = 0.5 * (cosTheta + 1.) * tMax; |
---|
208 | G4double twoS = 2. * s; |
---|
209 | |
---|
210 | // Define s-dependent stuff for omega-Term |
---|
211 | G4double brak1 = (twoS-m42) * (twoS-m42); |
---|
212 | G4double bOmega_3 = cOmega_3 * (-2. * cmOmega4 - 2. * cmOmega2 * twoS - brak1); |
---|
213 | G4double bOmega_2 = cOmega_2 * ( 2. * cmOmega2 * mOmega2 + sOmega1 * twoS + brak1); |
---|
214 | G4double bOmega_1 = cOmega_1 * (-4. * cmOmega2 * mOmega2 |
---|
215 | - 2. * mOmega2*mOmega2 |
---|
216 | - 2. * (cmOmega2 + 2 * mOmega2) * twoS |
---|
217 | - 3. * brak1); |
---|
218 | G4double bOmega_m = cOmega_m * (-2. * mOmega2*mOmega2 - 2. * mOmega2 * twoS - brak1); |
---|
219 | G4double bOmega_L = cOmega_L * (sOmega1 * mOmega2 + (cmOmega2 + 3. * mOmega2) * s + brak1); |
---|
220 | G4double bOmega_0 = -(bOmega_3 + bOmega_2 + bOmega_1 + bOmega_m); |
---|
221 | |
---|
222 | // Define s-dependent stuff for mix-Term |
---|
223 | G4double bMix_o1 = cMix_o1 * (dOmega1 - twoS); |
---|
224 | G4double bMix_s1 = cMix_s1 * (dSigma1 - twoS); |
---|
225 | G4double bMix_Omega = cMix_Omega * (dOmega2 - twoS); |
---|
226 | G4double bMix_sm = cMix_sm * (dSigma2 - twoS); |
---|
227 | G4double bMix_oL = cMix_oLc + cMix_oLs * s; |
---|
228 | G4double bMix_sL = cMix_sLc + cMix_sLs * s; |
---|
229 | |
---|
230 | G4double t1_Pion = 1. / (1. + tMax / cmPion2); |
---|
231 | G4double t2_Pion = 1. + tMax / mPion2; |
---|
232 | G4double t1_Sigma = 1. / (1. + tMax / cmSigma2); |
---|
233 | G4double t2_Sigma = 1. + tMax / mSigma2; |
---|
234 | G4double t1_Omega = 1. / (1. + tMax / cmOmega2); |
---|
235 | G4double t2_Omega = 1. + tMax / mOmega2; |
---|
236 | |
---|
237 | G4double norm = Cross(t1_Pion, t1_Sigma, t1_Omega, |
---|
238 | t2_Pion, t2_Sigma, t2_Omega, |
---|
239 | bMix_o1, bMix_s1, bMix_Omega, |
---|
240 | bMix_sm, bMix_oL, bMix_sL, |
---|
241 | bOmega_0, bOmega_1, bOmega_2, |
---|
242 | bOmega_3, bOmega_m, bOmega_L); |
---|
243 | |
---|
244 | t1_Pion = 1. / (1. + tp / cmPion2); |
---|
245 | t2_Pion = 1. + tp / mPion2; |
---|
246 | t1_Sigma = 1. / (1. + tp / cmSigma2); |
---|
247 | t2_Sigma = 1. + tp / mSigma2; |
---|
248 | t1_Omega = 1. / (1. + tp / cmOmega2); |
---|
249 | t2_Omega = 1. + tp / mOmega2; |
---|
250 | |
---|
251 | G4double dSigma; |
---|
252 | if (sym) |
---|
253 | { |
---|
254 | G4double to; |
---|
255 | norm = 2. * norm; |
---|
256 | to = tMax - tp; |
---|
257 | G4double t3_Pion = 1. / (1. + to / cmPion2); |
---|
258 | G4double t4_Pion = 1. + to / mPion2; |
---|
259 | G4double t3_Sigma = 1. / (1. + to / cmSigma2); |
---|
260 | G4double t4_Sigma = 1. + to / mSigma2; |
---|
261 | G4double t3_Omega = 1. / (1. + to / cmOmega2); |
---|
262 | G4double t4_Omega = 1. + to / mOmega2; |
---|
263 | |
---|
264 | dSigma = ( Cross(t1_Pion, t1_Sigma, t1_Omega, |
---|
265 | t2_Pion,t2_Sigma, t2_Omega, |
---|
266 | bMix_o1, bMix_s1, bMix_Omega, |
---|
267 | bMix_sm, bMix_oL, bMix_sL, |
---|
268 | bOmega_0, bOmega_1, bOmega_2, |
---|
269 | bOmega_3, bOmega_m, bOmega_L) - |
---|
270 | Cross(t3_Pion,t3_Sigma, t3_Omega, |
---|
271 | t4_Pion, t4_Sigma, t4_Omega, |
---|
272 | bMix_o1, bMix_s1, bMix_Omega, |
---|
273 | bMix_sm, bMix_oL, bMix_sL, |
---|
274 | bOmega_0, bOmega_1, bOmega_2, |
---|
275 | bOmega_3, bOmega_m, bOmega_L) ) |
---|
276 | / norm + 0.5; |
---|
277 | } |
---|
278 | else |
---|
279 | { |
---|
280 | dSigma = Cross(t1_Pion, t1_Sigma, t1_Omega, |
---|
281 | t2_Pion, t2_Sigma, t2_Omega, |
---|
282 | bMix_o1, bMix_s1, bMix_Omega, |
---|
283 | bMix_sm, bMix_oL, bMix_sL, |
---|
284 | bOmega_0, bOmega_1, bOmega_2, |
---|
285 | bOmega_3, bOmega_m, bOmega_L) |
---|
286 | / norm; |
---|
287 | } |
---|
288 | |
---|
289 | return dSigma; |
---|
290 | } |
---|
291 | |
---|
292 | |
---|
293 | G4double G4AngularDistribution::Cross(G4double tpPion, |
---|
294 | G4double tpSigma, |
---|
295 | G4double tpOmega, |
---|
296 | G4double tmPion, |
---|
297 | G4double tmSigma, |
---|
298 | G4double tmOmega, |
---|
299 | G4double bMix_o1, |
---|
300 | G4double bMix_s1, |
---|
301 | G4double bMix_Omega, |
---|
302 | G4double bMix_sm, |
---|
303 | G4double bMix_oL, |
---|
304 | G4double bMix_sL, |
---|
305 | G4double bOmega_0, |
---|
306 | G4double bOmega_1, |
---|
307 | G4double bOmega_2, |
---|
308 | G4double bOmega_3, |
---|
309 | G4double bOmega_m, |
---|
310 | G4double bOmega_L) const |
---|
311 | { |
---|
312 | G4double cross = 0; |
---|
313 | // Pion |
---|
314 | cross += ((cPion_3 * tpPion + cPion_2) * tpPion + cPion_1) * tpPion + cPion_m/tmPion + cPion_0 + cPion_L * std::log(tpPion*tmPion); |
---|
315 | // G4cout << "cross1 "<< cross<<G4endl; |
---|
316 | // Sigma |
---|
317 | cross += ((cSigma_3 * tpSigma + cSigma_2) * tpSigma + cSigma_1) * tpSigma + cSigma_m/tmSigma + cSigma_0 + cSigma_L * std::log(tpSigma*tmSigma); |
---|
318 | // G4cout << "cross2 "<< cross<<G4endl; |
---|
319 | // Omega |
---|
320 | cross += ((bOmega_3 * tpOmega + bOmega_2) * tpOmega + bOmega_1) * tpOmega + bOmega_m/tmOmega + bOmega_0 + bOmega_L * std::log(tpOmega*tmOmega) |
---|
321 | // Mix |
---|
322 | + bMix_o1 * (tpOmega - 1.) |
---|
323 | + bMix_s1 * (tpSigma - 1.) |
---|
324 | + bMix_Omega * std::log(tmOmega) |
---|
325 | + bMix_sm * std::log(tmSigma) |
---|
326 | + bMix_oL * std::log(tpOmega) |
---|
327 | + bMix_sL * std::log(tpSigma); |
---|
328 | /* G4cout << "cross3 "<< cross<<" " |
---|
329 | <<bMix_o1<<" " |
---|
330 | <<bMix_s1<<" " |
---|
331 | <<bMix_Omega<<" " |
---|
332 | <<bMix_sm<<" " |
---|
333 | <<bMix_oL<<" " |
---|
334 | <<bMix_sL<<" " |
---|
335 | <<tpOmega<<" " |
---|
336 | <<tpSigma<<" " |
---|
337 | <<tmOmega<<" " |
---|
338 | <<tmSigma<<" " |
---|
339 | <<tpOmega<<" " |
---|
340 | <<tpSigma |
---|
341 | <<G4endl; |
---|
342 | */ |
---|
343 | return cross; |
---|
344 | |
---|
345 | } |
---|