| 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 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4RDeBremsstrahlungSpectrum
33 //
34 // Author: V.Ivanchenko (Vladimir.Ivanc
35 //
36 // Creation date: 29 September 2001
37 //
38 // Modifications:
39 // 10.10.01 MGP Revision to improve code qua
40 // 15.11.01 VI Update spectrum model Bethe-
41 // 30.05.02 VI Update interpolation between
42 // parametrisation
43 // 21.02.03 V.Ivanchenko Energy bins are d
44 // 28.02.03 V.Ivanchenko Filename is defin
45 //
46 // -------------------------------------------
47
48 #include "G4RDeBremsstrahlungSpectrum.hh"
49 #include "G4RDBremsstrahlungParameters.hh"
50 #include "Randomize.hh"
51 #include "G4SystemOfUnits.hh"
52
53
54 G4RDeBremsstrahlungSpectrum::G4RDeBremsstrahlu
55 const G4String& name):G4RDVEnergySpectrum(),
56 lowestE(0.1*eV),
57 xp(bins)
58 {
59 length = xp.size();
60 theBRparam = new G4RDBremsstrahlungParameter
61
62 verbose = 0;
63 }
64
65
66 G4RDeBremsstrahlungSpectrum::~G4RDeBremsstrahl
67 {
68 delete theBRparam;
69 }
70
71
72 G4double G4RDeBremsstrahlungSpectrum::Probabil
73
74
75
76
77 const G
78 {
79 G4double tm = std::min(tmax, e);
80 G4double t0 = std::max(tmin, lowestE);
81 if(t0 >= tm) return 0.0;
82
83 t0 /= e;
84 tm /= e;
85
86 G4double z0 = lowestE/e;
87 G4DataVector p;
88
89 // Access parameters
90 for (size_t i=0; i<=length; i++) {
91 p.push_back(theBRparam->Parameter(i, Z, e)
92 }
93
94 G4double x = IntSpectrum(t0, tm, p);
95 G4double y = IntSpectrum(z0, 1.0, p);
96
97
98 if(1 < verbose) {
99 G4cout << "tcut(MeV)= " << tmin/MeV
100 << "; tMax(MeV)= " << tmax/MeV
101 << "; t0= " << t0
102 << "; tm= " << tm
103 << "; xp[0]= " << xp[0]
104 << "; z= " << z0
105 << "; val= " << x
106 << "; nor= " << y
107 << G4endl;
108 }
109 p.clear();
110
111 if(y > 0.0) x /= y;
112 else x = 0.0;
113 // if(x < 0.0) x = 0.0;
114
115 return x;
116 }
117
118
119 G4double G4RDeBremsstrahlungSpectrum::AverageE
120
121
122
123
124 const G4ParticleDefinition*) con
125 {
126 G4double tm = std::min(tmax, e);
127 G4double t0 = std::max(tmin, lowestE);
128 if(t0 >= tm) return 0.0;
129
130 t0 /= e;
131 tm /= e;
132
133 G4double z0 = lowestE/e;
134
135 G4DataVector p;
136
137 // Access parameters
138 for (size_t i=0; i<=length; i++) {
139 p.push_back(theBRparam->Parameter(i, Z,
140 }
141
142 G4double x = AverageValue(t0, tm, p);
143 G4double y = IntSpectrum(z0, 1.0, p);
144
145 // Add integrant over lowest energies
146 G4double zmin = tmin/e;
147 if(zmin < t0) {
148 G4double c = std::sqrt(theBRparam->Parame
149 x += p[0]*(t0 - zmin - c*(std::atan(t0/c)
150 }
151 x *= e;
152
153 if(1 < verbose) {
154 G4cout << "tcut(MeV)= " << tmin/MeV
155 << "; tMax(MeV)= " << tmax/MeV
156 << "; e(MeV)= " << e/MeV
157 << "; t0= " << t0
158 << "; tm= " << tm
159 << "; y= " << y
160 << "; x= " << x
161 << G4endl;
162 }
163 p.clear();
164
165 if(y > 0.0) x /= y;
166 else x = 0.0;
167 // if(x < 0.0) x = 0.0;
168
169 return x;
170 }
171
172
173 G4double G4RDeBremsstrahlungSpectrum::SampleEn
174
175
176
177
178 const G4ParticleDefinition*) cons
179 {
180 G4double tm = std::min(tmax, e);
181 G4double t0 = std::max(tmin, lowestE);
182 if(t0 >= tm) return 0.0;
183
184 t0 /= e;
185 tm /= e;
186
187 G4DataVector p;
188
189 for (size_t i=0; i<=length; i++) {
190 p.push_back(theBRparam->Parameter(i, Z, e)
191 }
192 G4double amaj = std::max(p[length], 1. - (p[
193
194 G4double amax = std::log(tm);
195 G4double amin = std::log(t0);
196 G4double tgam, q, fun;
197
198 do {
199 G4double x = amin + G4UniformRand()*(amax
200 tgam = std::exp(x);
201 fun = Function(tgam, p);
202
203 if(fun > amaj) {
204 G4cout << "WARNING in G4RDeBremsstra
205 << " Majoranta " << amaj
206 << " < " << fun
207 << G4endl;
208 }
209
210 q = amaj * G4UniformRand();
211 } while (q > fun);
212
213 tgam *= e;
214
215 p.clear();
216
217 return tgam;
218 }
219
220 G4double G4RDeBremsstrahlungSpectrum::IntSpect
221
222 const G4DataVector& p) const
223 {
224 G4double x1 = std::min(xMin, xp[0]);
225 G4double x2 = std::min(xMax, xp[0]);
226 G4double sum = 0.0;
227
228 if(x1 < x2) {
229 G4double k = (p[1] - p[0])/(xp[1] - xp[0])
230 sum += (1. - k*xp[0])*std::log(x2/x1) + k*
231 }
232
233 for (size_t i=0; i<length-1; i++) {
234 x1 = std::max(xMin, xp[i]);
235 x2 = std::min(xMax, xp[i+1]);
236 if(x1 < x2) {
237 G4double z1 = p[i];
238 G4double z2 = p[i+1];
239 sum += z2 - z1 + std::log(x2/x1)*(z1*x2
240 }
241 }
242 if(sum < 0.0) sum = 0.0;
243 return sum;
244 }
245
246 G4double G4RDeBremsstrahlungSpectrum::AverageV
247
248 const G4DataVector& p) const
249 {
250 G4double x1 = std::min(xMin, xp[0]);
251 G4double x2 = std::min(xMax, xp[0]);
252 G4double z1 = x1;
253 G4double z2 = x2;
254 G4double sum = 0.0;
255
256 if(x1 < x2) {
257 G4double k = (p[1] - p[0])/(xp[1] - xp[0])
258 sum += (z2 - z1)*(1. - k*xp[0]);
259 z1 *= x1;
260 z2 *= x2;
261 sum += 0.5*k*(z2 - z1);
262 }
263
264 for (size_t i=0; i<length-1; i++) {
265 x1 = std::max(xMin, xp[i]);
266 x2 = std::min(xMax, xp[i+1]);
267 if(x1 < x2) {
268 z1 = p[i];
269 z2 = p[i+1];
270 sum += 0.5*(z2 - z1)*(x2 + x1) + z1*x2 -
271 }
272 }
273 if(sum < 0.0) sum = 0.0;
274 return sum;
275 }
276
277 G4double G4RDeBremsstrahlungSpectrum::Function
278 const G4DataVector& p) const
279 {
280 G4double f = 0.0;
281
282 if(x <= xp[0]) {
283 f = p[0] + (p[1] - p[0])*(x - xp[0])/(xp[1
284
285 } else {
286
287 for (size_t i=0; i<length-1; i++) {
288
289 if(x <= xp[i+1]) {
290 f = p[i] + (p[i+1] - p[i])*(x - xp[i])
291 break;
292 }
293 }
294 }
295 return f;
296 }
297
298 void G4RDeBremsstrahlungSpectrum::PrintData()
299 { theBRparam->PrintData(); }
300
301 G4double G4RDeBremsstrahlungSpectrum::Excitati
302 {
303 return 0.0;
304 }
305
306 G4double G4RDeBremsstrahlungSpectrum::MaxEnerg
307 G4int, // Z,
308 const G4ParticleDefinition*)
309 {
310 return kineticEnergy;
311 }
312