| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
>> 23 // $Id: G4eIonisationSpectrum.cc,v 1.22 2003/06/16 17:00:35 gunter Exp $
>> 24 // GEANT4 tag $Name: geant4-05-02 $
26 // 25 //
27 // ------------------------------------------- 26 // -------------------------------------------------------------------
28 // 27 //
29 // GEANT4 Class file 28 // GEANT4 Class file
30 // 29 //
31 // 30 //
32 // File name: G4eIonisationSpectrum 31 // File name: G4eIonisationSpectrum
33 // 32 //
34 // Author: V.Ivanchenko (Vladimir.Ivanc 33 // Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
35 // << 34 //
36 // Creation date: 29 September 2001 35 // Creation date: 29 September 2001
37 // 36 //
38 // Modifications: << 37 // Modifications:
39 // 10.10.2001 MGP Revision to improv << 38 // 10.10.2001 MGP Revision to improve code quality and
40 // consistency with d 39 // consistency with design
41 // 02.11.2001 VI Optimize sampling << 40 // 02.11.2001 VI Optimize sampling of energy
42 // 29.11.2001 VI New parametrisatio << 41 // 29.11.2001 VI New parametrisation
43 // 19.04.2002 VI Add protection in << 42 // 19.04.2002 VI Add protection in case of energy below binding
44 // 30.05.2002 VI Update to 24-param 43 // 30.05.2002 VI Update to 24-parameters data
45 // 11.07.2002 VI Fix in integration 44 // 11.07.2002 VI Fix in integration over spectrum
46 // 23.03.2009 LP Added protection a <<
47 // faulty database fi <<
48 // 45 //
49 // ------------------------------------------- 46 // -------------------------------------------------------------------
50 // 47 //
51 48
52 #include "G4eIonisationSpectrum.hh" 49 #include "G4eIonisationSpectrum.hh"
53 #include "G4AtomicTransitionManager.hh" 50 #include "G4AtomicTransitionManager.hh"
54 #include "G4AtomicShell.hh" 51 #include "G4AtomicShell.hh"
55 #include "G4DataVector.hh" 52 #include "G4DataVector.hh"
56 #include "Randomize.hh" 53 #include "Randomize.hh"
57 #include "G4PhysicalConstants.hh" << 54
58 #include "G4SystemOfUnits.hh" <<
59 #include "G4Exp.hh" <<
60 55
61 G4eIonisationSpectrum::G4eIonisationSpectrum() 56 G4eIonisationSpectrum::G4eIonisationSpectrum():G4VEnergySpectrum(),
62 lowestE(0.1*eV), 57 lowestE(0.1*eV),
63 factor(1.3), 58 factor(1.3),
64 iMax(24), 59 iMax(24),
65 verbose(0) 60 verbose(0)
66 { 61 {
67 theParam = new G4eIonisationParameters(); 62 theParam = new G4eIonisationParameters();
68 } 63 }
69 64
70 65
71 G4eIonisationSpectrum::~G4eIonisationSpectrum( << 66 G4eIonisationSpectrum::~G4eIonisationSpectrum()
72 { 67 {
73 delete theParam; 68 delete theParam;
74 } 69 }
75 70
76 71
77 G4double G4eIonisationSpectrum::Probability(G4 << 72 G4double G4eIonisationSpectrum::Probability(G4int Z,
78 G4double tMin, << 73 G4double tMin,
79 G4double tMax, << 74 G4double tMax,
80 G4double e, 75 G4double e,
81 G4int shell, 76 G4int shell,
82 const G4ParticleDefinition* ) co 77 const G4ParticleDefinition* ) const
83 { 78 {
84 // Please comment what Probability does and << 79 // Please comment what Probability does and what are the three
85 // functions mentioned below 80 // functions mentioned below
86 // Describe the algorithms used 81 // Describe the algorithms used
87 82
88 G4double eMax = MaxEnergyOfSecondaries(e); 83 G4double eMax = MaxEnergyOfSecondaries(e);
89 G4double t0 = std::max(tMin, lowestE); 84 G4double t0 = std::max(tMin, lowestE);
90 G4double tm = std::min(tMax, eMax); 85 G4double tm = std::min(tMax, eMax);
91 if(t0 >= tm) return 0.0; 86 if(t0 >= tm) return 0.0;
92 87
93 G4double bindingEnergy = (G4AtomicTransition 88 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
94 Shell(Z, shell)->Bi 89 Shell(Z, shell)->BindingEnergy();
95 90
96 if(e <= bindingEnergy) return 0.0; 91 if(e <= bindingEnergy) return 0.0;
97 92
98 G4double energy = e + bindingEnergy; 93 G4double energy = e + bindingEnergy;
99 94
100 G4double x1 = std::min(0.5,(t0 + bindingEner 95 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
101 G4double x2 = std::min(0.5,(tm + bindingEner 96 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
102 97
103 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0. 98 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
104 G4cout << "G4eIonisationSpectrum::Probabil 99 G4cout << "G4eIonisationSpectrum::Probability: Z= " << Z
105 << "; shell= " << shell 100 << "; shell= " << shell
106 << "; E(keV)= " << e/keV 101 << "; E(keV)= " << e/keV
107 << "; Eb(keV)= " << bindingEnergy/k 102 << "; Eb(keV)= " << bindingEnergy/keV
108 << "; x1= " << x1 << 103 << "; x1= " << x1
109 << "; x2= " << x2 << 104 << "; x2= " << x2
110 << G4endl; 105 << G4endl;
111 << 106
112 } 107 }
113 108
114 G4DataVector p; 109 G4DataVector p;
115 110
116 // Access parameters 111 // Access parameters
117 for (G4int i=0; i<iMax; i++) << 112 for (G4int i=0; i<iMax; i++)
118 { 113 {
119 G4double x = theParam->Parameter(Z, shell, 114 G4double x = theParam->Parameter(Z, shell, i, e);
120 if(i<4) x /= energy; << 115 if(i<4) x /= energy;
121 p.push_back(x); << 116 p.push_back(x);
122 } 117 }
123 118
124 if(p[3] > 0.5) p[3] = 0.5; 119 if(p[3] > 0.5) p[3] = 0.5;
125 << 120
126 G4double gLocal = energy/electron_mass_c2 + << 121 G4double g = energy/electron_mass_c2 + 1.;
127 p.push_back((2.0*gLocal - 1.0)/(gLocal*gLoca << 122 p.push_back((2.0*g - 1.0)/(g*g));
128 << 123
129 //Add protection against division by zero: a << 124 p[iMax-1] = Function(p[3], p);
130 //parameter p[3] appears in the denominator. <<
131 if (p[3] > 0) <<
132 p[iMax-1] = Function(p[3], p); <<
133 else <<
134 { <<
135 G4cout << "WARNING: G4eIonisationSpectru <<
136 << "parameter p[3] <= 0. G4LEDATA dabat <<
137 << Z << ". Please check and/or update i <<
138 } <<
139 125
140 if(e >= 1. && e <= 0. && Z == 4) p.push_back 126 if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
141 127
142 << 128
143 G4double val = IntSpectrum(x1, x2, p); 129 G4double val = IntSpectrum(x1, x2, p);
144 G4double x0 = (lowestE + bindingEnergy)/ene 130 G4double x0 = (lowestE + bindingEnergy)/energy;
145 G4double nor = IntSpectrum(x0, 0.5, p); 131 G4double nor = IntSpectrum(x0, 0.5, p);
146 << 132
147 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0. 133 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
148 G4cout << "tcut= " << tMin << 134 G4cout << "tcut= " << tMin
149 << "; tMax= " << tMax << 135 << "; tMax= " << tMax
150 << "; x0= " << x0 << 136 << "; x0= " << x0
151 << "; x1= " << x1 << 137 << "; x1= " << x1
152 << "; x2= " << x2 << 138 << "; x2= " << x2
153 << "; val= " << val << 139 << "; val= " << val
154 << "; nor= " << nor << 140 << "; nor= " << nor
155 << "; sum= " << p[0] << 141 << "; sum= " << p[0]
156 << "; a= " << p[1] << 142 << "; a= " << p[1]
157 << "; b= " << p[2] << 143 << "; b= " << p[2]
158 << "; c= " << p[3] << 144 << "; c= " << p[3]
159 << G4endl; 145 << G4endl;
160 if(shell == 1) G4cout << "============" << << 146 if(shell == 1) G4cout << "============" << G4endl;
161 } 147 }
162 148
163 p.clear(); 149 p.clear();
164 150
165 if(nor > 0.0) val /= nor; 151 if(nor > 0.0) val /= nor;
166 else val = 0.0; 152 else val = 0.0;
167 153
168 return val; << 154 return val;
169 } 155 }
170 156
171 157
172 G4double G4eIonisationSpectrum::AverageEnergy( 158 G4double G4eIonisationSpectrum::AverageEnergy(G4int Z,
173 G4double tMin, << 159 G4double tMin,
174 G4double tMax, << 160 G4double tMax,
175 G4double e, 161 G4double e,
176 G4int shell, 162 G4int shell,
177 const G4ParticleDefinition* ) 163 const G4ParticleDefinition* ) const
178 { 164 {
179 // Please comment what AverageEnergy does an << 165 // Please comment what AverageEnergy does and what are the three
180 // functions mentioned below 166 // functions mentioned below
181 // Describe the algorithms used 167 // Describe the algorithms used
182 168
183 G4double eMax = MaxEnergyOfSecondaries(e); 169 G4double eMax = MaxEnergyOfSecondaries(e);
184 G4double t0 = std::max(tMin, lowestE); 170 G4double t0 = std::max(tMin, lowestE);
185 G4double tm = std::min(tMax, eMax); 171 G4double tm = std::min(tMax, eMax);
186 if(t0 >= tm) return 0.0; 172 if(t0 >= tm) return 0.0;
187 173
188 G4double bindingEnergy = (G4AtomicTransition 174 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
189 Shell(Z, shell)->Bi 175 Shell(Z, shell)->BindingEnergy();
190 176
191 if(e <= bindingEnergy) return 0.0; 177 if(e <= bindingEnergy) return 0.0;
192 178
193 G4double energy = e + bindingEnergy; 179 G4double energy = e + bindingEnergy;
194 180
195 G4double x1 = std::min(0.5,(t0 + bindingEner 181 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
196 G4double x2 = std::min(0.5,(tm + bindingEner 182 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
197 183
198 if(verbose > 1) { 184 if(verbose > 1) {
199 G4cout << "G4eIonisationSpectrum::AverageE 185 G4cout << "G4eIonisationSpectrum::AverageEnergy: Z= " << Z
200 << "; shell= " << shell 186 << "; shell= " << shell
201 << "; E(keV)= " << e/keV 187 << "; E(keV)= " << e/keV
202 << "; bindingE(keV)= " << bindingEn 188 << "; bindingE(keV)= " << bindingEnergy/keV
203 << "; x1= " << x1 << 189 << "; x1= " << x1
204 << "; x2= " << x2 << 190 << "; x2= " << x2
205 << G4endl; 191 << G4endl;
206 } 192 }
207 193
208 G4DataVector p; 194 G4DataVector p;
209 195
210 // Access parameters 196 // Access parameters
211 for (G4int i=0; i<iMax; i++) << 197 for (G4int i=0; i<iMax; i++)
212 { 198 {
213 G4double x = theParam->Parameter(Z, shell, 199 G4double x = theParam->Parameter(Z, shell, i, e);
214 if(i<4) x /= energy; << 200 if(i<4) x /= energy;
215 p.push_back(x); 201 p.push_back(x);
216 } 202 }
217 203
218 if(p[3] > 0.5) p[3] = 0.5; 204 if(p[3] > 0.5) p[3] = 0.5;
219 205
220 G4double gLocal2 = energy/electron_mass_c2 + << 206 G4double g = energy/electron_mass_c2 + 1.;
221 p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLo << 207 p.push_back((2.0*g - 1.0)/(g*g));
222 <<
223 <<
224 //Add protection against division by zero: a <<
225 //parameter p[3] appears in the denominator. <<
226 if (p[3] > 0) <<
227 p[iMax-1] = Function(p[3], p); <<
228 else <<
229 { <<
230 G4cout << "WARNING: G4eIonisationSpectru <<
231 << "parameter p[3] <= 0. G4LEDATA dabat <<
232 << Z << ". Please check and/or update i <<
233 } <<
234 208
>> 209 p[iMax-1] = Function(p[3], p);
>> 210
235 G4double val = AverageValue(x1, x2, p); 211 G4double val = AverageValue(x1, x2, p);
236 G4double x0 = (lowestE + bindingEnergy)/ene 212 G4double x0 = (lowestE + bindingEnergy)/energy;
237 G4double nor = IntSpectrum(x0, 0.5, p); 213 G4double nor = IntSpectrum(x0, 0.5, p);
238 val *= energy; 214 val *= energy;
239 215
240 if(verbose > 1) { 216 if(verbose > 1) {
241 G4cout << "tcut(MeV)= " << tMin/MeV << 217 G4cout << "tcut(MeV)= " << tMin/MeV
242 << "; tMax(MeV)= " << tMax/MeV << 218 << "; tMax(MeV)= " << tMax/MeV
243 << "; x0= " << x0 << 219 << "; x0= " << x0
244 << "; x1= " << x1 << 220 << "; x1= " << x1
245 << "; x2= " << x2 << 221 << "; x2= " << x2
246 << "; val= " << val << 222 << "; val= " << val
247 << "; nor= " << nor << 223 << "; nor= " << nor
248 << "; sum= " << p[0] << 224 << "; sum= " << p[0]
249 << "; a= " << p[1] << 225 << "; a= " << p[1]
250 << "; b= " << p[2] << 226 << "; b= " << p[2]
251 << "; c= " << p[3] << 227 << "; c= " << p[3]
252 << G4endl; 228 << G4endl;
253 } 229 }
254 230
255 p.clear(); 231 p.clear();
256 232
257 if(nor > 0.0) val /= nor; 233 if(nor > 0.0) val /= nor;
258 else val = 0.0; 234 else val = 0.0;
259 235
260 return val; << 236 return val;
261 } 237 }
262 238
263 239
264 G4double G4eIonisationSpectrum::SampleEnergy(G 240 G4double G4eIonisationSpectrum::SampleEnergy(G4int Z,
265 G4double tMin, << 241 G4double tMin,
266 G4double tMax, << 242 G4double tMax,
267 G4double e, 243 G4double e,
268 G4int shell, 244 G4int shell,
269 const G4ParticleDefinition* ) c 245 const G4ParticleDefinition* ) const
270 { 246 {
271 // Please comment what SampleEnergy does 247 // Please comment what SampleEnergy does
272 G4double tDelta = 0.0; 248 G4double tDelta = 0.0;
273 G4double t0 = std::max(tMin, lowestE); 249 G4double t0 = std::max(tMin, lowestE);
274 G4double tm = std::min(tMax, MaxEnergyOfSeco 250 G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
275 if(t0 > tm) return tDelta; 251 if(t0 > tm) return tDelta;
276 252
277 G4double bindingEnergy = (G4AtomicTransition 253 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
278 Shell(Z, shell)->Bi 254 Shell(Z, shell)->BindingEnergy();
279 255
280 if(e <= bindingEnergy) return 0.0; 256 if(e <= bindingEnergy) return 0.0;
281 257
282 G4double energy = e + bindingEnergy; 258 G4double energy = e + bindingEnergy;
283 259
284 G4double x1 = std::min(0.5,(t0 + bindingEner 260 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
285 G4double x2 = std::min(0.5,(tm + bindingEner 261 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
286 if(x1 >= x2) return tDelta; 262 if(x1 >= x2) return tDelta;
287 263
288 if(verbose > 1) { 264 if(verbose > 1) {
289 G4cout << "G4eIonisationSpectrum::SampleEn 265 G4cout << "G4eIonisationSpectrum::SampleEnergy: Z= " << Z
290 << "; shell= " << shell 266 << "; shell= " << shell
291 << "; E(keV)= " << e/keV 267 << "; E(keV)= " << e/keV
292 << G4endl; 268 << G4endl;
293 } 269 }
294 270
295 // Access parameters 271 // Access parameters
296 G4DataVector p; 272 G4DataVector p;
297 273
298 // Access parameters 274 // Access parameters
299 for (G4int i=0; i<iMax; i++) << 275 for (G4int i=0; i<iMax; i++)
300 { 276 {
301 G4double x = theParam->Parameter(Z, shell, 277 G4double x = theParam->Parameter(Z, shell, i, e);
302 if(i<4) x /= energy; << 278 if(i<4) x /= energy;
303 p.push_back(x); 279 p.push_back(x);
304 } 280 }
305 281
306 if(p[3] > 0.5) p[3] = 0.5; 282 if(p[3] > 0.5) p[3] = 0.5;
307 283
308 G4double gLocal3 = energy/electron_mass_c2 + << 284 G4double g = energy/electron_mass_c2 + 1.;
309 p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLo << 285 p.push_back((2.0*g - 1.0)/(g*g));
310 <<
311 286
312 //Add protection against division by zero: a << 287 p[iMax-1] = Function(p[3], p);
313 //parameter p[3] appears in the denominator. <<
314 if (p[3] > 0) <<
315 p[iMax-1] = Function(p[3], p); <<
316 else <<
317 { <<
318 G4cout << "WARNING: G4eIonisationSpectru <<
319 << "parameter p[3] <= 0. G4LEDATA dabat <<
320 << Z << ". Please check and/or update i <<
321 } <<
322 288
323 G4double aria1 = 0.0; 289 G4double aria1 = 0.0;
324 G4double a1 = std::max(x1,p[1]); 290 G4double a1 = std::max(x1,p[1]);
325 G4double a2 = std::min(x2,p[3]); 291 G4double a2 = std::min(x2,p[3]);
326 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p); 292 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
327 G4double aria2 = 0.0; 293 G4double aria2 = 0.0;
328 G4double a3 = std::max(x1,p[3]); 294 G4double a3 = std::max(x1,p[3]);
329 G4double a4 = x2; 295 G4double a4 = x2;
330 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p); 296 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
331 297
332 G4double aria = (aria1 + aria2)*G4UniformRan 298 G4double aria = (aria1 + aria2)*G4UniformRand();
333 G4double amaj, fun, q, x, z1, z2, dx, dx1; 299 G4double amaj, fun, q, x, z1, z2, dx, dx1;
334 300
335 //======= First aria to sample ===== 301 //======= First aria to sample =====
336 302
337 if(aria <= aria1) { << 303 if(aria <= aria1) {
338 304
339 amaj = p[4]; 305 amaj = p[4];
340 for (G4int j=5; j<iMax; j++) { 306 for (G4int j=5; j<iMax; j++) {
341 if(p[j] > amaj) amaj = p[j]; 307 if(p[j] > amaj) amaj = p[j];
342 } 308 }
343 309
344 a1 = 1./a1; 310 a1 = 1./a1;
345 a2 = 1./a2; 311 a2 = 1./a2;
346 312
347 G4int i; 313 G4int i;
348 do { 314 do {
349 315
350 x = 1./(a2 + G4UniformRand()*(a1 - a2)); 316 x = 1./(a2 + G4UniformRand()*(a1 - a2));
351 z1 = p[1]; 317 z1 = p[1];
352 z2 = p[3]; 318 z2 = p[3];
353 dx = (p[2] - p[1]) / 3.0; 319 dx = (p[2] - p[1]) / 3.0;
354 dx1= G4Exp(std::log(p[3]/p[2]) / 16.0); << 320 dx1= exp(log(p[3]/p[2]) / 16.0);
355 for (i=4; i<iMax-1; i++) { 321 for (i=4; i<iMax-1; i++) {
356 322
357 if (i < 7) { 323 if (i < 7) {
358 z2 = z1 + dx; 324 z2 = z1 + dx;
359 } else if(iMax-2 == i) { 325 } else if(iMax-2 == i) {
360 z2 = p[3]; 326 z2 = p[3];
361 break; 327 break;
362 } else { 328 } else {
363 z2 = z1*dx1; 329 z2 = z1*dx1;
364 } 330 }
365 if(x >= z1 && x <= z2) break; 331 if(x >= z1 && x <= z2) break;
366 z1 = z2; 332 z1 = z2;
367 } 333 }
368 fun = p[i] + (x - z1) * (p[i+1] - p[i])/ 334 fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
369 335
370 if(fun > amaj) { 336 if(fun > amaj) {
371 G4cout << "WARNING in G4eIonisationS << 337 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
372 << " Majoranta " << amaj << 338 << " Majoranta " << amaj
373 << " < " << fun 339 << " < " << fun
374 << " in the first aria at x= 340 << " in the first aria at x= " << x
375 << G4endl; 341 << G4endl;
376 } 342 }
377 343
378 q = amaj*G4UniformRand(); 344 q = amaj*G4UniformRand();
379 345
380 } while (q >= fun); 346 } while (q >= fun);
381 347
382 //======= Second aria to sample ===== 348 //======= Second aria to sample =====
383 349
384 } else { 350 } else {
385 351
386 amaj = std::max(p[iMax-1], Function(0.5, p 352 amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
387 a1 = 1./a3; 353 a1 = 1./a3;
388 a2 = 1./a4; 354 a2 = 1./a4;
389 355
390 do { 356 do {
391 357
392 x = 1./(a2 + G4UniformRand()*(a1 - a2)); 358 x = 1./(a2 + G4UniformRand()*(a1 - a2));
393 fun = Function(x, p); 359 fun = Function(x, p);
394 360
395 if(fun > amaj) { 361 if(fun > amaj) {
396 G4cout << "WARNING in G4eIonisationS << 362 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
397 << " Majoranta " << amaj << 363 << " Majoranta " << amaj
398 << " < " << fun 364 << " < " << fun
399 << " in the second aria at x= 365 << " in the second aria at x= " << x
400 << G4endl; 366 << G4endl;
401 } 367 }
402 368
403 q = amaj*G4UniformRand(); 369 q = amaj*G4UniformRand();
404 370
405 } while (q >= fun); 371 } while (q >= fun);
406 372
407 } 373 }
408 374
409 p.clear(); 375 p.clear();
410 376
411 tDelta = x*energy - bindingEnergy; 377 tDelta = x*energy - bindingEnergy;
412 378
413 if(verbose > 1) { 379 if(verbose > 1) {
414 G4cout << "tcut(MeV)= " << tMin/MeV << 380 G4cout << "tcut(MeV)= " << tMin/MeV
415 << "; tMax(MeV)= " << tMax/MeV << 381 << "; tMax(MeV)= " << tMax/MeV
416 << "; x1= " << x1 << 382 << "; x1= " << x1
417 << "; x2= " << x2 << 383 << "; x2= " << x2
418 << "; a1= " << a1 << 384 << "; a1= " << a1
419 << "; a2= " << a2 << 385 << "; a2= " << a2
420 << "; x= " << x << 386 << "; x= " << x
421 << "; be= " << bindingEnergy << 387 << "; be= " << bindingEnergy
422 << "; e= " << e << 388 << "; e= " << e
423 << "; tDelta= " << tDelta << 389 << "; tDelta= " << tDelta
424 << G4endl; 390 << G4endl;
425 } 391 }
426 392
427 393
428 return tDelta; << 394 return tDelta;
429 } 395 }
430 396
431 397
432 G4double G4eIonisationSpectrum::IntSpectrum(G4 << 398 G4double G4eIonisationSpectrum::IntSpectrum(G4double xMin,
433 G4double xMax, 399 G4double xMax,
434 const G4DataVector& p) const 400 const G4DataVector& p) const
435 { 401 {
436 // Please comment what IntSpectrum does 402 // Please comment what IntSpectrum does
437 G4double sum = 0.0; 403 G4double sum = 0.0;
438 if(xMin >= xMax) return sum; 404 if(xMin >= xMax) return sum;
439 405
440 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, 406 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;
441 407
442 // Integral over interpolation aria 408 // Integral over interpolation aria
443 if(xMin < p[3]) { 409 if(xMin < p[3]) {
444 410
445 x1 = p[1]; 411 x1 = p[1];
446 y1 = p[4]; 412 y1 = p[4];
447 413
448 G4double dx = (p[2] - p[1]) / 3.0; 414 G4double dx = (p[2] - p[1]) / 3.0;
449 G4double dx1= G4Exp(std::log(p[3]/p[2]) / << 415 G4double dx1= exp(log(p[3]/p[2]) / 16.0);
450 416
451 for (size_t i=0; i<19; i++) { 417 for (size_t i=0; i<19; i++) {
452 418
453 q = 0.0; 419 q = 0.0;
454 if (i < 3) { 420 if (i < 3) {
455 x2 = x1 + dx; 421 x2 = x1 + dx;
456 } else if(18 == i) { 422 } else if(18 == i) {
457 x2 = p[3]; 423 x2 = p[3];
458 } else { 424 } else {
459 x2 = x1*dx1; 425 x2 = x1*dx1;
460 } 426 }
461 427
462 y2 = p[5 + i]; 428 y2 = p[5 + i];
463 429
464 if (xMax <= x1) { 430 if (xMax <= x1) {
465 break; 431 break;
466 } else if (xMin < x2) { 432 } else if (xMin < x2) {
467 433
468 xs1 = x1; 434 xs1 = x1;
469 xs2 = x2; 435 xs2 = x2;
470 ys1 = y1; 436 ys1 = y1;
471 ys2 = y2; 437 ys2 = y2;
472 438
473 if (x2 > x1) { 439 if (x2 > x1) {
474 if (xMin > x1) { 440 if (xMin > x1) {
475 xs1 = xMin; 441 xs1 = xMin;
476 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - 442 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
477 } << 443 }
478 if (xMax < x2) { 444 if (xMax < x2) {
479 xs2 = xMax; 445 xs2 = xMax;
480 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - 446 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
481 } << 447 }
482 if (xs2 > xs1) { 448 if (xs2 > xs1) {
483 q = (ys1*xs2 - ys2*xs1)/(xs1*xs2) << 449 q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
484 + std::log(xs2/xs1)*(ys2 - ys1) << 450 + log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
485 sum += q; 451 sum += q;
486 if(p.size() == 26) G4cout << "i= " 452 if(p.size() == 26) G4cout << "i= " << i << " q= " << q << " sum= " << sum << G4endl;
487 } 453 }
488 } << 454 }
489 } 455 }
490 x1 = x2; 456 x1 = x2;
491 y1 = y2; 457 y1 = y2;
492 } 458 }
493 } 459 }
494 460
495 // Integral over aria with parametrised form << 461 // Integral over aria with parametrised formula
496 462
497 x1 = std::max(xMin, p[3]); 463 x1 = std::max(xMin, p[3]);
498 if(x1 >= xMax) return sum; 464 if(x1 >= xMax) return sum;
499 x2 = xMax; 465 x2 = xMax;
500 466
501 xs1 = 1./x1; 467 xs1 = 1./x1;
502 xs2 = 1./x2; 468 xs2 = 1./x2;
503 q = (xs1 - xs2)*(1.0 - p[0]) << 469 q = (xs1 - xs2)*(1.0 - p[0])
504 - p[iMax]*std::log(x2/x1) << 470 - p[iMax]*log(x2/x1)
505 + (1. - p[iMax])*(x2 - x1) 471 + (1. - p[iMax])*(x2 - x1)
506 + 1./(1. - x2) - 1./(1. - x1) << 472 + 1./(1. - x2) - 1./(1. - x1)
507 + p[iMax]*std::log((1. - x2)/(1. - x1)) << 473 + p[iMax]*log((1. - x2)/(1. - x1))
508 + 0.25*p[0]*(xs1*xs1 - xs2*xs2); 474 + 0.25*p[0]*(xs1*xs1 - xs2*xs2);
509 sum += q; 475 sum += q;
510 if(p.size() == 26) G4cout << "param... q= " 476 if(p.size() == 26) G4cout << "param... q= " << q << " sum= " << sum << G4endl;
511 477
512 return sum; 478 return sum;
513 } << 479 }
514 480
515 481
516 G4double G4eIonisationSpectrum::AverageValue(G << 482 G4double G4eIonisationSpectrum::AverageValue(G4double xMin,
517 G4double xMax, 483 G4double xMax,
518 const G4DataVector& p) const 484 const G4DataVector& p) const
519 { 485 {
520 G4double sum = 0.0; 486 G4double sum = 0.0;
521 if(xMin >= xMax) return sum; 487 if(xMin >= xMax) return sum;
522 488
523 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2; 489 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;
524 490
525 // Integral over interpolation aria 491 // Integral over interpolation aria
526 if(xMin < p[3]) { 492 if(xMin < p[3]) {
527 493
528 x1 = p[1]; 494 x1 = p[1];
529 y1 = p[4]; 495 y1 = p[4];
530 496
531 G4double dx = (p[2] - p[1]) / 3.0; 497 G4double dx = (p[2] - p[1]) / 3.0;
532 G4double dx1= G4Exp(std::log(p[3]/p[2]) / << 498 G4double dx1= exp(log(p[3]/p[2]) / 16.0);
533 499
534 for (size_t i=0; i<19; i++) { 500 for (size_t i=0; i<19; i++) {
535 501
536 if (i < 3) { 502 if (i < 3) {
537 x2 = x1 + dx; 503 x2 = x1 + dx;
538 } else if(18 == i) { 504 } else if(18 == i) {
539 x2 = p[3]; 505 x2 = p[3];
540 } else { 506 } else {
541 x2 = x1*dx1; 507 x2 = x1*dx1;
542 } 508 }
543 509
544 y2 = p[5 + i]; 510 y2 = p[5 + i];
545 511
546 if (xMax <= x1) { 512 if (xMax <= x1) {
547 break; 513 break;
548 } else if (xMin < x2) { 514 } else if (xMin < x2) {
549 515
550 xs1 = x1; 516 xs1 = x1;
551 xs2 = x2; 517 xs2 = x2;
552 ys1 = y1; 518 ys1 = y1;
553 ys2 = y2; 519 ys2 = y2;
554 520
555 if (x2 > x1) { 521 if (x2 > x1) {
556 if (xMin > x1) { 522 if (xMin > x1) {
557 xs1 = xMin; 523 xs1 = xMin;
558 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - 524 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
559 } << 525 }
560 if (xMax < x2) { 526 if (xMax < x2) {
561 xs2 = xMax; 527 xs2 = xMax;
562 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - 528 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
563 } << 529 }
564 if (xs2 > xs1) { 530 if (xs2 > xs1) {
565 sum += std::log(xs2/xs1)*(ys1*xs2 << 531 sum += log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
566 + ys2 - ys1; 532 + ys2 - ys1;
567 } 533 }
568 } << 534 }
569 } 535 }
570 x1 = x2; 536 x1 = x2;
571 y1 = y2; 537 y1 = y2;
572 538
573 } 539 }
574 } 540 }
575 541
576 // Integral over aria with parametrised form << 542 // Integral over aria with parametrised formula
577 543
578 x1 = std::max(xMin, p[3]); 544 x1 = std::max(xMin, p[3]);
579 if(x1 >= xMax) return sum; 545 if(x1 >= xMax) return sum;
580 x2 = xMax; 546 x2 = xMax;
581 547
582 xs1 = 1./x1; 548 xs1 = 1./x1;
583 xs2 = 1./x2; 549 xs2 = 1./x2;
584 550
585 sum += std::log(x2/x1)*(1.0 - p[0]) << 551 sum += log(x2/x1)*(1.0 - p[0])
586 + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1) 552 + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)
587 + 1./(1. - x2) - 1./(1. - x1) << 553 + 1./(1. - x2) - 1./(1. - x1)
588 + (1. + p[iMax])*std::log((1. - x2)/(1 << 554 + (1. + p[iMax])*log((1. - x2)/(1. - x1))
589 + 0.5*p[0]*(xs1 - xs2); 555 + 0.5*p[0]*(xs1 - xs2);
590 556
591 return sum; 557 return sum;
592 } << 558 }
593 559
594 560
595 void G4eIonisationSpectrum::PrintData() const << 561 void G4eIonisationSpectrum::PrintData() const
596 { 562 {
597 theParam->PrintData(); 563 theParam->PrintData();
598 } 564 }
599 565
600 G4double G4eIonisationSpectrum::MaxEnergyOfSec 566 G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
601 G4int, // Z = 0, 567 G4int, // Z = 0,
602 const G4ParticleDefinition* 568 const G4ParticleDefinition* ) const
603 { << 569 { return 0.5 * kineticEnergy; };
604 return 0.5 * kineticEnergy; << 570
605 } << 571
>> 572
>> 573
606 574