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 | // |
---|
27 | // $Id: G4InitXscPAI.cc,v 1.9 2006/06/29 19:53:00 gunter Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-02-ref-02 $ |
---|
29 | // |
---|
30 | // |
---|
31 | // G4InitXscPAI.cc -- class implementation file |
---|
32 | // |
---|
33 | // GEANT 4 class implementation file |
---|
34 | // |
---|
35 | // For information related to this code, please, contact |
---|
36 | // the Geant4 Collaboration. |
---|
37 | // |
---|
38 | // R&D: Vladimir.Grichine@cern.ch |
---|
39 | // |
---|
40 | // History: |
---|
41 | // |
---|
42 | |
---|
43 | |
---|
44 | |
---|
45 | #include "G4InitXscPAI.hh" |
---|
46 | |
---|
47 | #include "globals.hh" |
---|
48 | #include "G4ios.hh" |
---|
49 | #include "G4Poisson.hh" |
---|
50 | #include "G4Integrator.hh" |
---|
51 | #include "G4Material.hh" |
---|
52 | #include "G4MaterialCutsCouple.hh" |
---|
53 | #include "G4SandiaTable.hh" |
---|
54 | |
---|
55 | |
---|
56 | |
---|
57 | // Local class constants |
---|
58 | |
---|
59 | const G4double G4InitXscPAI::fDelta = 0.005 ; // energy shift from interval border |
---|
60 | const G4int G4InitXscPAI::fPAIbin = 100 ; // size of energy transfer vectors |
---|
61 | const G4double G4InitXscPAI::fSolidDensity = 0.05*g/cm3 ; // ~gas-solid border |
---|
62 | |
---|
63 | ////////////////////////////////////////////////////////////////// |
---|
64 | // |
---|
65 | // Constructor |
---|
66 | // |
---|
67 | |
---|
68 | using namespace std; |
---|
69 | |
---|
70 | G4InitXscPAI::G4InitXscPAI( const G4MaterialCutsCouple* matCC) |
---|
71 | : fPAIxscVector(NULL), |
---|
72 | fPAIdEdxVector(NULL), |
---|
73 | fPAIphotonVector(NULL), |
---|
74 | fPAIelectronVector(NULL), |
---|
75 | fChCosSqVector(NULL), |
---|
76 | fChWidthVector(NULL) |
---|
77 | { |
---|
78 | G4int i, j, matIndex; |
---|
79 | |
---|
80 | fDensity = matCC->GetMaterial()->GetDensity(); |
---|
81 | fElectronDensity = matCC->GetMaterial()->GetElectronDensity(); |
---|
82 | matIndex = matCC->GetMaterial()->GetIndex(); |
---|
83 | |
---|
84 | fSandia = new G4SandiaTable(matIndex); |
---|
85 | fIntervalNumber = fSandia->GetMaxInterval()-1; |
---|
86 | |
---|
87 | fMatSandiaMatrix = new G4OrderedTable(); |
---|
88 | |
---|
89 | for (i = 0; i < fIntervalNumber; i++) |
---|
90 | { |
---|
91 | fMatSandiaMatrix->push_back(new G4DataVector(5,0.)); |
---|
92 | } |
---|
93 | for (G4int i = 0; i < fIntervalNumber; i++) |
---|
94 | { |
---|
95 | (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0); |
---|
96 | |
---|
97 | for(j = 1; j < 5 ; j++) |
---|
98 | { |
---|
99 | (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity; |
---|
100 | } |
---|
101 | } |
---|
102 | KillCloseIntervals(); |
---|
103 | Normalisation(); |
---|
104 | |
---|
105 | } |
---|
106 | |
---|
107 | |
---|
108 | |
---|
109 | |
---|
110 | //////////////////////////////////////////////////////////////////////////// |
---|
111 | // |
---|
112 | // Destructor |
---|
113 | |
---|
114 | G4InitXscPAI::~G4InitXscPAI() |
---|
115 | { |
---|
116 | if(fPAIxscVector) delete fPAIxscVector; |
---|
117 | if(fPAIdEdxVector) delete fPAIdEdxVector; |
---|
118 | if(fPAIphotonVector) delete fPAIphotonVector; |
---|
119 | if(fPAIelectronVector) delete fPAIelectronVector; |
---|
120 | if(fChCosSqVector) delete fChCosSqVector; |
---|
121 | if(fChWidthVector) delete fChWidthVector; |
---|
122 | } |
---|
123 | |
---|
124 | //////////////////////////////////////////////////////////////////////// |
---|
125 | // |
---|
126 | // Kill close intervals, recalculate fIntervalNumber |
---|
127 | |
---|
128 | void G4InitXscPAI::KillCloseIntervals() |
---|
129 | { |
---|
130 | G4int i, j, k; |
---|
131 | G4double energy1, energy2; |
---|
132 | |
---|
133 | for( i = 0 ; i < fIntervalNumber - 1 ; i++ ) |
---|
134 | { |
---|
135 | energy1 = (*(*fMatSandiaMatrix)[i])[0]; |
---|
136 | energy2 = (*(*fMatSandiaMatrix)[i+1])[0]; |
---|
137 | |
---|
138 | if( energy2 - energy1 > 1.5*fDelta*(energy1 + energy2) ) continue ; |
---|
139 | else |
---|
140 | { |
---|
141 | for(j = i; j < fIntervalNumber-1; j++) |
---|
142 | { |
---|
143 | for( k = 0; k < 5; k++ ) |
---|
144 | { |
---|
145 | (*(*fMatSandiaMatrix)[j])[k] = (*(*fMatSandiaMatrix)[j+1])[k]; |
---|
146 | } |
---|
147 | } |
---|
148 | fIntervalNumber-- ; |
---|
149 | i-- ; |
---|
150 | } |
---|
151 | } |
---|
152 | |
---|
153 | } |
---|
154 | |
---|
155 | //////////////////////////////////////////////////////////////////////// |
---|
156 | // |
---|
157 | // Kill close intervals, recalculate fIntervalNumber |
---|
158 | |
---|
159 | void G4InitXscPAI::Normalisation() |
---|
160 | { |
---|
161 | G4int i, j; |
---|
162 | G4double energy1, energy2, delta, cof; // , shift; |
---|
163 | |
---|
164 | energy1 = (*(*fMatSandiaMatrix)[fIntervalNumber-1])[0]; |
---|
165 | energy2 = 2.*(*(*fMatSandiaMatrix)[fIntervalNumber-1])[0]; |
---|
166 | |
---|
167 | |
---|
168 | cof = RutherfordIntegral(fIntervalNumber-1,energy1,energy2); |
---|
169 | |
---|
170 | for( i = fIntervalNumber-2; i >= 0; i-- ) |
---|
171 | { |
---|
172 | energy1 = (*(*fMatSandiaMatrix)[i])[0]; |
---|
173 | energy2 = (*(*fMatSandiaMatrix)[i+1])[0]; |
---|
174 | |
---|
175 | cof += RutherfordIntegral(i,energy1,energy2); |
---|
176 | // G4cout<<"norm. cof = "<<cof<<G4endl; |
---|
177 | } |
---|
178 | fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ; |
---|
179 | fNormalizationCof *= fElectronDensity; |
---|
180 | delta = fNormalizationCof - cof; |
---|
181 | fNormalizationCof /= cof; |
---|
182 | // G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof |
---|
183 | // <<"; at delta ="<<delta<<G4endl ; |
---|
184 | |
---|
185 | for (G4int i = 0; i < fIntervalNumber; i++) // renormalisation on QM sum rule |
---|
186 | { |
---|
187 | for(j = 1; j < 5 ; j++) |
---|
188 | { |
---|
189 | (*(*fMatSandiaMatrix)[i])[j] *= fNormalizationCof; |
---|
190 | } |
---|
191 | } |
---|
192 | /* |
---|
193 | if(delta > 0) // shift the first energy interval |
---|
194 | { |
---|
195 | for(i=1;i<100;i++) |
---|
196 | { |
---|
197 | energy1 = (1.-i/100.)*(*(*fMatSandiaMatrix)[0])[0]; |
---|
198 | energy2 = (*(*fMatSandiaMatrix)[0])[0]; |
---|
199 | shift = RutherfordIntegral(0,energy1,energy2); |
---|
200 | G4cout<<shift<<"\t"; |
---|
201 | if(shift >= delta) break; |
---|
202 | } |
---|
203 | (*(*fMatSandiaMatrix)[0])[0] = energy1; |
---|
204 | cof += shift; |
---|
205 | } |
---|
206 | else if(delta < 0) |
---|
207 | { |
---|
208 | for(i=1;i<100;i++) |
---|
209 | { |
---|
210 | energy1 = (*(*fMatSandiaMatrix)[0])[0]; |
---|
211 | energy2 = (*(*fMatSandiaMatrix)[0])[0] + |
---|
212 | ( (*(*fMatSandiaMatrix)[0])[0] - (*(*fMatSandiaMatrix)[0])[0] )*i/100.; |
---|
213 | shift = RutherfordIntegral(0,energy1,energy2); |
---|
214 | if( shift >= std::abs(delta) ) break; |
---|
215 | } |
---|
216 | (*(*fMatSandiaMatrix)[0])[0] = energy2; |
---|
217 | cof -= shift; |
---|
218 | } |
---|
219 | G4cout<<G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof/cof |
---|
220 | <<"; at delta ="<<delta<<" and i = "<<i<<G4endl ; |
---|
221 | */ |
---|
222 | } |
---|
223 | |
---|
224 | |
---|
225 | |
---|
226 | |
---|
227 | |
---|
228 | //////////////////////////////////////////////////////////////////// |
---|
229 | // |
---|
230 | // Integration over electrons that could be considered |
---|
231 | // quasi-free at energy transfer of interest |
---|
232 | |
---|
233 | G4double G4InitXscPAI::RutherfordIntegral( G4int k, |
---|
234 | G4double x1, |
---|
235 | G4double x2 ) |
---|
236 | { |
---|
237 | G4double c1, c2, c3, a1, a2, a3, a4 ; |
---|
238 | |
---|
239 | a1 = (*(*fMatSandiaMatrix)[k])[1]; |
---|
240 | a2 = (*(*fMatSandiaMatrix)[k])[2]; |
---|
241 | a3 = (*(*fMatSandiaMatrix)[k])[3]; |
---|
242 | a4 = (*(*fMatSandiaMatrix)[k])[4]; |
---|
243 | // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; |
---|
244 | c1 = (x2 - x1)/x1/x2 ; |
---|
245 | c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ; |
---|
246 | c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ; |
---|
247 | // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl; |
---|
248 | |
---|
249 | return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ; |
---|
250 | |
---|
251 | } // end of RutherfordIntegral |
---|
252 | |
---|
253 | /////////////////////////////////////////////////////////////// |
---|
254 | // |
---|
255 | // Integrate photo-absorption cross-section from I1 up to omega |
---|
256 | |
---|
257 | G4double G4InitXscPAI::IntegralTerm(G4double omega) |
---|
258 | { |
---|
259 | G4int i; |
---|
260 | G4double energy1, energy2, result = 0.; |
---|
261 | |
---|
262 | for( i = 0; i <= fIntervalTmax; i++ ) |
---|
263 | { |
---|
264 | if(i == fIntervalTmax) |
---|
265 | { |
---|
266 | energy1 = (*(*fMatSandiaMatrix)[i])[0]; |
---|
267 | result += RutherfordIntegral(i,energy1,omega); |
---|
268 | } |
---|
269 | else |
---|
270 | { |
---|
271 | if( omega <= (*(*fMatSandiaMatrix)[i+1])[0]) |
---|
272 | { |
---|
273 | energy1 = (*(*fMatSandiaMatrix)[i])[0]; |
---|
274 | result += RutherfordIntegral(i,energy1,omega); |
---|
275 | break; |
---|
276 | } |
---|
277 | else |
---|
278 | { |
---|
279 | energy1 = (*(*fMatSandiaMatrix)[i])[0]; |
---|
280 | energy2 = (*(*fMatSandiaMatrix)[i+1])[0]; |
---|
281 | result += RutherfordIntegral(i,energy1,energy2); |
---|
282 | } |
---|
283 | } |
---|
284 | // G4cout<<"IntegralTerm<<"("<<omega<<")"<<" = "<<result<<G4endl; |
---|
285 | } |
---|
286 | return result; |
---|
287 | } |
---|
288 | |
---|
289 | |
---|
290 | //////////////////////////////////////////////////////////////// |
---|
291 | // |
---|
292 | // Imaginary part of dielectric constant |
---|
293 | // (G4int k - interval number, G4double en1 - energy point) |
---|
294 | |
---|
295 | G4double G4InitXscPAI::ImPartDielectricConst( G4int k , |
---|
296 | G4double energy1 ) |
---|
297 | { |
---|
298 | G4double energy2,energy3,energy4,a1,a2,a3,a4,result; |
---|
299 | |
---|
300 | a1 = (*(*fMatSandiaMatrix)[k])[1]; |
---|
301 | a2 = (*(*fMatSandiaMatrix)[k])[2]; |
---|
302 | a3 = (*(*fMatSandiaMatrix)[k])[3]; |
---|
303 | a4 = (*(*fMatSandiaMatrix)[k])[4]; |
---|
304 | |
---|
305 | energy2 = energy1*energy1; |
---|
306 | energy3 = energy2*energy1; |
---|
307 | energy4 = energy3*energy1; |
---|
308 | |
---|
309 | result = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ; |
---|
310 | result *= hbarc/energy1 ; |
---|
311 | |
---|
312 | return result ; |
---|
313 | |
---|
314 | } // end of ImPartDielectricConst |
---|
315 | |
---|
316 | //////////////////////////////////////////////////////////////// |
---|
317 | // |
---|
318 | // Modulus squared of dielectric constant |
---|
319 | // (G4int k - interval number, G4double omega - energy point) |
---|
320 | |
---|
321 | G4double G4InitXscPAI::ModuleSqDielectricConst( G4int k , |
---|
322 | G4double omega ) |
---|
323 | { |
---|
324 | G4double eIm2, eRe2, result; |
---|
325 | |
---|
326 | result = ImPartDielectricConst(k,omega); |
---|
327 | eIm2 = result*result; |
---|
328 | |
---|
329 | result = RePartDielectricConst(omega); |
---|
330 | eRe2 = result*result; |
---|
331 | |
---|
332 | result = eIm2 + eRe2; |
---|
333 | |
---|
334 | return result ; |
---|
335 | } |
---|
336 | |
---|
337 | |
---|
338 | ////////////////////////////////////////////////////////////////////////////// |
---|
339 | // |
---|
340 | // Real part of dielectric constant minus unit: epsilon_1 - 1 |
---|
341 | // (G4double enb - energy point) |
---|
342 | // |
---|
343 | |
---|
344 | G4double G4InitXscPAI::RePartDielectricConst(G4double enb) |
---|
345 | { |
---|
346 | G4int i; |
---|
347 | G4double x0, x02, x03, x04, x05, x1, x2, a1,a2,a3,a4,xx1 ,xx2 , xx12, |
---|
348 | c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ; |
---|
349 | |
---|
350 | x0 = enb ; |
---|
351 | result = 0 ; |
---|
352 | |
---|
353 | for( i = 0; i < fIntervalNumber-1; i++) |
---|
354 | { |
---|
355 | x1 = (*(*fMatSandiaMatrix)[i])[0]; |
---|
356 | x2 = (*(*fMatSandiaMatrix)[i+1])[0] ; |
---|
357 | |
---|
358 | a1 = (*(*fMatSandiaMatrix)[i])[1]; |
---|
359 | a2 = (*(*fMatSandiaMatrix)[i])[2]; |
---|
360 | a3 = (*(*fMatSandiaMatrix)[i])[3]; |
---|
361 | a4 = (*(*fMatSandiaMatrix)[i])[4]; |
---|
362 | |
---|
363 | if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta ) |
---|
364 | { |
---|
365 | if(x0 >= x1) x0 = x1*(1+fDelta); |
---|
366 | else x0 = x1*(1-fDelta); |
---|
367 | } |
---|
368 | if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta ) |
---|
369 | { |
---|
370 | if(x0 >= x2) x0 = x2*(1+fDelta); |
---|
371 | else x0 = x2*(1-fDelta); |
---|
372 | } |
---|
373 | xx1 = x1 - x0 ; |
---|
374 | xx2 = x2 - x0 ; |
---|
375 | xx12 = xx2/xx1 ; |
---|
376 | |
---|
377 | if( xx12 < 0 ) xx12 = -xx12; |
---|
378 | |
---|
379 | xln1 = log(x2/x1) ; |
---|
380 | xln2 = log(xx12) ; |
---|
381 | xln3 = log((x2 + x0)/(x1 + x0)) ; |
---|
382 | |
---|
383 | x02 = x0*x0 ; |
---|
384 | x03 = x02*x0 ; |
---|
385 | x04 = x03*x0 ; |
---|
386 | x05 = x04*x0; |
---|
387 | |
---|
388 | c1 = (x2 - x1)/x1/x2 ; |
---|
389 | c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ; |
---|
390 | c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ; |
---|
391 | |
---|
392 | result -= (a1/x02 + a3/x04)*xln1 ; |
---|
393 | result -= (a2/x02 + a4/x04)*c1 ; |
---|
394 | result -= a3*c2/2/x02 ; |
---|
395 | result -= a4*c3/3/x02 ; |
---|
396 | |
---|
397 | cof1 = a1/x02 + a3/x04 ; |
---|
398 | cof2 = a2/x03 + a4/x05 ; |
---|
399 | |
---|
400 | result += 0.5*(cof1 +cof2)*xln2 ; |
---|
401 | result += 0.5*(cof1 - cof2)*xln3 ; |
---|
402 | } |
---|
403 | result *= 2*hbarc/pi ; |
---|
404 | |
---|
405 | return result ; |
---|
406 | |
---|
407 | } // end of RePartDielectricConst |
---|
408 | |
---|
409 | ////////////////////////////////////////////////////////////////////// |
---|
410 | // |
---|
411 | // PAI differential cross-section in terms of |
---|
412 | // simplified Allison's equation |
---|
413 | // |
---|
414 | |
---|
415 | G4double G4InitXscPAI::DifPAIxSection( G4double omega ) |
---|
416 | { |
---|
417 | G4int i = fCurrentInterval; |
---|
418 | G4double betaGammaSq = fBetaGammaSq; |
---|
419 | G4double integralTerm = IntegralTerm(omega); |
---|
420 | G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ; |
---|
421 | G4double epsilonRe = RePartDielectricConst(omega); |
---|
422 | G4double epsilonIm = ImPartDielectricConst(i,omega); |
---|
423 | G4double be4 ; |
---|
424 | G4double betaBohr2 = fine_structure_const*fine_structure_const ; |
---|
425 | G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ; |
---|
426 | be2 = betaGammaSq/(1 + betaGammaSq) ; |
---|
427 | be4 = be2*be2 ; |
---|
428 | |
---|
429 | cof = 1 ; |
---|
430 | x1 = log(2*electron_mass_c2/omega) ; |
---|
431 | |
---|
432 | if( betaGammaSq < 0.01 ) x2 = log(be2) ; |
---|
433 | else |
---|
434 | { |
---|
435 | x2 = -log( (1/betaGammaSq - epsilonRe)* |
---|
436 | (1/betaGammaSq - epsilonRe) + |
---|
437 | epsilonIm*epsilonIm )/2 ; |
---|
438 | } |
---|
439 | if( epsilonIm == 0.0 || betaGammaSq < 0.01 ) |
---|
440 | { |
---|
441 | x6=0 ; |
---|
442 | } |
---|
443 | else |
---|
444 | { |
---|
445 | x3 = -epsilonRe + 1/betaGammaSq ; |
---|
446 | x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) + |
---|
447 | epsilonIm*epsilonIm) ; |
---|
448 | |
---|
449 | x7 = atan2(epsilonIm,x3) ; |
---|
450 | x6 = x5 * x7 ; |
---|
451 | } |
---|
452 | // if(fImPartDielectricConst[i] == 0) x6 = 0 ; |
---|
453 | |
---|
454 | x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ; |
---|
455 | // if( x4 < 0.0 ) x4 = 0.0 ; |
---|
456 | x8 = (1 + epsilonRe)*(1 + epsilonRe) + |
---|
457 | epsilonIm*epsilonIm; |
---|
458 | |
---|
459 | result = (x4 + cof*integralTerm/omega/omega) ; |
---|
460 | if(result < 1.0e-8) result = 1.0e-8 ; |
---|
461 | result *= fine_structure_const/be2/pi ; |
---|
462 | // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ; |
---|
463 | // result *= (1-exp(-be2/betaBohr2)) ; |
---|
464 | result *= (1-exp(-be4/betaBohr4)) ; |
---|
465 | if(fDensity >= fSolidDensity) |
---|
466 | { |
---|
467 | result /= x8 ; |
---|
468 | } |
---|
469 | return result ; |
---|
470 | |
---|
471 | } // end of DifPAIxSection |
---|
472 | |
---|
473 | ////////////////////////////////////////////////////////////////////// |
---|
474 | // |
---|
475 | // Differential PAI dEdx(omega)=omega*dNdx(omega) |
---|
476 | // |
---|
477 | |
---|
478 | G4double G4InitXscPAI::DifPAIdEdx( G4double omega ) |
---|
479 | { |
---|
480 | G4double dEdx = omega*DifPAIxSection(omega); |
---|
481 | return dEdx; |
---|
482 | } |
---|
483 | |
---|
484 | ////////////////////////////////////////////////////////////////////////// |
---|
485 | // |
---|
486 | // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons |
---|
487 | |
---|
488 | G4double G4InitXscPAI::PAIdNdxCherenkov( G4double omega ) |
---|
489 | { |
---|
490 | G4int i = fCurrentInterval; |
---|
491 | G4double betaGammaSq = fBetaGammaSq; |
---|
492 | G4double epsilonRe = RePartDielectricConst(omega); |
---|
493 | G4double epsilonIm = ImPartDielectricConst(i,omega); |
---|
494 | |
---|
495 | G4double cof, logarithm, x3, x5, argument, modul2, dNdxC ; |
---|
496 | G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr ; |
---|
497 | |
---|
498 | cof = 1.0 ; |
---|
499 | cofBetaBohr = 4.0 ; |
---|
500 | betaBohr2 = fine_structure_const*fine_structure_const ; |
---|
501 | betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ; |
---|
502 | |
---|
503 | be2 = betaGammaSq/(1 + betaGammaSq) ; |
---|
504 | be4 = be2*be2 ; |
---|
505 | |
---|
506 | if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ; |
---|
507 | else |
---|
508 | { |
---|
509 | logarithm = -log( (1/betaGammaSq - epsilonRe)* |
---|
510 | (1/betaGammaSq - epsilonRe) + |
---|
511 | epsilonIm*epsilonIm )*0.5 ; |
---|
512 | logarithm += log(1+1.0/betaGammaSq) ; |
---|
513 | } |
---|
514 | |
---|
515 | if( epsilonIm == 0.0 || betaGammaSq < 0.01 ) |
---|
516 | { |
---|
517 | argument = 0.0 ; |
---|
518 | } |
---|
519 | else |
---|
520 | { |
---|
521 | x3 = -epsilonRe + 1.0/betaGammaSq ; |
---|
522 | x5 = -1.0 - epsilonRe + |
---|
523 | be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) + |
---|
524 | epsilonIm*epsilonIm) ; |
---|
525 | if( x3 == 0.0 ) argument = 0.5*pi; |
---|
526 | else argument = atan2(epsilonIm,x3) ; |
---|
527 | argument *= x5 ; |
---|
528 | } |
---|
529 | dNdxC = ( logarithm*epsilonIm + argument )/hbarc ; |
---|
530 | |
---|
531 | if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ; |
---|
532 | |
---|
533 | dNdxC *= fine_structure_const/be2/pi ; |
---|
534 | |
---|
535 | dNdxC *= (1-exp(-be4/betaBohr4)) ; |
---|
536 | |
---|
537 | if(fDensity >= fSolidDensity) |
---|
538 | { |
---|
539 | modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) + |
---|
540 | epsilonIm*epsilonIm; |
---|
541 | dNdxC /= modul2 ; |
---|
542 | } |
---|
543 | return dNdxC ; |
---|
544 | |
---|
545 | } // end of PAIdNdxCerenkov |
---|
546 | |
---|
547 | ////////////////////////////////////////////////////////////////////////// |
---|
548 | // |
---|
549 | // Calculation od dN/dx of collisions with creation of longitudinal EM |
---|
550 | // excitations (plasmons, delta-electrons) |
---|
551 | |
---|
552 | G4double G4InitXscPAI::PAIdNdxPlasmon( G4double omega ) |
---|
553 | { |
---|
554 | G4int i = fCurrentInterval; |
---|
555 | G4double betaGammaSq = fBetaGammaSq; |
---|
556 | G4double integralTerm = IntegralTerm(omega); |
---|
557 | G4double epsilonRe = RePartDielectricConst(omega); |
---|
558 | G4double epsilonIm = ImPartDielectricConst(i,omega); |
---|
559 | |
---|
560 | G4double cof, resonance, modul2, dNdxP ; |
---|
561 | G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr ; |
---|
562 | |
---|
563 | cof = 1 ; |
---|
564 | cofBetaBohr = 4.0 ; |
---|
565 | betaBohr2 = fine_structure_const*fine_structure_const ; |
---|
566 | betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ; |
---|
567 | |
---|
568 | be2 = betaGammaSq/(1 + betaGammaSq) ; |
---|
569 | be4 = be2*be2 ; |
---|
570 | |
---|
571 | resonance = log(2*electron_mass_c2*be2/omega) ; |
---|
572 | resonance *= epsilonIm/hbarc ; |
---|
573 | |
---|
574 | |
---|
575 | dNdxP = ( resonance + cof*integralTerm/omega/omega ) ; |
---|
576 | |
---|
577 | if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ; |
---|
578 | |
---|
579 | dNdxP *= fine_structure_const/be2/pi ; |
---|
580 | dNdxP *= (1-exp(-be4/betaBohr4)) ; |
---|
581 | |
---|
582 | if( fDensity >= fSolidDensity ) |
---|
583 | { |
---|
584 | modul2 = (1 + epsilonRe)*(1 + epsilonRe) + |
---|
585 | epsilonIm*epsilonIm; |
---|
586 | dNdxP /= modul2 ; |
---|
587 | } |
---|
588 | return dNdxP ; |
---|
589 | |
---|
590 | } // end of PAIdNdxPlasmon |
---|
591 | |
---|
592 | //////////////////////////////////////////////////////////////////////// |
---|
593 | // |
---|
594 | // Calculation of the PAI integral cross-section |
---|
595 | // = specific primary ionisation, 1/cm |
---|
596 | // |
---|
597 | |
---|
598 | void G4InitXscPAI::IntegralPAIxSection(G4double bg2, G4double Tmax) |
---|
599 | { |
---|
600 | G4int i,k,i1,i2; |
---|
601 | G4double energy1, energy2, result = 0.; |
---|
602 | |
---|
603 | fBetaGammaSq = bg2; |
---|
604 | fTmax = Tmax; |
---|
605 | |
---|
606 | if(fPAIxscVector) delete fPAIxscVector; |
---|
607 | |
---|
608 | fPAIxscVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); |
---|
609 | fPAIxscVector->PutValue(fPAIbin-1,result); |
---|
610 | |
---|
611 | for( i = fIntervalNumber - 1; i >= 0; i-- ) |
---|
612 | { |
---|
613 | if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
614 | } |
---|
615 | if (i < 0) i = 0; // Tmax should be more than |
---|
616 | // first ionisation potential |
---|
617 | fIntervalTmax = i; |
---|
618 | |
---|
619 | G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral; |
---|
620 | |
---|
621 | for( k = fPAIbin - 2; k >= 0; k-- ) |
---|
622 | { |
---|
623 | energy1 = fPAIxscVector->GetLowEdgeEnergy(k); |
---|
624 | energy2 = fPAIxscVector->GetLowEdgeEnergy(k+1); |
---|
625 | |
---|
626 | for( i = fIntervalTmax; i >= 0; i-- ) |
---|
627 | { |
---|
628 | if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
629 | } |
---|
630 | if(i < 0) i = 0; |
---|
631 | i2 = i; |
---|
632 | |
---|
633 | for( i = fIntervalTmax; i >= 0; i-- ) |
---|
634 | { |
---|
635 | if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
636 | } |
---|
637 | if(i < 0) i = 0; |
---|
638 | i1 = i; |
---|
639 | |
---|
640 | if( i1 == i2 ) |
---|
641 | { |
---|
642 | fCurrentInterval = i1; |
---|
643 | result += integral.Legendre10(this,&G4InitXscPAI::DifPAIxSection, |
---|
644 | energy1,energy2); |
---|
645 | fPAIxscVector->PutValue(k,result); |
---|
646 | } |
---|
647 | else |
---|
648 | { |
---|
649 | for( i = i2; i >= i1; i-- ) |
---|
650 | { |
---|
651 | fCurrentInterval = i; |
---|
652 | |
---|
653 | if( i==i2 ) result += integral.Legendre10(this, |
---|
654 | &G4InitXscPAI::DifPAIxSection, |
---|
655 | (*(*fMatSandiaMatrix)[i])[0] ,energy2); |
---|
656 | |
---|
657 | else if( i == i1 ) result += integral.Legendre10(this, |
---|
658 | &G4InitXscPAI::DifPAIxSection,energy1, |
---|
659 | (*(*fMatSandiaMatrix)[i+1])[0]); |
---|
660 | |
---|
661 | else result += integral.Legendre10(this, |
---|
662 | &G4InitXscPAI::DifPAIxSection, |
---|
663 | (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]); |
---|
664 | } |
---|
665 | fPAIxscVector->PutValue(k,result); |
---|
666 | } |
---|
667 | // G4cout<<k<<"\t"<<result<<G4endl; |
---|
668 | } |
---|
669 | return ; |
---|
670 | } |
---|
671 | |
---|
672 | |
---|
673 | //////////////////////////////////////////////////////////////////////// |
---|
674 | // |
---|
675 | // Calculation of the PAI integral dEdx |
---|
676 | // = mean energy loss per unit length, keV/cm |
---|
677 | // |
---|
678 | |
---|
679 | void G4InitXscPAI::IntegralPAIdEdx(G4double bg2, G4double Tmax) |
---|
680 | { |
---|
681 | G4int i,k,i1,i2; |
---|
682 | G4double energy1, energy2, result = 0.; |
---|
683 | |
---|
684 | fBetaGammaSq = bg2; |
---|
685 | fTmax = Tmax; |
---|
686 | |
---|
687 | if(fPAIdEdxVector) delete fPAIdEdxVector; |
---|
688 | |
---|
689 | fPAIdEdxVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); |
---|
690 | fPAIdEdxVector->PutValue(fPAIbin-1,result); |
---|
691 | |
---|
692 | for( i = fIntervalNumber - 1; i >= 0; i-- ) |
---|
693 | { |
---|
694 | if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
695 | } |
---|
696 | if (i < 0) i = 0; // Tmax should be more than |
---|
697 | // first ionisation potential |
---|
698 | fIntervalTmax = i; |
---|
699 | |
---|
700 | G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral; |
---|
701 | |
---|
702 | for( k = fPAIbin - 2; k >= 0; k-- ) |
---|
703 | { |
---|
704 | energy1 = fPAIdEdxVector->GetLowEdgeEnergy(k); |
---|
705 | energy2 = fPAIdEdxVector->GetLowEdgeEnergy(k+1); |
---|
706 | |
---|
707 | for( i = fIntervalTmax; i >= 0; i-- ) |
---|
708 | { |
---|
709 | if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
710 | } |
---|
711 | if(i < 0) i = 0; |
---|
712 | i2 = i; |
---|
713 | |
---|
714 | for( i = fIntervalTmax; i >= 0; i-- ) |
---|
715 | { |
---|
716 | if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
717 | } |
---|
718 | if(i < 0) i = 0; |
---|
719 | i1 = i; |
---|
720 | |
---|
721 | if( i1 == i2 ) |
---|
722 | { |
---|
723 | fCurrentInterval = i1; |
---|
724 | result += integral.Legendre10(this,&G4InitXscPAI::DifPAIdEdx, |
---|
725 | energy1,energy2); |
---|
726 | fPAIdEdxVector->PutValue(k,result); |
---|
727 | } |
---|
728 | else |
---|
729 | { |
---|
730 | for( i = i2; i >= i1; i-- ) |
---|
731 | { |
---|
732 | fCurrentInterval = i; |
---|
733 | |
---|
734 | if( i==i2 ) result += integral.Legendre10(this, |
---|
735 | &G4InitXscPAI::DifPAIdEdx, |
---|
736 | (*(*fMatSandiaMatrix)[i])[0] ,energy2); |
---|
737 | |
---|
738 | else if( i == i1 ) result += integral.Legendre10(this, |
---|
739 | &G4InitXscPAI::DifPAIdEdx,energy1, |
---|
740 | (*(*fMatSandiaMatrix)[i+1])[0]); |
---|
741 | |
---|
742 | else result += integral.Legendre10(this, |
---|
743 | &G4InitXscPAI::DifPAIdEdx, |
---|
744 | (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]); |
---|
745 | } |
---|
746 | fPAIdEdxVector->PutValue(k,result); |
---|
747 | } |
---|
748 | // G4cout<<k<<"\t"<<result<<G4endl; |
---|
749 | } |
---|
750 | return ; |
---|
751 | } |
---|
752 | |
---|
753 | //////////////////////////////////////////////////////////////////////// |
---|
754 | // |
---|
755 | // Calculation of the PAI Cerenkov integral cross-section |
---|
756 | // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm |
---|
757 | // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm |
---|
758 | |
---|
759 | void G4InitXscPAI::IntegralCherenkov(G4double bg2, G4double Tmax) |
---|
760 | { |
---|
761 | G4int i,k,i1,i2; |
---|
762 | G4double energy1, energy2, beta2, module2, cos2, width, result = 0.; |
---|
763 | |
---|
764 | fBetaGammaSq = bg2; |
---|
765 | fTmax = Tmax; |
---|
766 | beta2 = bg2/(1+bg2); |
---|
767 | |
---|
768 | if(fPAIphotonVector) delete fPAIphotonVector; |
---|
769 | if(fChCosSqVector) delete fChCosSqVector; |
---|
770 | if(fChWidthVector) delete fChWidthVector; |
---|
771 | |
---|
772 | fPAIphotonVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); |
---|
773 | fChCosSqVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); |
---|
774 | fChWidthVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); |
---|
775 | |
---|
776 | fPAIphotonVector->PutValue(fPAIbin-1,result); |
---|
777 | fChCosSqVector->PutValue(fPAIbin-1,1.); |
---|
778 | fChWidthVector->PutValue(fPAIbin-1,1e-7); |
---|
779 | |
---|
780 | for( i = fIntervalNumber - 1; i >= 0; i-- ) |
---|
781 | { |
---|
782 | if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
783 | } |
---|
784 | if (i < 0) i = 0; // Tmax should be more than |
---|
785 | // first ionisation potential |
---|
786 | fIntervalTmax = i; |
---|
787 | |
---|
788 | G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral; |
---|
789 | |
---|
790 | for( k = fPAIbin - 2; k >= 0; k-- ) |
---|
791 | { |
---|
792 | energy1 = fPAIphotonVector->GetLowEdgeEnergy(k); |
---|
793 | energy2 = fPAIphotonVector->GetLowEdgeEnergy(k+1); |
---|
794 | |
---|
795 | for( i = fIntervalTmax; i >= 0; i-- ) |
---|
796 | { |
---|
797 | if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
798 | } |
---|
799 | if(i < 0) i = 0; |
---|
800 | i2 = i; |
---|
801 | |
---|
802 | for( i = fIntervalTmax; i >= 0; i-- ) |
---|
803 | { |
---|
804 | if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
805 | } |
---|
806 | if(i < 0) i = 0; |
---|
807 | i1 = i; |
---|
808 | |
---|
809 | module2 = ModuleSqDielectricConst(i1,energy1); |
---|
810 | cos2 = RePartDielectricConst(energy1)/module2/beta2; |
---|
811 | width = ImPartDielectricConst(i1,energy1)/module2/beta2; |
---|
812 | |
---|
813 | fChCosSqVector->PutValue(k,cos2); |
---|
814 | fChWidthVector->PutValue(k,width); |
---|
815 | |
---|
816 | if( i1 == i2 ) |
---|
817 | { |
---|
818 | fCurrentInterval = i1; |
---|
819 | result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxCherenkov, |
---|
820 | energy1,energy2); |
---|
821 | fPAIphotonVector->PutValue(k,result); |
---|
822 | |
---|
823 | } |
---|
824 | else |
---|
825 | { |
---|
826 | for( i = i2; i >= i1; i-- ) |
---|
827 | { |
---|
828 | fCurrentInterval = i; |
---|
829 | |
---|
830 | if( i==i2 ) result += integral.Legendre10(this, |
---|
831 | &G4InitXscPAI::PAIdNdxCherenkov, |
---|
832 | (*(*fMatSandiaMatrix)[i])[0] ,energy2); |
---|
833 | |
---|
834 | else if( i == i1 ) result += integral.Legendre10(this, |
---|
835 | &G4InitXscPAI::PAIdNdxCherenkov,energy1, |
---|
836 | (*(*fMatSandiaMatrix)[i+1])[0]); |
---|
837 | |
---|
838 | else result += integral.Legendre10(this, |
---|
839 | &G4InitXscPAI::PAIdNdxCherenkov, |
---|
840 | (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]); |
---|
841 | } |
---|
842 | fPAIphotonVector->PutValue(k,result); |
---|
843 | } |
---|
844 | // G4cout<<k<<"\t"<<result<<G4endl; |
---|
845 | } |
---|
846 | return; |
---|
847 | } // end of IntegralCerenkov |
---|
848 | |
---|
849 | //////////////////////////////////////////////////////////////////////// |
---|
850 | // |
---|
851 | // Calculation of the PAI Plasmon integral cross-section |
---|
852 | // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm |
---|
853 | // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm |
---|
854 | |
---|
855 | void G4InitXscPAI::IntegralPlasmon(G4double bg2, G4double Tmax) |
---|
856 | { |
---|
857 | G4int i,k,i1,i2; |
---|
858 | G4double energy1, energy2, result = 0.; |
---|
859 | |
---|
860 | fBetaGammaSq = bg2; |
---|
861 | fTmax = Tmax; |
---|
862 | |
---|
863 | if(fPAIelectronVector) delete fPAIelectronVector; |
---|
864 | |
---|
865 | fPAIelectronVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); |
---|
866 | fPAIelectronVector->PutValue(fPAIbin-1,result); |
---|
867 | |
---|
868 | for( i = fIntervalNumber - 1; i >= 0; i-- ) |
---|
869 | { |
---|
870 | if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
871 | } |
---|
872 | if (i < 0) i = 0; // Tmax should be more than |
---|
873 | // first ionisation potential |
---|
874 | fIntervalTmax = i; |
---|
875 | |
---|
876 | G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral; |
---|
877 | |
---|
878 | for( k = fPAIbin - 2; k >= 0; k-- ) |
---|
879 | { |
---|
880 | energy1 = fPAIelectronVector->GetLowEdgeEnergy(k); |
---|
881 | energy2 = fPAIelectronVector->GetLowEdgeEnergy(k+1); |
---|
882 | |
---|
883 | for( i = fIntervalTmax; i >= 0; i-- ) |
---|
884 | { |
---|
885 | if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
886 | } |
---|
887 | if(i < 0) i = 0; |
---|
888 | i2 = i; |
---|
889 | |
---|
890 | for( i = fIntervalTmax; i >= 0; i-- ) |
---|
891 | { |
---|
892 | if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break; |
---|
893 | } |
---|
894 | if(i < 0) i = 0; |
---|
895 | i1 = i; |
---|
896 | |
---|
897 | if( i1 == i2 ) |
---|
898 | { |
---|
899 | fCurrentInterval = i1; |
---|
900 | result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxPlasmon, |
---|
901 | energy1,energy2); |
---|
902 | fPAIelectronVector->PutValue(k,result); |
---|
903 | } |
---|
904 | else |
---|
905 | { |
---|
906 | for( i = i2; i >= i1; i-- ) |
---|
907 | { |
---|
908 | fCurrentInterval = i; |
---|
909 | |
---|
910 | if( i==i2 ) result += integral.Legendre10(this, |
---|
911 | &G4InitXscPAI::PAIdNdxPlasmon, |
---|
912 | (*(*fMatSandiaMatrix)[i])[0] ,energy2); |
---|
913 | |
---|
914 | else if( i == i1 ) result += integral.Legendre10(this, |
---|
915 | &G4InitXscPAI::PAIdNdxPlasmon,energy1, |
---|
916 | (*(*fMatSandiaMatrix)[i+1])[0]); |
---|
917 | |
---|
918 | else result += integral.Legendre10(this, |
---|
919 | &G4InitXscPAI::PAIdNdxPlasmon, |
---|
920 | (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]); |
---|
921 | } |
---|
922 | fPAIelectronVector->PutValue(k,result); |
---|
923 | } |
---|
924 | // G4cout<<k<<"\t"<<result<<G4endl; |
---|
925 | } |
---|
926 | return; |
---|
927 | } // end of IntegralPlasmon |
---|
928 | |
---|
929 | |
---|
930 | ///////////////////////////////////////////////////////////////////////// |
---|
931 | // |
---|
932 | // |
---|
933 | |
---|
934 | G4double G4InitXscPAI::GetPhotonLambda( G4double omega ) |
---|
935 | { |
---|
936 | G4int i ; |
---|
937 | G4double omega2, omega3, omega4, a1, a2, a3, a4, lambda ; |
---|
938 | |
---|
939 | omega2 = omega*omega ; |
---|
940 | omega3 = omega2*omega ; |
---|
941 | omega4 = omega2*omega2 ; |
---|
942 | |
---|
943 | for(i = 0; i < fIntervalNumber;i++) |
---|
944 | { |
---|
945 | if( omega < (*(*fMatSandiaMatrix)[i])[0] ) break ; |
---|
946 | } |
---|
947 | if( i == 0 ) |
---|
948 | { |
---|
949 | G4cout<<"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<G4endl; |
---|
950 | } |
---|
951 | else i-- ; |
---|
952 | |
---|
953 | a1 = (*(*fMatSandiaMatrix)[i])[1]; |
---|
954 | a2 = (*(*fMatSandiaMatrix)[i])[2]; |
---|
955 | a3 = (*(*fMatSandiaMatrix)[i])[3]; |
---|
956 | a4 = (*(*fMatSandiaMatrix)[i])[4]; |
---|
957 | |
---|
958 | lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4); |
---|
959 | |
---|
960 | return lambda ; |
---|
961 | } |
---|
962 | |
---|
963 | ///////////////////////////////////////////////////////////////////////// |
---|
964 | // |
---|
965 | // |
---|
966 | |
---|
967 | ///////////////////////////////////////////////////////////////////////// |
---|
968 | // |
---|
969 | // |
---|
970 | |
---|
971 | G4double G4InitXscPAI::GetStepEnergyLoss( G4double step ) |
---|
972 | { |
---|
973 | G4double loss = 0.0 ; |
---|
974 | loss *= step; |
---|
975 | |
---|
976 | return loss ; |
---|
977 | } |
---|
978 | |
---|
979 | ///////////////////////////////////////////////////////////////////////// |
---|
980 | // |
---|
981 | // |
---|
982 | |
---|
983 | G4double G4InitXscPAI::GetStepCerenkovLoss( G4double step ) |
---|
984 | { |
---|
985 | G4double loss = 0.0 ; |
---|
986 | loss *= step; |
---|
987 | |
---|
988 | return loss ; |
---|
989 | } |
---|
990 | |
---|
991 | ///////////////////////////////////////////////////////////////////////// |
---|
992 | // |
---|
993 | // |
---|
994 | |
---|
995 | G4double G4InitXscPAI::GetStepPlasmonLoss( G4double step ) |
---|
996 | { |
---|
997 | |
---|
998 | |
---|
999 | G4double loss = 0.0 ; |
---|
1000 | loss *= step; |
---|
1001 | return loss ; |
---|
1002 | } |
---|
1003 | |
---|
1004 | |
---|
1005 | // |
---|
1006 | // end of G4InitXscPAI implementation file |
---|
1007 | // |
---|
1008 | //////////////////////////////////////////////////////////////////////////// |
---|
1009 | |
---|