| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4RDeBremsstrahlungSpectrum.cc,v 1.15 2006/06/29 19:41:58 gunter Exp $
>> 27 // GEANT4 tag $Name: geant4-09-01-ref-00 $
26 // 28 //
27 // ------------------------------------------- 29 // -------------------------------------------------------------------
28 // 30 //
29 // GEANT4 Class file 31 // GEANT4 Class file
30 // 32 //
31 // 33 //
32 // File name: G4RDeBremsstrahlungSpectrum 34 // File name: G4RDeBremsstrahlungSpectrum
33 // 35 //
34 // Author: V.Ivanchenko (Vladimir.Ivanc 36 // Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
35 // 37 //
36 // Creation date: 29 September 2001 38 // Creation date: 29 September 2001
37 // 39 //
38 // Modifications: 40 // Modifications:
39 // 10.10.01 MGP Revision to improve code qua 41 // 10.10.01 MGP Revision to improve code quality and consistency with design
40 // 15.11.01 VI Update spectrum model Bethe- 42 // 15.11.01 VI Update spectrum model Bethe-Haitler spectrum at high energy
41 // 30.05.02 VI Update interpolation between 43 // 30.05.02 VI Update interpolation between 2 last energy points in the
42 // parametrisation 44 // parametrisation
43 // 21.02.03 V.Ivanchenko Energy bins are d 45 // 21.02.03 V.Ivanchenko Energy bins are defined in the constructor
44 // 28.02.03 V.Ivanchenko Filename is defin 46 // 28.02.03 V.Ivanchenko Filename is defined in the constructor
45 // 47 //
46 // ------------------------------------------- 48 // -------------------------------------------------------------------
47 49
48 #include "G4RDeBremsstrahlungSpectrum.hh" 50 #include "G4RDeBremsstrahlungSpectrum.hh"
49 #include "G4RDBremsstrahlungParameters.hh" 51 #include "G4RDBremsstrahlungParameters.hh"
50 #include "Randomize.hh" 52 #include "Randomize.hh"
51 #include "G4SystemOfUnits.hh" <<
52 53
53 54
54 G4RDeBremsstrahlungSpectrum::G4RDeBremsstrahlu 55 G4RDeBremsstrahlungSpectrum::G4RDeBremsstrahlungSpectrum(const G4DataVector& bins,
55 const G4String& name):G4RDVEnergySpectrum(), 56 const G4String& name):G4RDVEnergySpectrum(),
56 lowestE(0.1*eV), 57 lowestE(0.1*eV),
57 xp(bins) 58 xp(bins)
58 { 59 {
59 length = xp.size(); 60 length = xp.size();
60 theBRparam = new G4RDBremsstrahlungParameter 61 theBRparam = new G4RDBremsstrahlungParameters(name,length+1);
61 62
62 verbose = 0; 63 verbose = 0;
63 } 64 }
64 65
65 66
66 G4RDeBremsstrahlungSpectrum::~G4RDeBremsstrahl 67 G4RDeBremsstrahlungSpectrum::~G4RDeBremsstrahlungSpectrum()
67 { 68 {
68 delete theBRparam; 69 delete theBRparam;
69 } 70 }
70 71
71 72
72 G4double G4RDeBremsstrahlungSpectrum::Probabil 73 G4double G4RDeBremsstrahlungSpectrum::Probability(G4int Z,
73 74 G4double tmin,
74 75 G4double tmax,
75 76 G4double e,
76 77 G4int,
77 const G 78 const G4ParticleDefinition*) const
78 { 79 {
79 G4double tm = std::min(tmax, e); 80 G4double tm = std::min(tmax, e);
80 G4double t0 = std::max(tmin, lowestE); 81 G4double t0 = std::max(tmin, lowestE);
81 if(t0 >= tm) return 0.0; 82 if(t0 >= tm) return 0.0;
82 83
83 t0 /= e; 84 t0 /= e;
84 tm /= e; 85 tm /= e;
85 86
86 G4double z0 = lowestE/e; 87 G4double z0 = lowestE/e;
87 G4DataVector p; 88 G4DataVector p;
88 89
89 // Access parameters 90 // Access parameters
90 for (size_t i=0; i<=length; i++) { 91 for (size_t i=0; i<=length; i++) {
91 p.push_back(theBRparam->Parameter(i, Z, e) 92 p.push_back(theBRparam->Parameter(i, Z, e));
92 } 93 }
93 94
94 G4double x = IntSpectrum(t0, tm, p); 95 G4double x = IntSpectrum(t0, tm, p);
95 G4double y = IntSpectrum(z0, 1.0, p); 96 G4double y = IntSpectrum(z0, 1.0, p);
96 97
97 98
98 if(1 < verbose) { 99 if(1 < verbose) {
99 G4cout << "tcut(MeV)= " << tmin/MeV 100 G4cout << "tcut(MeV)= " << tmin/MeV
100 << "; tMax(MeV)= " << tmax/MeV 101 << "; tMax(MeV)= " << tmax/MeV
101 << "; t0= " << t0 102 << "; t0= " << t0
102 << "; tm= " << tm 103 << "; tm= " << tm
103 << "; xp[0]= " << xp[0] 104 << "; xp[0]= " << xp[0]
104 << "; z= " << z0 105 << "; z= " << z0
105 << "; val= " << x 106 << "; val= " << x
106 << "; nor= " << y 107 << "; nor= " << y
107 << G4endl; 108 << G4endl;
108 } 109 }
109 p.clear(); 110 p.clear();
110 111
111 if(y > 0.0) x /= y; 112 if(y > 0.0) x /= y;
112 else x = 0.0; 113 else x = 0.0;
113 // if(x < 0.0) x = 0.0; 114 // if(x < 0.0) x = 0.0;
114 115
115 return x; 116 return x;
116 } 117 }
117 118
118 119
119 G4double G4RDeBremsstrahlungSpectrum::AverageE 120 G4double G4RDeBremsstrahlungSpectrum::AverageEnergy(G4int Z,
120 121 G4double tmin,
121 122 G4double tmax,
122 123 G4double e,
123 124 G4int,
124 const G4ParticleDefinition*) con 125 const G4ParticleDefinition*) const
125 { 126 {
126 G4double tm = std::min(tmax, e); 127 G4double tm = std::min(tmax, e);
127 G4double t0 = std::max(tmin, lowestE); 128 G4double t0 = std::max(tmin, lowestE);
128 if(t0 >= tm) return 0.0; 129 if(t0 >= tm) return 0.0;
129 130
130 t0 /= e; 131 t0 /= e;
131 tm /= e; 132 tm /= e;
132 133
133 G4double z0 = lowestE/e; 134 G4double z0 = lowestE/e;
134 135
135 G4DataVector p; 136 G4DataVector p;
136 137
137 // Access parameters 138 // Access parameters
138 for (size_t i=0; i<=length; i++) { 139 for (size_t i=0; i<=length; i++) {
139 p.push_back(theBRparam->Parameter(i, Z, 140 p.push_back(theBRparam->Parameter(i, Z, e));
140 } 141 }
141 142
142 G4double x = AverageValue(t0, tm, p); 143 G4double x = AverageValue(t0, tm, p);
143 G4double y = IntSpectrum(z0, 1.0, p); 144 G4double y = IntSpectrum(z0, 1.0, p);
144 145
145 // Add integrant over lowest energies 146 // Add integrant over lowest energies
146 G4double zmin = tmin/e; 147 G4double zmin = tmin/e;
147 if(zmin < t0) { 148 if(zmin < t0) {
148 G4double c = std::sqrt(theBRparam->Parame 149 G4double c = std::sqrt(theBRparam->ParameterC(Z));
149 x += p[0]*(t0 - zmin - c*(std::atan(t0/c) 150 x += p[0]*(t0 - zmin - c*(std::atan(t0/c) - std::atan(zmin/c)));
150 } 151 }
151 x *= e; 152 x *= e;
152 153
153 if(1 < verbose) { 154 if(1 < verbose) {
154 G4cout << "tcut(MeV)= " << tmin/MeV 155 G4cout << "tcut(MeV)= " << tmin/MeV
155 << "; tMax(MeV)= " << tmax/MeV 156 << "; tMax(MeV)= " << tmax/MeV
156 << "; e(MeV)= " << e/MeV 157 << "; e(MeV)= " << e/MeV
157 << "; t0= " << t0 158 << "; t0= " << t0
158 << "; tm= " << tm 159 << "; tm= " << tm
159 << "; y= " << y 160 << "; y= " << y
160 << "; x= " << x 161 << "; x= " << x
161 << G4endl; 162 << G4endl;
162 } 163 }
163 p.clear(); 164 p.clear();
164 165
165 if(y > 0.0) x /= y; 166 if(y > 0.0) x /= y;
166 else x = 0.0; 167 else x = 0.0;
167 // if(x < 0.0) x = 0.0; 168 // if(x < 0.0) x = 0.0;
168 169
169 return x; 170 return x;
170 } 171 }
171 172
172 173
173 G4double G4RDeBremsstrahlungSpectrum::SampleEn 174 G4double G4RDeBremsstrahlungSpectrum::SampleEnergy(G4int Z,
174 175 G4double tmin,
175 176 G4double tmax,
176 177 G4double e,
177 178 G4int,
178 const G4ParticleDefinition*) cons 179 const G4ParticleDefinition*) const
179 { 180 {
180 G4double tm = std::min(tmax, e); 181 G4double tm = std::min(tmax, e);
181 G4double t0 = std::max(tmin, lowestE); 182 G4double t0 = std::max(tmin, lowestE);
182 if(t0 >= tm) return 0.0; 183 if(t0 >= tm) return 0.0;
183 184
184 t0 /= e; 185 t0 /= e;
185 tm /= e; 186 tm /= e;
186 187
187 G4DataVector p; 188 G4DataVector p;
188 189
189 for (size_t i=0; i<=length; i++) { 190 for (size_t i=0; i<=length; i++) {
190 p.push_back(theBRparam->Parameter(i, Z, e) 191 p.push_back(theBRparam->Parameter(i, Z, e));
191 } 192 }
192 G4double amaj = std::max(p[length], 1. - (p[ 193 G4double amaj = std::max(p[length], 1. - (p[1] - p[0])*xp[0]/(xp[1] - xp[0]) );
193 194
194 G4double amax = std::log(tm); 195 G4double amax = std::log(tm);
195 G4double amin = std::log(t0); 196 G4double amin = std::log(t0);
196 G4double tgam, q, fun; 197 G4double tgam, q, fun;
197 198
198 do { 199 do {
199 G4double x = amin + G4UniformRand()*(amax 200 G4double x = amin + G4UniformRand()*(amax - amin);
200 tgam = std::exp(x); 201 tgam = std::exp(x);
201 fun = Function(tgam, p); 202 fun = Function(tgam, p);
202 203
203 if(fun > amaj) { 204 if(fun > amaj) {
204 G4cout << "WARNING in G4RDeBremsstra 205 G4cout << "WARNING in G4RDeBremsstrahlungSpectrum::SampleEnergy:"
205 << " Majoranta " << amaj 206 << " Majoranta " << amaj
206 << " < " << fun 207 << " < " << fun
207 << G4endl; 208 << G4endl;
208 } 209 }
209 210
210 q = amaj * G4UniformRand(); 211 q = amaj * G4UniformRand();
211 } while (q > fun); 212 } while (q > fun);
212 213
213 tgam *= e; 214 tgam *= e;
214 215
215 p.clear(); 216 p.clear();
216 217
217 return tgam; 218 return tgam;
218 } 219 }
219 220
220 G4double G4RDeBremsstrahlungSpectrum::IntSpect 221 G4double G4RDeBremsstrahlungSpectrum::IntSpectrum(G4double xMin,
221 222 G4double xMax,
222 const G4DataVector& p) const 223 const G4DataVector& p) const
223 { 224 {
224 G4double x1 = std::min(xMin, xp[0]); 225 G4double x1 = std::min(xMin, xp[0]);
225 G4double x2 = std::min(xMax, xp[0]); 226 G4double x2 = std::min(xMax, xp[0]);
226 G4double sum = 0.0; 227 G4double sum = 0.0;
227 228
228 if(x1 < x2) { 229 if(x1 < x2) {
229 G4double k = (p[1] - p[0])/(xp[1] - xp[0]) 230 G4double k = (p[1] - p[0])/(xp[1] - xp[0]);
230 sum += (1. - k*xp[0])*std::log(x2/x1) + k* 231 sum += (1. - k*xp[0])*std::log(x2/x1) + k*(x2 - x1);
231 } 232 }
232 233
233 for (size_t i=0; i<length-1; i++) { 234 for (size_t i=0; i<length-1; i++) {
234 x1 = std::max(xMin, xp[i]); 235 x1 = std::max(xMin, xp[i]);
235 x2 = std::min(xMax, xp[i+1]); 236 x2 = std::min(xMax, xp[i+1]);
236 if(x1 < x2) { 237 if(x1 < x2) {
237 G4double z1 = p[i]; 238 G4double z1 = p[i];
238 G4double z2 = p[i+1]; 239 G4double z2 = p[i+1];
239 sum += z2 - z1 + std::log(x2/x1)*(z1*x2 240 sum += z2 - z1 + std::log(x2/x1)*(z1*x2 - z2*x1)/(x2 - x1);
240 } 241 }
241 } 242 }
242 if(sum < 0.0) sum = 0.0; 243 if(sum < 0.0) sum = 0.0;
243 return sum; 244 return sum;
244 } 245 }
245 246
246 G4double G4RDeBremsstrahlungSpectrum::AverageV 247 G4double G4RDeBremsstrahlungSpectrum::AverageValue(G4double xMin,
247 248 G4double xMax,
248 const G4DataVector& p) const 249 const G4DataVector& p) const
249 { 250 {
250 G4double x1 = std::min(xMin, xp[0]); 251 G4double x1 = std::min(xMin, xp[0]);
251 G4double x2 = std::min(xMax, xp[0]); 252 G4double x2 = std::min(xMax, xp[0]);
252 G4double z1 = x1; 253 G4double z1 = x1;
253 G4double z2 = x2; 254 G4double z2 = x2;
254 G4double sum = 0.0; 255 G4double sum = 0.0;
255 256
256 if(x1 < x2) { 257 if(x1 < x2) {
257 G4double k = (p[1] - p[0])/(xp[1] - xp[0]) 258 G4double k = (p[1] - p[0])/(xp[1] - xp[0]);
258 sum += (z2 - z1)*(1. - k*xp[0]); 259 sum += (z2 - z1)*(1. - k*xp[0]);
259 z1 *= x1; 260 z1 *= x1;
260 z2 *= x2; 261 z2 *= x2;
261 sum += 0.5*k*(z2 - z1); 262 sum += 0.5*k*(z2 - z1);
262 } 263 }
263 264
264 for (size_t i=0; i<length-1; i++) { 265 for (size_t i=0; i<length-1; i++) {
265 x1 = std::max(xMin, xp[i]); 266 x1 = std::max(xMin, xp[i]);
266 x2 = std::min(xMax, xp[i+1]); 267 x2 = std::min(xMax, xp[i+1]);
267 if(x1 < x2) { 268 if(x1 < x2) {
268 z1 = p[i]; 269 z1 = p[i];
269 z2 = p[i+1]; 270 z2 = p[i+1];
270 sum += 0.5*(z2 - z1)*(x2 + x1) + z1*x2 - 271 sum += 0.5*(z2 - z1)*(x2 + x1) + z1*x2 - z2*x1;
271 } 272 }
272 } 273 }
273 if(sum < 0.0) sum = 0.0; 274 if(sum < 0.0) sum = 0.0;
274 return sum; 275 return sum;
275 } 276 }
276 277
277 G4double G4RDeBremsstrahlungSpectrum::Function 278 G4double G4RDeBremsstrahlungSpectrum::Function(G4double x,
278 const G4DataVector& p) const 279 const G4DataVector& p) const
279 { 280 {
280 G4double f = 0.0; 281 G4double f = 0.0;
281 282
282 if(x <= xp[0]) { 283 if(x <= xp[0]) {
283 f = p[0] + (p[1] - p[0])*(x - xp[0])/(xp[1 284 f = p[0] + (p[1] - p[0])*(x - xp[0])/(xp[1] - xp[0]);
284 285
285 } else { 286 } else {
286 287
287 for (size_t i=0; i<length-1; i++) { 288 for (size_t i=0; i<length-1; i++) {
288 289
289 if(x <= xp[i+1]) { 290 if(x <= xp[i+1]) {
290 f = p[i] + (p[i+1] - p[i])*(x - xp[i]) 291 f = p[i] + (p[i+1] - p[i])*(x - xp[i])/(xp[i+1] - xp[i]);
291 break; 292 break;
292 } 293 }
293 } 294 }
294 } 295 }
295 return f; 296 return f;
296 } 297 }
297 298
298 void G4RDeBremsstrahlungSpectrum::PrintData() 299 void G4RDeBremsstrahlungSpectrum::PrintData() const
299 { theBRparam->PrintData(); } 300 { theBRparam->PrintData(); }
300 301
301 G4double G4RDeBremsstrahlungSpectrum::Excitati 302 G4double G4RDeBremsstrahlungSpectrum::Excitation(G4int , G4double ) const
302 { 303 {
303 return 0.0; 304 return 0.0;
304 } 305 }
305 306
306 G4double G4RDeBremsstrahlungSpectrum::MaxEnerg 307 G4double G4RDeBremsstrahlungSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
307 G4int, // Z, 308 G4int, // Z,
308 const G4ParticleDefinition*) 309 const G4ParticleDefinition*) const
309 { 310 {
310 return kineticEnergy; 311 return kineticEnergy;
311 } 312 }
312 313