| 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$
>> 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: G4RDeIonisationSpectrum 34 // File name: G4RDeIonisationSpectrum
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.2001 MGP Revision to improv 41 // 10.10.2001 MGP Revision to improve code quality and
40 // consistency with d 42 // consistency with design
41 // 02.11.2001 VI Optimize sampling 43 // 02.11.2001 VI Optimize sampling of energy
42 // 29.11.2001 VI New parametrisatio 44 // 29.11.2001 VI New parametrisation
43 // 19.04.2002 VI Add protection in 45 // 19.04.2002 VI Add protection in case of energy below binding
44 // 30.05.2002 VI Update to 24-param 46 // 30.05.2002 VI Update to 24-parameters data
45 // 11.07.2002 VI Fix in integration 47 // 11.07.2002 VI Fix in integration over spectrum
46 // 48 //
47 // ------------------------------------------- 49 // -------------------------------------------------------------------
48 // 50 //
49 51
50 #include "G4RDeIonisationSpectrum.hh" 52 #include "G4RDeIonisationSpectrum.hh"
51 #include "G4RDAtomicTransitionManager.hh" 53 #include "G4RDAtomicTransitionManager.hh"
52 #include "G4RDAtomicShell.hh" 54 #include "G4RDAtomicShell.hh"
53 #include "G4PhysicalConstants.hh" 55 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh" 56 #include "G4SystemOfUnits.hh"
55 #include "G4DataVector.hh" 57 #include "G4DataVector.hh"
56 #include "Randomize.hh" 58 #include "Randomize.hh"
57 59
58 60
59 G4RDeIonisationSpectrum::G4RDeIonisationSpectr 61 G4RDeIonisationSpectrum::G4RDeIonisationSpectrum():G4RDVEnergySpectrum(),
60 lowestE(0.1*eV), 62 lowestE(0.1*eV),
61 factor(1.3), 63 factor(1.3),
62 iMax(24), 64 iMax(24),
63 verbose(0) 65 verbose(0)
64 { 66 {
65 theParam = new G4RDeIonisationParameters(); 67 theParam = new G4RDeIonisationParameters();
66 } 68 }
67 69
68 70
69 G4RDeIonisationSpectrum::~G4RDeIonisationSpect 71 G4RDeIonisationSpectrum::~G4RDeIonisationSpectrum()
70 { 72 {
71 delete theParam; 73 delete theParam;
72 } 74 }
73 75
74 76
75 G4double G4RDeIonisationSpectrum::Probability( 77 G4double G4RDeIonisationSpectrum::Probability(G4int Z,
76 G4double tMin, 78 G4double tMin,
77 G4double tMax, 79 G4double tMax,
78 G4double e, 80 G4double e,
79 G4int shell, 81 G4int shell,
80 const G4ParticleDefinition* ) co 82 const G4ParticleDefinition* ) const
81 { 83 {
82 // Please comment what Probability does and 84 // Please comment what Probability does and what are the three
83 // functions mentioned below 85 // functions mentioned below
84 // Describe the algorithms used 86 // Describe the algorithms used
85 87
86 G4double eMax = MaxEnergyOfSecondaries(e); 88 G4double eMax = MaxEnergyOfSecondaries(e);
87 G4double t0 = std::max(tMin, lowestE); 89 G4double t0 = std::max(tMin, lowestE);
88 G4double tm = std::min(tMax, eMax); 90 G4double tm = std::min(tMax, eMax);
89 if(t0 >= tm) return 0.0; 91 if(t0 >= tm) return 0.0;
90 92
91 G4double bindingEnergy = (G4RDAtomicTransiti 93 G4double bindingEnergy = (G4RDAtomicTransitionManager::Instance())->
92 Shell(Z, shell)->Bi 94 Shell(Z, shell)->BindingEnergy();
93 95
94 if(e <= bindingEnergy) return 0.0; 96 if(e <= bindingEnergy) return 0.0;
95 97
96 G4double energy = e + bindingEnergy; 98 G4double energy = e + bindingEnergy;
97 99
98 G4double x1 = std::min(0.5,(t0 + bindingEner 100 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
99 G4double x2 = std::min(0.5,(tm + bindingEner 101 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
100 102
101 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0. 103 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
102 G4cout << "G4RDeIonisationSpectrum::Probab 104 G4cout << "G4RDeIonisationSpectrum::Probability: Z= " << Z
103 << "; shell= " << shell 105 << "; shell= " << shell
104 << "; E(keV)= " << e/keV 106 << "; E(keV)= " << e/keV
105 << "; Eb(keV)= " << bindingEnergy/k 107 << "; Eb(keV)= " << bindingEnergy/keV
106 << "; x1= " << x1 108 << "; x1= " << x1
107 << "; x2= " << x2 109 << "; x2= " << x2
108 << G4endl; 110 << G4endl;
109 111
110 } 112 }
111 113
112 G4DataVector p; 114 G4DataVector p;
113 115
114 // Access parameters 116 // Access parameters
115 for (G4int i=0; i<iMax; i++) 117 for (G4int i=0; i<iMax; i++)
116 { 118 {
117 G4double x = theParam->Parameter(Z, shell, 119 G4double x = theParam->Parameter(Z, shell, i, e);
118 if(i<4) x /= energy; 120 if(i<4) x /= energy;
119 p.push_back(x); 121 p.push_back(x);
120 } 122 }
121 123
122 if(p[3] > 0.5) p[3] = 0.5; 124 if(p[3] > 0.5) p[3] = 0.5;
123 125
124 G4double g = energy/electron_mass_c2 + 1.; 126 G4double g = energy/electron_mass_c2 + 1.;
125 p.push_back((2.0*g - 1.0)/(g*g)); 127 p.push_back((2.0*g - 1.0)/(g*g));
126 128
127 p[iMax-1] = Function(p[3], p); 129 p[iMax-1] = Function(p[3], p);
128 130
129 if(e >= 1. && e <= 0. && Z == 4) p.push_back 131 if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
130 132
131 133
132 G4double val = IntSpectrum(x1, x2, p); 134 G4double val = IntSpectrum(x1, x2, p);
133 G4double x0 = (lowestE + bindingEnergy)/ene 135 G4double x0 = (lowestE + bindingEnergy)/energy;
134 G4double nor = IntSpectrum(x0, 0.5, p); 136 G4double nor = IntSpectrum(x0, 0.5, p);
135 137
136 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0. 138 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
137 G4cout << "tcut= " << tMin 139 G4cout << "tcut= " << tMin
138 << "; tMax= " << tMax 140 << "; tMax= " << tMax
139 << "; x0= " << x0 141 << "; x0= " << x0
140 << "; x1= " << x1 142 << "; x1= " << x1
141 << "; x2= " << x2 143 << "; x2= " << x2
142 << "; val= " << val 144 << "; val= " << val
143 << "; nor= " << nor 145 << "; nor= " << nor
144 << "; sum= " << p[0] 146 << "; sum= " << p[0]
145 << "; a= " << p[1] 147 << "; a= " << p[1]
146 << "; b= " << p[2] 148 << "; b= " << p[2]
147 << "; c= " << p[3] 149 << "; c= " << p[3]
148 << G4endl; 150 << G4endl;
149 if(shell == 1) G4cout << "============" << 151 if(shell == 1) G4cout << "============" << G4endl;
150 } 152 }
151 153
152 p.clear(); 154 p.clear();
153 155
154 if(nor > 0.0) val /= nor; 156 if(nor > 0.0) val /= nor;
155 else val = 0.0; 157 else val = 0.0;
156 158
157 return val; 159 return val;
158 } 160 }
159 161
160 162
161 G4double G4RDeIonisationSpectrum::AverageEnerg 163 G4double G4RDeIonisationSpectrum::AverageEnergy(G4int Z,
162 G4double tMin, 164 G4double tMin,
163 G4double tMax, 165 G4double tMax,
164 G4double e, 166 G4double e,
165 G4int shell, 167 G4int shell,
166 const G4ParticleDefinition* ) 168 const G4ParticleDefinition* ) const
167 { 169 {
168 // Please comment what AverageEnergy does an 170 // Please comment what AverageEnergy does and what are the three
169 // functions mentioned below 171 // functions mentioned below
170 // Describe the algorithms used 172 // Describe the algorithms used
171 173
172 G4double eMax = MaxEnergyOfSecondaries(e); 174 G4double eMax = MaxEnergyOfSecondaries(e);
173 G4double t0 = std::max(tMin, lowestE); 175 G4double t0 = std::max(tMin, lowestE);
174 G4double tm = std::min(tMax, eMax); 176 G4double tm = std::min(tMax, eMax);
175 if(t0 >= tm) return 0.0; 177 if(t0 >= tm) return 0.0;
176 178
177 G4double bindingEnergy = (G4RDAtomicTransiti 179 G4double bindingEnergy = (G4RDAtomicTransitionManager::Instance())->
178 Shell(Z, shell)->Bi 180 Shell(Z, shell)->BindingEnergy();
179 181
180 if(e <= bindingEnergy) return 0.0; 182 if(e <= bindingEnergy) return 0.0;
181 183
182 G4double energy = e + bindingEnergy; 184 G4double energy = e + bindingEnergy;
183 185
184 G4double x1 = std::min(0.5,(t0 + bindingEner 186 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
185 G4double x2 = std::min(0.5,(tm + bindingEner 187 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
186 188
187 if(verbose > 1) { 189 if(verbose > 1) {
188 G4cout << "G4RDeIonisationSpectrum::Averag 190 G4cout << "G4RDeIonisationSpectrum::AverageEnergy: Z= " << Z
189 << "; shell= " << shell 191 << "; shell= " << shell
190 << "; E(keV)= " << e/keV 192 << "; E(keV)= " << e/keV
191 << "; bindingE(keV)= " << bindingEn 193 << "; bindingE(keV)= " << bindingEnergy/keV
192 << "; x1= " << x1 194 << "; x1= " << x1
193 << "; x2= " << x2 195 << "; x2= " << x2
194 << G4endl; 196 << G4endl;
195 } 197 }
196 198
197 G4DataVector p; 199 G4DataVector p;
198 200
199 // Access parameters 201 // Access parameters
200 for (G4int i=0; i<iMax; i++) 202 for (G4int i=0; i<iMax; i++)
201 { 203 {
202 G4double x = theParam->Parameter(Z, shell, 204 G4double x = theParam->Parameter(Z, shell, i, e);
203 if(i<4) x /= energy; 205 if(i<4) x /= energy;
204 p.push_back(x); 206 p.push_back(x);
205 } 207 }
206 208
207 if(p[3] > 0.5) p[3] = 0.5; 209 if(p[3] > 0.5) p[3] = 0.5;
208 210
209 G4double g = energy/electron_mass_c2 + 1.; 211 G4double g = energy/electron_mass_c2 + 1.;
210 p.push_back((2.0*g - 1.0)/(g*g)); 212 p.push_back((2.0*g - 1.0)/(g*g));
211 213
212 p[iMax-1] = Function(p[3], p); 214 p[iMax-1] = Function(p[3], p);
213 215
214 G4double val = AverageValue(x1, x2, p); 216 G4double val = AverageValue(x1, x2, p);
215 G4double x0 = (lowestE + bindingEnergy)/ene 217 G4double x0 = (lowestE + bindingEnergy)/energy;
216 G4double nor = IntSpectrum(x0, 0.5, p); 218 G4double nor = IntSpectrum(x0, 0.5, p);
217 val *= energy; 219 val *= energy;
218 220
219 if(verbose > 1) { 221 if(verbose > 1) {
220 G4cout << "tcut(MeV)= " << tMin/MeV 222 G4cout << "tcut(MeV)= " << tMin/MeV
221 << "; tMax(MeV)= " << tMax/MeV 223 << "; tMax(MeV)= " << tMax/MeV
222 << "; x0= " << x0 224 << "; x0= " << x0
223 << "; x1= " << x1 225 << "; x1= " << x1
224 << "; x2= " << x2 226 << "; x2= " << x2
225 << "; val= " << val 227 << "; val= " << val
226 << "; nor= " << nor 228 << "; nor= " << nor
227 << "; sum= " << p[0] 229 << "; sum= " << p[0]
228 << "; a= " << p[1] 230 << "; a= " << p[1]
229 << "; b= " << p[2] 231 << "; b= " << p[2]
230 << "; c= " << p[3] 232 << "; c= " << p[3]
231 << G4endl; 233 << G4endl;
232 } 234 }
233 235
234 p.clear(); 236 p.clear();
235 237
236 if(nor > 0.0) val /= nor; 238 if(nor > 0.0) val /= nor;
237 else val = 0.0; 239 else val = 0.0;
238 240
239 return val; 241 return val;
240 } 242 }
241 243
242 244
243 G4double G4RDeIonisationSpectrum::SampleEnergy 245 G4double G4RDeIonisationSpectrum::SampleEnergy(G4int Z,
244 G4double tMin, 246 G4double tMin,
245 G4double tMax, 247 G4double tMax,
246 G4double e, 248 G4double e,
247 G4int shell, 249 G4int shell,
248 const G4ParticleDefinition* ) c 250 const G4ParticleDefinition* ) const
249 { 251 {
250 // Please comment what SampleEnergy does 252 // Please comment what SampleEnergy does
251 G4double tDelta = 0.0; 253 G4double tDelta = 0.0;
252 G4double t0 = std::max(tMin, lowestE); 254 G4double t0 = std::max(tMin, lowestE);
253 G4double tm = std::min(tMax, MaxEnergyOfSeco 255 G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
254 if(t0 > tm) return tDelta; 256 if(t0 > tm) return tDelta;
255 257
256 G4double bindingEnergy = (G4RDAtomicTransiti 258 G4double bindingEnergy = (G4RDAtomicTransitionManager::Instance())->
257 Shell(Z, shell)->Bi 259 Shell(Z, shell)->BindingEnergy();
258 260
259 if(e <= bindingEnergy) return 0.0; 261 if(e <= bindingEnergy) return 0.0;
260 262
261 G4double energy = e + bindingEnergy; 263 G4double energy = e + bindingEnergy;
262 264
263 G4double x1 = std::min(0.5,(t0 + bindingEner 265 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
264 G4double x2 = std::min(0.5,(tm + bindingEner 266 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
265 if(x1 >= x2) return tDelta; 267 if(x1 >= x2) return tDelta;
266 268
267 if(verbose > 1) { 269 if(verbose > 1) {
268 G4cout << "G4RDeIonisationSpectrum::Sample 270 G4cout << "G4RDeIonisationSpectrum::SampleEnergy: Z= " << Z
269 << "; shell= " << shell 271 << "; shell= " << shell
270 << "; E(keV)= " << e/keV 272 << "; E(keV)= " << e/keV
271 << G4endl; 273 << G4endl;
272 } 274 }
273 275
274 // Access parameters 276 // Access parameters
275 G4DataVector p; 277 G4DataVector p;
276 278
277 // Access parameters 279 // Access parameters
278 for (G4int i=0; i<iMax; i++) 280 for (G4int i=0; i<iMax; i++)
279 { 281 {
280 G4double x = theParam->Parameter(Z, shell, 282 G4double x = theParam->Parameter(Z, shell, i, e);
281 if(i<4) x /= energy; 283 if(i<4) x /= energy;
282 p.push_back(x); 284 p.push_back(x);
283 } 285 }
284 286
285 if(p[3] > 0.5) p[3] = 0.5; 287 if(p[3] > 0.5) p[3] = 0.5;
286 288
287 G4double g = energy/electron_mass_c2 + 1.; 289 G4double g = energy/electron_mass_c2 + 1.;
288 p.push_back((2.0*g - 1.0)/(g*g)); 290 p.push_back((2.0*g - 1.0)/(g*g));
289 291
290 p[iMax-1] = Function(p[3], p); 292 p[iMax-1] = Function(p[3], p);
291 293
292 G4double aria1 = 0.0; 294 G4double aria1 = 0.0;
293 G4double a1 = std::max(x1,p[1]); 295 G4double a1 = std::max(x1,p[1]);
294 G4double a2 = std::min(x2,p[3]); 296 G4double a2 = std::min(x2,p[3]);
295 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p); 297 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
296 G4double aria2 = 0.0; 298 G4double aria2 = 0.0;
297 G4double a3 = std::max(x1,p[3]); 299 G4double a3 = std::max(x1,p[3]);
298 G4double a4 = x2; 300 G4double a4 = x2;
299 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p); 301 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
300 302
301 G4double aria = (aria1 + aria2)*G4UniformRan 303 G4double aria = (aria1 + aria2)*G4UniformRand();
302 G4double amaj, fun, q, x, z1, z2, dx, dx1; 304 G4double amaj, fun, q, x, z1, z2, dx, dx1;
303 305
304 //======= First aria to sample ===== 306 //======= First aria to sample =====
305 307
306 if(aria <= aria1) { 308 if(aria <= aria1) {
307 309
308 amaj = p[4]; 310 amaj = p[4];
309 for (G4int j=5; j<iMax; j++) { 311 for (G4int j=5; j<iMax; j++) {
310 if(p[j] > amaj) amaj = p[j]; 312 if(p[j] > amaj) amaj = p[j];
311 } 313 }
312 314
313 a1 = 1./a1; 315 a1 = 1./a1;
314 a2 = 1./a2; 316 a2 = 1./a2;
315 317
316 G4int i; 318 G4int i;
317 do { 319 do {
318 320
319 x = 1./(a2 + G4UniformRand()*(a1 - a2)); 321 x = 1./(a2 + G4UniformRand()*(a1 - a2));
320 z1 = p[1]; 322 z1 = p[1];
321 z2 = p[3]; 323 z2 = p[3];
322 dx = (p[2] - p[1]) / 3.0; 324 dx = (p[2] - p[1]) / 3.0;
323 dx1= std::exp(std::log(p[3]/p[2]) / 16.0 325 dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
324 for (i=4; i<iMax-1; i++) { 326 for (i=4; i<iMax-1; i++) {
325 327
326 if (i < 7) { 328 if (i < 7) {
327 z2 = z1 + dx; 329 z2 = z1 + dx;
328 } else if(iMax-2 == i) { 330 } else if(iMax-2 == i) {
329 z2 = p[3]; 331 z2 = p[3];
330 break; 332 break;
331 } else { 333 } else {
332 z2 = z1*dx1; 334 z2 = z1*dx1;
333 } 335 }
334 if(x >= z1 && x <= z2) break; 336 if(x >= z1 && x <= z2) break;
335 z1 = z2; 337 z1 = z2;
336 } 338 }
337 fun = p[i] + (x - z1) * (p[i+1] - p[i])/ 339 fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
338 340
339 if(fun > amaj) { 341 if(fun > amaj) {
340 G4cout << "WARNING in G4RDeIonisatio 342 G4cout << "WARNING in G4RDeIonisationSpectrum::SampleEnergy:"
341 << " Majoranta " << amaj 343 << " Majoranta " << amaj
342 << " < " << fun 344 << " < " << fun
343 << " in the first aria at x= 345 << " in the first aria at x= " << x
344 << G4endl; 346 << G4endl;
345 } 347 }
346 348
347 q = amaj*G4UniformRand(); 349 q = amaj*G4UniformRand();
348 350
349 } while (q >= fun); 351 } while (q >= fun);
350 352
351 //======= Second aria to sample ===== 353 //======= Second aria to sample =====
352 354
353 } else { 355 } else {
354 356
355 amaj = std::max(p[iMax-1], Function(0.5, p 357 amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
356 a1 = 1./a3; 358 a1 = 1./a3;
357 a2 = 1./a4; 359 a2 = 1./a4;
358 360
359 do { 361 do {
360 362
361 x = 1./(a2 + G4UniformRand()*(a1 - a2)); 363 x = 1./(a2 + G4UniformRand()*(a1 - a2));
362 fun = Function(x, p); 364 fun = Function(x, p);
363 365
364 if(fun > amaj) { 366 if(fun > amaj) {
365 G4cout << "WARNING in G4RDeIonisatio 367 G4cout << "WARNING in G4RDeIonisationSpectrum::SampleEnergy:"
366 << " Majoranta " << amaj 368 << " Majoranta " << amaj
367 << " < " << fun 369 << " < " << fun
368 << " in the second aria at x= 370 << " in the second aria at x= " << x
369 << G4endl; 371 << G4endl;
370 } 372 }
371 373
372 q = amaj*G4UniformRand(); 374 q = amaj*G4UniformRand();
373 375
374 } while (q >= fun); 376 } while (q >= fun);
375 377
376 } 378 }
377 379
378 p.clear(); 380 p.clear();
379 381
380 tDelta = x*energy - bindingEnergy; 382 tDelta = x*energy - bindingEnergy;
381 383
382 if(verbose > 1) { 384 if(verbose > 1) {
383 G4cout << "tcut(MeV)= " << tMin/MeV 385 G4cout << "tcut(MeV)= " << tMin/MeV
384 << "; tMax(MeV)= " << tMax/MeV 386 << "; tMax(MeV)= " << tMax/MeV
385 << "; x1= " << x1 387 << "; x1= " << x1
386 << "; x2= " << x2 388 << "; x2= " << x2
387 << "; a1= " << a1 389 << "; a1= " << a1
388 << "; a2= " << a2 390 << "; a2= " << a2
389 << "; x= " << x 391 << "; x= " << x
390 << "; be= " << bindingEnergy 392 << "; be= " << bindingEnergy
391 << "; e= " << e 393 << "; e= " << e
392 << "; tDelta= " << tDelta 394 << "; tDelta= " << tDelta
393 << G4endl; 395 << G4endl;
394 } 396 }
395 397
396 398
397 return tDelta; 399 return tDelta;
398 } 400 }
399 401
400 402
401 G4double G4RDeIonisationSpectrum::IntSpectrum( 403 G4double G4RDeIonisationSpectrum::IntSpectrum(G4double xMin,
402 G4double xMax, 404 G4double xMax,
403 const G4DataVector& p) const 405 const G4DataVector& p) const
404 { 406 {
405 // Please comment what IntSpectrum does 407 // Please comment what IntSpectrum does
406 G4double sum = 0.0; 408 G4double sum = 0.0;
407 if(xMin >= xMax) return sum; 409 if(xMin >= xMax) return sum;
408 410
409 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, 411 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;
410 412
411 // Integral over interpolation aria 413 // Integral over interpolation aria
412 if(xMin < p[3]) { 414 if(xMin < p[3]) {
413 415
414 x1 = p[1]; 416 x1 = p[1];
415 y1 = p[4]; 417 y1 = p[4];
416 418
417 G4double dx = (p[2] - p[1]) / 3.0; 419 G4double dx = (p[2] - p[1]) / 3.0;
418 G4double dx1= std::exp(std::log(p[3]/p[2]) 420 G4double dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
419 421
420 for (size_t i=0; i<19; i++) { 422 for (size_t i=0; i<19; i++) {
421 423
422 q = 0.0; 424 q = 0.0;
423 if (i < 3) { 425 if (i < 3) {
424 x2 = x1 + dx; 426 x2 = x1 + dx;
425 } else if(18 == i) { 427 } else if(18 == i) {
426 x2 = p[3]; 428 x2 = p[3];
427 } else { 429 } else {
428 x2 = x1*dx1; 430 x2 = x1*dx1;
429 } 431 }
430 432
431 y2 = p[5 + i]; 433 y2 = p[5 + i];
432 434
433 if (xMax <= x1) { 435 if (xMax <= x1) {
434 break; 436 break;
435 } else if (xMin < x2) { 437 } else if (xMin < x2) {
436 438
437 xs1 = x1; 439 xs1 = x1;
438 xs2 = x2; 440 xs2 = x2;
439 ys1 = y1; 441 ys1 = y1;
440 ys2 = y2; 442 ys2 = y2;
441 443
442 if (x2 > x1) { 444 if (x2 > x1) {
443 if (xMin > x1) { 445 if (xMin > x1) {
444 xs1 = xMin; 446 xs1 = xMin;
445 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - 447 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
446 } 448 }
447 if (xMax < x2) { 449 if (xMax < x2) {
448 xs2 = xMax; 450 xs2 = xMax;
449 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - 451 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
450 } 452 }
451 if (xs2 > xs1) { 453 if (xs2 > xs1) {
452 q = (ys1*xs2 - ys2*xs1)/(xs1*xs2) 454 q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
453 + std::log(xs2/xs1)*(ys2 - ys1) 455 + std::log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
454 sum += q; 456 sum += q;
455 if(p.size() == 26) G4cout << "i= " 457 if(p.size() == 26) G4cout << "i= " << i << " q= " << q << " sum= " << sum << G4endl;
456 } 458 }
457 } 459 }
458 } 460 }
459 x1 = x2; 461 x1 = x2;
460 y1 = y2; 462 y1 = y2;
461 } 463 }
462 } 464 }
463 465
464 // Integral over aria with parametrised form 466 // Integral over aria with parametrised formula
465 467
466 x1 = std::max(xMin, p[3]); 468 x1 = std::max(xMin, p[3]);
467 if(x1 >= xMax) return sum; 469 if(x1 >= xMax) return sum;
468 x2 = xMax; 470 x2 = xMax;
469 471
470 xs1 = 1./x1; 472 xs1 = 1./x1;
471 xs2 = 1./x2; 473 xs2 = 1./x2;
472 q = (xs1 - xs2)*(1.0 - p[0]) 474 q = (xs1 - xs2)*(1.0 - p[0])
473 - p[iMax]*std::log(x2/x1) 475 - p[iMax]*std::log(x2/x1)
474 + (1. - p[iMax])*(x2 - x1) 476 + (1. - p[iMax])*(x2 - x1)
475 + 1./(1. - x2) - 1./(1. - x1) 477 + 1./(1. - x2) - 1./(1. - x1)
476 + p[iMax]*std::log((1. - x2)/(1. - x1)) 478 + p[iMax]*std::log((1. - x2)/(1. - x1))
477 + 0.25*p[0]*(xs1*xs1 - xs2*xs2); 479 + 0.25*p[0]*(xs1*xs1 - xs2*xs2);
478 sum += q; 480 sum += q;
479 if(p.size() == 26) G4cout << "param... q= " 481 if(p.size() == 26) G4cout << "param... q= " << q << " sum= " << sum << G4endl;
480 482
481 return sum; 483 return sum;
482 } 484 }
483 485
484 486
485 G4double G4RDeIonisationSpectrum::AverageValue 487 G4double G4RDeIonisationSpectrum::AverageValue(G4double xMin,
486 G4double xMax, 488 G4double xMax,
487 const G4DataVector& p) const 489 const G4DataVector& p) const
488 { 490 {
489 G4double sum = 0.0; 491 G4double sum = 0.0;
490 if(xMin >= xMax) return sum; 492 if(xMin >= xMax) return sum;
491 493
492 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2; 494 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;
493 495
494 // Integral over interpolation aria 496 // Integral over interpolation aria
495 if(xMin < p[3]) { 497 if(xMin < p[3]) {
496 498
497 x1 = p[1]; 499 x1 = p[1];
498 y1 = p[4]; 500 y1 = p[4];
499 501
500 G4double dx = (p[2] - p[1]) / 3.0; 502 G4double dx = (p[2] - p[1]) / 3.0;
501 G4double dx1= std::exp(std::log(p[3]/p[2]) 503 G4double dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
502 504
503 for (size_t i=0; i<19; i++) { 505 for (size_t i=0; i<19; i++) {
504 506
505 if (i < 3) { 507 if (i < 3) {
506 x2 = x1 + dx; 508 x2 = x1 + dx;
507 } else if(18 == i) { 509 } else if(18 == i) {
508 x2 = p[3]; 510 x2 = p[3];
509 } else { 511 } else {
510 x2 = x1*dx1; 512 x2 = x1*dx1;
511 } 513 }
512 514
513 y2 = p[5 + i]; 515 y2 = p[5 + i];
514 516
515 if (xMax <= x1) { 517 if (xMax <= x1) {
516 break; 518 break;
517 } else if (xMin < x2) { 519 } else if (xMin < x2) {
518 520
519 xs1 = x1; 521 xs1 = x1;
520 xs2 = x2; 522 xs2 = x2;
521 ys1 = y1; 523 ys1 = y1;
522 ys2 = y2; 524 ys2 = y2;
523 525
524 if (x2 > x1) { 526 if (x2 > x1) {
525 if (xMin > x1) { 527 if (xMin > x1) {
526 xs1 = xMin; 528 xs1 = xMin;
527 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - 529 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
528 } 530 }
529 if (xMax < x2) { 531 if (xMax < x2) {
530 xs2 = xMax; 532 xs2 = xMax;
531 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - 533 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
532 } 534 }
533 if (xs2 > xs1) { 535 if (xs2 > xs1) {
534 sum += std::log(xs2/xs1)*(ys1*xs2 536 sum += std::log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
535 + ys2 - ys1; 537 + ys2 - ys1;
536 } 538 }
537 } 539 }
538 } 540 }
539 x1 = x2; 541 x1 = x2;
540 y1 = y2; 542 y1 = y2;
541 543
542 } 544 }
543 } 545 }
544 546
545 // Integral over aria with parametrised form 547 // Integral over aria with parametrised formula
546 548
547 x1 = std::max(xMin, p[3]); 549 x1 = std::max(xMin, p[3]);
548 if(x1 >= xMax) return sum; 550 if(x1 >= xMax) return sum;
549 x2 = xMax; 551 x2 = xMax;
550 552
551 xs1 = 1./x1; 553 xs1 = 1./x1;
552 xs2 = 1./x2; 554 xs2 = 1./x2;
553 555
554 sum += std::log(x2/x1)*(1.0 - p[0]) 556 sum += std::log(x2/x1)*(1.0 - p[0])
555 + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1) 557 + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)
556 + 1./(1. - x2) - 1./(1. - x1) 558 + 1./(1. - x2) - 1./(1. - x1)
557 + (1. + p[iMax])*std::log((1. - x2)/(1 559 + (1. + p[iMax])*std::log((1. - x2)/(1. - x1))
558 + 0.5*p[0]*(xs1 - xs2); 560 + 0.5*p[0]*(xs1 - xs2);
559 561
560 return sum; 562 return sum;
561 } 563 }
562 564
563 565
564 void G4RDeIonisationSpectrum::PrintData() cons 566 void G4RDeIonisationSpectrum::PrintData() const
565 { 567 {
566 theParam->PrintData(); 568 theParam->PrintData();
567 } 569 }
568 570
569 G4double G4RDeIonisationSpectrum::MaxEnergyOfS 571 G4double G4RDeIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
570 G4int, // Z = 0, 572 G4int, // Z = 0,
571 const G4ParticleDefinition* 573 const G4ParticleDefinition* ) const
572 { 574 {
573 return 0.5 * kineticEnergy; 575 return 0.5 * kineticEnergy;
574 } 576 }
575 577