source: trunk/source/processes/electromagnetic/standard/src/G4PAIySection.cc@ 1344

Last change on this file since 1344 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 36.3 KB
Line 
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// $Id: G4PAIySection.cc,v 1.6 2010/11/04 17:30:32 vnivanch Exp $
27// GEANT4 tag $Name: emstand-V09-03-25 $
28//
29//
30// G4PAIySection.cc -- class implementation file
31//
32// GEANT 4 class implementation file
33//
34// For information related to this code, please, contact
35// the Geant4 Collaboration.
36//
37// R&D: Vladimir.Grichine@cern.ch
38//
39// History:
40//
41// 01.10.07, V.Ivanchenko create using V.Grichine G4PAIxSection class
42// 26.07.09, V.Ivanchenko added protection for mumerical exceptions for
43// low-density materials
44//
45
46#include "G4PAIySection.hh"
47
48#include "globals.hh"
49#include "G4ios.hh"
50#include "G4Poisson.hh"
51#include "G4Material.hh"
52#include "G4MaterialCutsCouple.hh"
53#include "G4SandiaTable.hh"
54
55using namespace std;
56
57// Local class constants
58
59const G4double G4PAIySection::fDelta = 0.005 ; // energy shift from interval border
60const G4double G4PAIySection::fError = 0.005 ; // error in lin-log approximation
61
62const G4int G4PAIySection::fMaxSplineSize = 500 ; // Max size of output spline
63 // arrays
64
65//////////////////////////////////////////////////////////////////
66//
67// Constructor
68//
69
70G4PAIySection::G4PAIySection()
71{
72 fSandia = 0;
73 fDensity = fElectronDensity = fNormalizationCof = 0.0;
74 fIntervalNumber = fSplineNumber = 0;
75}
76
77////////////////////////////////////////////////////////////////////////////
78//
79// Destructor
80
81G4PAIySection::~G4PAIySection()
82{}
83
84////////////////////////////////////////////////////////////////////////
85//
86// Test Constructor with beta*gamma square value
87
88void G4PAIySection::Initialize( const G4Material* material,
89 G4double maxEnergyTransfer,
90 G4double betaGammaSq)
91{
92 G4int i, j, numberOfElements ;
93
94 fDensity = material->GetDensity();
95 fElectronDensity = material->GetElectronDensity() ;
96 numberOfElements = material->GetNumberOfElements() ;
97
98 fSandia = material->GetSandiaTable();
99
100 fIntervalNumber = fSandia->GetMaxInterval();
101
102 fIntervalNumber--;
103
104 for(i=1;i<=fIntervalNumber;i++)
105 {
106 G4double e = fSandia->GetSandiaMatTablePAI(i,0);
107 if(e >= maxEnergyTransfer || i > fIntervalNumber)
108 {
109 fEnergyInterval[i] = maxEnergyTransfer ;
110 fIntervalNumber = i ;
111 break;
112 }
113 fEnergyInterval[i] = e;
114 fA1[i] = fSandia->GetSandiaMatTablePAI(i,1);
115 fA2[i] = fSandia->GetSandiaMatTablePAI(i,2);
116 fA3[i] = fSandia->GetSandiaMatTablePAI(i,3);
117 fA4[i] = fSandia->GetSandiaMatTablePAI(i,4);
118
119 }
120 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
121 {
122 fIntervalNumber++;
123 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer ;
124 fA1[fIntervalNumber] = fA1[fIntervalNumber-1] ;
125 fA2[fIntervalNumber] = fA2[fIntervalNumber-1] ;
126 fA3[fIntervalNumber] = fA3[fIntervalNumber-1] ;
127 fA4[fIntervalNumber] = fA4[fIntervalNumber-1] ;
128 }
129
130 // Now checking, if two borders are too close together
131 for(i=1;i<fIntervalNumber;i++)
132 {
133 // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
134 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
135 if(fEnergyInterval[i+1]-fEnergyInterval[i] <
136 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
137 {
138 for(j=i;j<fIntervalNumber;j++)
139 {
140 fEnergyInterval[j] = fEnergyInterval[j+1] ;
141 fA1[j] = fA1[j+1] ;
142 fA2[j] = fA2[j+1] ;
143 fA3[j] = fA3[j+1] ;
144 fA4[j] = fA4[j+1] ;
145 }
146 fIntervalNumber-- ;
147 i-- ;
148 }
149 }
150
151 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
152
153 G4double betaGammaSqRef =
154 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
155
156 NormShift(betaGammaSqRef) ;
157 SplainPAI(betaGammaSqRef) ;
158
159 // Preparation of integral PAI cross section for input betaGammaSq
160
161 for(i = 1 ; i <= fSplineNumber ; i++)
162 {
163 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
164 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
165 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
166 }
167 IntegralPAIySection() ;
168 IntegralCerenkov() ;
169 IntegralPlasmon() ;
170}
171
172/////////////////////////////////////////////////////////////////////////
173//
174// General control function for class G4PAIySection
175//
176
177void G4PAIySection::InitPAI()
178{
179 G4int i ;
180 G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
181 fLorentzFactor[fRefGammaNumber] - 1;
182
183 // Preparation of integral PAI cross section for reference gamma
184
185 NormShift(betaGammaSq) ;
186 SplainPAI(betaGammaSq) ;
187
188 IntegralPAIySection() ;
189 IntegralCerenkov() ;
190 IntegralPlasmon() ;
191
192 for(i = 0 ; i<=fSplineNumber ; i++)
193 {
194 fPAItable[i][fRefGammaNumber] = fIntegralPAIySection[i] ;
195 if(i != 0)
196 {
197 fPAItable[i][0] = fSplineEnergy[i] ;
198 }
199 }
200 fPAItable[0][0] = fSplineNumber ;
201
202 for(G4int j = 1 ; j < 112 ; j++) // for other gammas
203 {
204 if( j == fRefGammaNumber ) continue ;
205
206 betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1 ;
207
208 for(i = 1 ; i <= fSplineNumber ; i++)
209 {
210 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
211 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
212 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
213 }
214 IntegralPAIySection() ;
215 IntegralCerenkov() ;
216 IntegralPlasmon() ;
217
218 for(i = 0 ; i <= fSplineNumber ; i++)
219 {
220 fPAItable[i][j] = fIntegralPAIySection[i] ;
221 }
222 }
223
224}
225
226///////////////////////////////////////////////////////////////////////
227//
228// Shifting from borders to intervals Creation of first energy points
229//
230
231void G4PAIySection::NormShift(G4double betaGammaSq)
232{
233 G4int i, j ;
234
235 for( i = 1 ; i <= fIntervalNumber-1 ; i++ )
236 {
237 for( j = 1 ; j <= 2 ; j++ )
238 {
239 fSplineNumber = (i-1)*2 + j ;
240
241 if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
242 else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
243 // G4cout<<"cn = "<<fSplineNumber<<"; "<<"energy = "
244 // <<fSplineEnergy[fSplineNumber]<<G4endl;
245 }
246 }
247 fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
248
249 j = 1 ;
250
251 for(i=2;i<=fSplineNumber;i++)
252 {
253 if(fSplineEnergy[i]<fEnergyInterval[j+1])
254 {
255 fIntegralTerm[i] = fIntegralTerm[i-1] +
256 RutherfordIntegral(j,fSplineEnergy[i-1],
257 fSplineEnergy[i] ) ;
258 }
259 else
260 {
261 G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
262 fEnergyInterval[j+1] ) ;
263 j++;
264 fIntegralTerm[i] = fIntegralTerm[i-1] + x +
265 RutherfordIntegral(j,fEnergyInterval[j],
266 fSplineEnergy[i] ) ;
267 }
268 // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t"<<fIntegralTerm[i]<<"\n"<<G4endl;
269 }
270 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
271 fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber] ;
272
273 // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl ;
274
275 // Calculation of PAI differrential cross-section (1/(keV*cm))
276 // in the energy points near borders of energy intervals
277
278 for(G4int k=1;k<=fIntervalNumber-1;k++)
279 {
280 for(j=1;j<=2;j++)
281 {
282 i = (k-1)*2 + j ;
283 fImPartDielectricConst[i] = fNormalizationCof*
284 ImPartDielectricConst(k,fSplineEnergy[i]);
285 fRePartDielectricConst[i] = fNormalizationCof*
286 RePartDielectricConst(fSplineEnergy[i]);
287 fIntegralTerm[i] *= fNormalizationCof;
288
289 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
290 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
291 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
292 }
293 }
294
295} // end of NormShift
296
297/////////////////////////////////////////////////////////////////////////
298//
299// Creation of new energy points as geometrical mean of existing
300// one, calculation PAI_cs for them, while the error of logarithmic
301// linear approximation would be smaller than 'fError'
302
303void G4PAIySection::SplainPAI(G4double betaGammaSq)
304{
305 G4int k = 1 ;
306 G4int i = 1 ;
307
308 while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
309 {
310 if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
311 {
312 k++ ; // Here next energy point is in next energy interval
313 i++;
314 continue;
315 }
316 // Shifting of arrayes for inserting the geometrical
317 // average of 'i' and 'i+1' energy points to 'i+1' place
318 fSplineNumber++;
319
320 for(G4int j = fSplineNumber; j >= i+2 ; j-- )
321 {
322 fSplineEnergy[j] = fSplineEnergy[j-1];
323 fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
324 fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
325 fIntegralTerm[j] = fIntegralTerm[j-1];
326
327 fDifPAIySection[j] = fDifPAIySection[j-1];
328 fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
329 fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
330 }
331 G4double x1 = fSplineEnergy[i];
332 G4double x2 = fSplineEnergy[i+1];
333 G4double yy1 = fDifPAIySection[i];
334 G4double y2 = fDifPAIySection[i+1];
335
336 G4double en1 = sqrt(x1*x2);
337 fSplineEnergy[i+1] = en1;
338
339 // Calculation of logarithmic linear approximation
340 // in this (enr) energy point, which number is 'i+1' now
341
342 G4double a = log10(y2/yy1)/log10(x2/x1);
343 G4double b = log10(yy1) - a*log10(x1);
344 G4double y = a*log10(en1) + b ;
345 y = pow(10.,y);
346
347 // Calculation of the PAI dif. cross-section at this point
348
349 fImPartDielectricConst[i+1] = fNormalizationCof*
350 ImPartDielectricConst(k,fSplineEnergy[i+1]);
351 fRePartDielectricConst[i+1] = fNormalizationCof*
352 RePartDielectricConst(fSplineEnergy[i+1]);
353 fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
354 RutherfordIntegral(k,fSplineEnergy[i],
355 fSplineEnergy[i+1]);
356
357 fDifPAIySection[i+1] = DifPAIySection(i+1,betaGammaSq);
358 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
359 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
360
361 // Condition for next division of this segment or to pass
362 // to higher energies
363
364 G4double x = 2*(fDifPAIySection[i+1] - y)/(fDifPAIySection[i+1] + y);
365
366 if( x < 0 )
367 {
368 x = -x ;
369 }
370 if( x > fError && fSplineNumber < fMaxSplineSize-1 )
371 {
372 continue; // next division
373 }
374 i += 2; // pass to next segment
375
376 } // close 'while'
377
378} // end of SplainPAI
379
380
381////////////////////////////////////////////////////////////////////
382//
383// Integration over electrons that could be considered
384// quasi-free at energy transfer of interest
385
386G4double G4PAIySection::RutherfordIntegral( G4int k,
387 G4double x1,
388 G4double x2 )
389{
390 G4double c1, c2, c3 ;
391 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
392 c1 = (x2 - x1)/x1/x2 ;
393 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
394 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
395 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
396
397 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3 ;
398
399} // end of RutherfordIntegral
400
401
402/////////////////////////////////////////////////////////////////
403//
404// Imaginary part of dielectric constant
405// (G4int k - interval number, G4double en1 - energy point)
406
407G4double G4PAIySection::ImPartDielectricConst( G4int k ,
408 G4double energy1 )
409{
410 G4double energy2,energy3,energy4,result;
411
412 energy2 = energy1*energy1;
413 energy3 = energy2*energy1;
414 energy4 = energy3*energy1;
415
416 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4 ;
417 result *=hbarc/energy1 ;
418
419 return result ;
420
421} // end of ImPartDielectricConst
422
423
424//////////////////////////////////////////////////////////////////////////////
425//
426// Real part of dielectric constant minus unit: epsilon_1 - 1
427// (G4double enb - energy point)
428//
429
430G4double G4PAIySection::RePartDielectricConst(G4double enb)
431{
432 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
433 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
434
435 x0 = enb ;
436 result = 0 ;
437
438 for(G4int i=1;i<=fIntervalNumber-1;i++)
439 {
440 x1 = fEnergyInterval[i] ;
441 x2 = fEnergyInterval[i+1] ;
442 xx1 = x1 - x0 ;
443 xx2 = x2 - x0 ;
444 xx12 = xx2/xx1 ;
445
446 if(xx12<0)
447 {
448 xx12 = -xx12;
449 }
450 xln1 = log(x2/x1) ;
451 xln2 = log(xx12) ;
452 xln3 = log((x2 + x0)/(x1 + x0)) ;
453 x02 = x0*x0 ;
454 x03 = x02*x0 ;
455 x04 = x03*x0 ;
456 x05 = x04*x0;
457 c1 = (x2 - x1)/x1/x2 ;
458 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
459 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
460
461 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1 ;
462 result -= (fA2[i]/x02 + fA4[i]/x04)*c1 ;
463 result -= fA3[i]*c2/2/x02 ;
464 result -= fA4[i]*c3/3/x02 ;
465
466 cof1 = fA1[i]/x02 + fA3[i]/x04 ;
467 cof2 = fA2[i]/x03 + fA4[i]/x05 ;
468
469 result += 0.5*(cof1 +cof2)*xln2 ;
470 result += 0.5*(cof1 - cof2)*xln3 ;
471 }
472 result *= 2*hbarc/pi ;
473
474 return result ;
475
476} // end of RePartDielectricConst
477
478//////////////////////////////////////////////////////////////////////
479//
480// PAI differential cross-section in terms of
481// simplified Allison's equation
482//
483
484G4double G4PAIySection::DifPAIySection( G4int i ,
485 G4double betaGammaSq )
486{
487 G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
488 //G4double beta, be4 ;
489 G4double be4 ;
490 G4double betaBohr2 = fine_structure_const*fine_structure_const ;
491 G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
492 be2 = betaGammaSq/(1 + betaGammaSq) ;
493 be4 = be2*be2 ;
494 // beta = sqrt(be2) ;
495 cof = 1 ;
496 x1 = log(2*electron_mass_c2/fSplineEnergy[i]) ;
497
498 if( betaGammaSq < 0.01 ) x2 = log(be2) ;
499 else
500 {
501 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
502 (1/betaGammaSq - fRePartDielectricConst[i]) +
503 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2 ;
504 }
505 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
506 {
507 x6=0 ;
508 }
509 else
510 {
511 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq ;
512 x5 = -1 - fRePartDielectricConst[i] +
513 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
514 fImPartDielectricConst[i]*fImPartDielectricConst[i]) ;
515
516 x7 = atan2(fImPartDielectricConst[i],x3) ;
517 x6 = x5 * x7 ;
518 }
519 // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
520
521 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc ;
522 // if( x4 < 0.0 ) x4 = 0.0 ;
523 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
524 fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
525
526 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]) ;
527 if(result < 1.0e-8) result = 1.0e-8 ;
528 result *= fine_structure_const/be2/pi ;
529 // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
530 // result *= (1-exp(-be2/betaBohr2)) ;
531 result *= (1-exp(-be4/betaBohr4)) ;
532 // if(fDensity >= 0.1)
533 if(x8 > 0.)
534 {
535 result /= x8 ;
536 }
537 return result ;
538
539} // end of DifPAIySection
540
541//////////////////////////////////////////////////////////////////////////
542//
543// Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
544
545G4double G4PAIySection::PAIdNdxCerenkov( G4int i ,
546 G4double betaGammaSq )
547{
548 G4double cof, logarithm, x3, x5, argument, modul2, dNdxC ;
549 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr ;
550
551 cof = 1.0 ;
552 cofBetaBohr = 4.0 ;
553 betaBohr2 = fine_structure_const*fine_structure_const ;
554 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
555
556 be2 = betaGammaSq/(1 + betaGammaSq) ;
557 be4 = be2*be2 ;
558
559 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
560 else
561 {
562 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
563 (1/betaGammaSq - fRePartDielectricConst[i]) +
564 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5 ;
565 logarithm += log(1+1.0/betaGammaSq) ;
566 }
567
568 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
569 {
570 argument = 0.0 ;
571 }
572 else
573 {
574 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq ;
575 x5 = -1.0 - fRePartDielectricConst[i] +
576 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
577 fImPartDielectricConst[i]*fImPartDielectricConst[i]) ;
578 if( x3 == 0.0 ) argument = 0.5*pi;
579 else argument = atan2(fImPartDielectricConst[i],x3) ;
580 argument *= x5 ;
581 }
582 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc ;
583
584 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
585
586 dNdxC *= fine_structure_const/be2/pi ;
587
588 dNdxC *= (1-exp(-be4/betaBohr4)) ;
589
590 // if(fDensity >= 0.1)
591 // {
592 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
593 fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
594 if(modul2 > 0.)
595 {
596 dNdxC /= modul2 ;
597 }
598 return dNdxC ;
599
600} // end of PAIdNdxCerenkov
601
602//////////////////////////////////////////////////////////////////////////
603//
604// Calculation od dN/dx of collisions with creation of longitudinal EM
605// excitations (plasmons, delta-electrons)
606
607G4double G4PAIySection::PAIdNdxPlasmon( G4int i ,
608 G4double betaGammaSq )
609{
610 G4double cof, resonance, modul2, dNdxP ;
611 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr ;
612
613 cof = 1 ;
614 cofBetaBohr = 4.0 ;
615 betaBohr2 = fine_structure_const*fine_structure_const ;
616 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
617
618 be2 = betaGammaSq/(1 + betaGammaSq) ;
619 be4 = be2*be2 ;
620
621 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]) ;
622 resonance *= fImPartDielectricConst[i]/hbarc ;
623
624
625 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] ) ;
626
627 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
628
629 dNdxP *= fine_structure_const/be2/pi ;
630 dNdxP *= (1-exp(-be4/betaBohr4)) ;
631
632// if( fDensity >= 0.1 )
633// {
634 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
635 fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
636 if(modul2 > 0.)
637 {
638 dNdxP /= modul2 ;
639 }
640 return dNdxP ;
641
642} // end of PAIdNdxPlasmon
643
644////////////////////////////////////////////////////////////////////////
645//
646// Calculation of the PAI integral cross-section
647// fIntegralPAIySection[1] = specific primary ionisation, 1/cm
648// and fIntegralPAIySection[0] = mean energy loss per cm in keV/cm
649
650void G4PAIySection::IntegralPAIySection()
651{
652 fIntegralPAIySection[fSplineNumber] = 0 ;
653 fIntegralPAIdEdx[fSplineNumber] = 0 ;
654 fIntegralPAIySection[0] = 0 ;
655 G4int k = fIntervalNumber -1 ;
656
657 for(G4int i = fSplineNumber-1 ; i >= 1 ; i--)
658 {
659 if(fSplineEnergy[i] >= fEnergyInterval[k])
660 {
661 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + SumOverInterval(i) ;
662 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i) ;
663 }
664 else
665 {
666 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] +
667 SumOverBorder(i+1,fEnergyInterval[k]) ;
668 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
669 SumOverBorderdEdx(i+1,fEnergyInterval[k]) ;
670 k-- ;
671 }
672 }
673} // end of IntegralPAIySection
674
675////////////////////////////////////////////////////////////////////////
676//
677// Calculation of the PAI Cerenkov integral cross-section
678// fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
679// and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
680
681void G4PAIySection::IntegralCerenkov()
682{
683 G4int i, k ;
684 fIntegralCerenkov[fSplineNumber] = 0 ;
685 fIntegralCerenkov[0] = 0 ;
686 k = fIntervalNumber -1 ;
687
688 for( i = fSplineNumber-1 ; i >= 1 ; i-- )
689 {
690 if(fSplineEnergy[i] >= fEnergyInterval[k])
691 {
692 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i) ;
693 // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
694 }
695 else
696 {
697 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
698 SumOverBordCerenkov(i+1,fEnergyInterval[k]) ;
699 k-- ;
700 // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
701 }
702 }
703
704} // end of IntegralCerenkov
705
706////////////////////////////////////////////////////////////////////////
707//
708// Calculation of the PAI Plasmon integral cross-section
709// fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
710// and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
711
712void G4PAIySection::IntegralPlasmon()
713{
714 fIntegralPlasmon[fSplineNumber] = 0 ;
715 fIntegralPlasmon[0] = 0 ;
716 G4int k = fIntervalNumber -1 ;
717 for(G4int i=fSplineNumber-1;i>=1;i--)
718 {
719 if(fSplineEnergy[i] >= fEnergyInterval[k])
720 {
721 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i) ;
722 }
723 else
724 {
725 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
726 SumOverBordPlasmon(i+1,fEnergyInterval[k]) ;
727 k-- ;
728 }
729 }
730
731} // end of IntegralPlasmon
732
733//////////////////////////////////////////////////////////////////////
734//
735// Calculation the PAI integral cross-section inside
736// of interval of continuous values of photo-ionisation
737// cross-section. Parameter 'i' is the number of interval.
738
739G4double G4PAIySection::SumOverInterval( G4int i )
740{
741 G4double x0,x1,y0,yy1,a,b,c,result ;
742
743 x0 = fSplineEnergy[i] ;
744 x1 = fSplineEnergy[i+1] ;
745 y0 = fDifPAIySection[i] ;
746 yy1 = fDifPAIySection[i+1];
747 c = x1/x0;
748 a = log10(yy1/y0)/log10(c) ;
749 // b = log10(y0) - a*log10(x0) ;
750 b = y0/pow(x0,a) ;
751 a += 1 ;
752 if(a == 0)
753 {
754 result = b*log(x1/x0) ;
755 }
756 else
757 {
758 result = y0*(x1*pow(c,a-1) - x0)/a ;
759 }
760 a++;
761 if(a == 0)
762 {
763 fIntegralPAIySection[0] += b*log(x1/x0) ;
764 }
765 else
766 {
767 fIntegralPAIySection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
768 }
769 return result ;
770
771} // end of SumOverInterval
772
773/////////////////////////////////
774
775G4double G4PAIySection::SumOverIntervaldEdx( G4int i )
776{
777 G4double x0,x1,y0,yy1,a,b,c,result ;
778
779 x0 = fSplineEnergy[i] ;
780 x1 = fSplineEnergy[i+1] ;
781 y0 = fDifPAIySection[i] ;
782 yy1 = fDifPAIySection[i+1];
783 c = x1/x0;
784 a = log10(yy1/y0)/log10(c) ;
785 // b = log10(y0) - a*log10(x0) ;
786 b = y0/pow(x0,a) ;
787 a += 2 ;
788 if(a == 0)
789 {
790 result = b*log(x1/x0) ;
791 }
792 else
793 {
794 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
795 }
796 return result ;
797
798} // end of SumOverInterval
799
800//////////////////////////////////////////////////////////////////////
801//
802// Calculation the PAI Cerenkov integral cross-section inside
803// of interval of continuous values of photo-ionisation Cerenkov
804// cross-section. Parameter 'i' is the number of interval.
805
806G4double G4PAIySection::SumOverInterCerenkov( G4int i )
807{
808 G4double x0,x1,y0,yy1,a,c,result ;
809
810 x0 = fSplineEnergy[i] ;
811 x1 = fSplineEnergy[i+1] ;
812 y0 = fdNdxCerenkov[i] ;
813 yy1 = fdNdxCerenkov[i+1];
814 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
815 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
816
817 c = x1/x0;
818 a = log10(yy1/y0)/log10(c) ;
819 G4double b = 0.0;
820 if(a < 20.) b = y0/pow(x0,a) ;
821
822 a += 1.0 ;
823 if(a == 0) result = b*log(c) ;
824 else result = y0*(x1*pow(c,a-1) - x0)/a ;
825 a += 1.0 ;
826
827 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0) ;
828 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
829 // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
830 return result ;
831
832} // end of SumOverInterCerenkov
833
834//////////////////////////////////////////////////////////////////////
835//
836// Calculation the PAI Plasmon integral cross-section inside
837// of interval of continuous values of photo-ionisation Plasmon
838// cross-section. Parameter 'i' is the number of interval.
839
840G4double G4PAIySection::SumOverInterPlasmon( G4int i )
841{
842 G4double x0,x1,y0,yy1,a,c,result ;
843
844 x0 = fSplineEnergy[i] ;
845 x1 = fSplineEnergy[i+1] ;
846 y0 = fdNdxPlasmon[i] ;
847 yy1 = fdNdxPlasmon[i+1];
848 c =x1/x0;
849 a = log10(yy1/y0)/log10(c) ;
850
851 G4double b = 0.0;
852 if(a < 20.) b = y0/pow(x0,a) ;
853
854 a += 1.0 ;
855 if(a == 0) result = b*log(x1/x0) ;
856 else result = y0*(x1*pow(c,a-1) - x0)/a ;
857 a += 1.0 ;
858
859 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0) ;
860 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
861
862 return result ;
863
864} // end of SumOverInterPlasmon
865
866///////////////////////////////////////////////////////////////////////////////
867//
868// Integration of PAI cross-section for the case of
869// passing across border between intervals
870
871G4double G4PAIySection::SumOverBorder( G4int i ,
872 G4double en0 )
873{
874 G4double x0,x1,y0,yy1,a,c,d,e0,result ;
875
876 e0 = en0 ;
877 x0 = fSplineEnergy[i] ;
878 x1 = fSplineEnergy[i+1] ;
879 y0 = fDifPAIySection[i] ;
880 yy1 = fDifPAIySection[i+1] ;
881
882 c = x1/x0;
883 d = e0/x0;
884 a = log10(yy1/y0)/log10(x1/x0) ;
885
886 G4double b = 0.0;
887 if(a < 20.) b = y0/pow(x0,a) ;
888
889 a += 1 ;
890 if(a == 0)
891 {
892 result = b*log(x0/e0) ;
893 }
894 else
895 {
896 result = y0*(x0 - e0*pow(d,a-1))/a ;
897 }
898 a++ ;
899 if(a == 0)
900 {
901 fIntegralPAIySection[0] += b*log(x0/e0) ;
902 }
903 else
904 {
905 fIntegralPAIySection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
906 }
907 x0 = fSplineEnergy[i - 1] ;
908 x1 = fSplineEnergy[i - 2] ;
909 y0 = fDifPAIySection[i - 1] ;
910 yy1 = fDifPAIySection[i - 2] ;
911
912 c = x1/x0;
913 d = e0/x0;
914 a = log10(yy1/y0)/log10(x1/x0) ;
915 // b0 = log10(y0) - a*log10(x0) ;
916 b = y0/pow(x0,a) ;
917 a += 1 ;
918 if(a == 0)
919 {
920 result += b*log(e0/x0) ;
921 }
922 else
923 {
924 result += y0*(e0*pow(d,a-1) - x0)/a ;
925 }
926 a++ ;
927 if(a == 0)
928 {
929 fIntegralPAIySection[0] += b*log(e0/x0) ;
930 }
931 else
932 {
933 fIntegralPAIySection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
934 }
935 return result ;
936
937}
938
939///////////////////////////////////////////////////////////////////////
940
941G4double G4PAIySection::SumOverBorderdEdx( G4int i ,
942 G4double en0 )
943{
944 G4double x0,x1,y0,yy1,a,c,d,e0,result ;
945
946 e0 = en0 ;
947 x0 = fSplineEnergy[i] ;
948 x1 = fSplineEnergy[i+1] ;
949 y0 = fDifPAIySection[i] ;
950 yy1 = fDifPAIySection[i+1] ;
951
952 c = x1/x0;
953 d = e0/x0;
954 a = log10(yy1/y0)/log10(x1/x0) ;
955
956 G4double b = 0.0;
957 if(a < 20.) b = y0/pow(x0,a) ;
958
959 a += 2 ;
960 if(a == 0)
961 {
962 result = b*log(x0/e0) ;
963 }
964 else
965 {
966 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
967 }
968 x0 = fSplineEnergy[i - 1] ;
969 x1 = fSplineEnergy[i - 2] ;
970 y0 = fDifPAIySection[i - 1] ;
971 yy1 = fDifPAIySection[i - 2] ;
972
973 c = x1/x0;
974 d = e0/x0;
975 a = log10(yy1/y0)/log10(x1/x0) ;
976
977 if(a < 20.) b = y0/pow(x0,a) ;
978
979 a += 2 ;
980 if(a == 0)
981 {
982 result += b*log(e0/x0) ;
983 }
984 else
985 {
986 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
987 }
988 return result ;
989
990}
991
992///////////////////////////////////////////////////////////////////////////////
993//
994// Integration of Cerenkov cross-section for the case of
995// passing across border between intervals
996
997G4double G4PAIySection::SumOverBordCerenkov( G4int i ,
998 G4double en0 )
999{
1000 G4double x0,x1,y0,yy1,a,e0,c,d,result ;
1001
1002 e0 = en0 ;
1003 x0 = fSplineEnergy[i] ;
1004 x1 = fSplineEnergy[i+1] ;
1005 y0 = fdNdxCerenkov[i] ;
1006 yy1 = fdNdxCerenkov[i+1] ;
1007
1008 // G4cout<<G4endl;
1009 //G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1010 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1011 c = x1/x0 ;
1012 d = e0/x0 ;
1013 a = log10(yy1/y0)/log10(c) ;
1014
1015 G4double b = 0.0;
1016 if(a < 20.) b = y0/pow(x0,a) ;
1017
1018 a += 1.0 ;
1019 if( a == 0 ) result = b*log(x0/e0) ;
1020 else result = y0*(x0 - e0*pow(d,a-1))/a ;
1021 a += 1.0 ;
1022
1023 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0) ;
1024 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1025
1026 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1027
1028 x0 = fSplineEnergy[i - 1] ;
1029 x1 = fSplineEnergy[i - 2] ;
1030 y0 = fdNdxCerenkov[i - 1] ;
1031 yy1 = fdNdxCerenkov[i - 2] ;
1032
1033 //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1034 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1035
1036 c = x1/x0 ;
1037 d = e0/x0 ;
1038 a = log10(yy1/y0)/log10(x1/x0) ;
1039
1040 // G4cout << "a= " << a << G4endl;
1041 if(a < 20.) b = y0/pow(x0,a) ;
1042
1043 if(a > 20.0) b = 0.0;
1044 else b = y0/pow(x0,a); // pow(10.,b0) ;
1045
1046 //G4cout << "b= " << b << G4endl;
1047
1048 a += 1.0 ;
1049 if( a == 0 ) result += b*log(e0/x0) ;
1050 else result += y0*(e0*pow(d,a-1) - x0 )/a ;
1051 a += 1.0 ;
1052 //G4cout << "result= " << result << G4endl;
1053
1054 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0) ;
1055 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1056
1057 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1058
1059 return result ;
1060
1061}
1062
1063///////////////////////////////////////////////////////////////////////////////
1064//
1065// Integration of Plasmon cross-section for the case of
1066// passing across border between intervals
1067
1068G4double G4PAIySection::SumOverBordPlasmon( G4int i ,
1069 G4double en0 )
1070{
1071 G4double x0,x1,y0,yy1,a,c,d,e0,result ;
1072
1073 e0 = en0 ;
1074 x0 = fSplineEnergy[i] ;
1075 x1 = fSplineEnergy[i+1] ;
1076 y0 = fdNdxPlasmon[i] ;
1077 yy1 = fdNdxPlasmon[i+1] ;
1078
1079 c = x1/x0 ;
1080 d = e0/x0 ;
1081 a = log10(yy1/y0)/log10(c) ;
1082
1083 G4double b = 0.0;
1084 if(a < 20.) b = y0/pow(x0,a) ;
1085
1086 a += 1.0 ;
1087 if( a == 0 ) result = b*log(x0/e0) ;
1088 else result = y0*(x0 - e0*pow(d,a-1))/a ;
1089 a += 1.0 ;
1090
1091 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0) ;
1092 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1093
1094 x0 = fSplineEnergy[i - 1] ;
1095 x1 = fSplineEnergy[i - 2] ;
1096 y0 = fdNdxPlasmon[i - 1] ;
1097 yy1 = fdNdxPlasmon[i - 2] ;
1098
1099 c = x1/x0 ;
1100 d = e0/x0 ;
1101 a = log10(yy1/y0)/log10(c) ;
1102
1103 if(a < 20.) b = y0/pow(x0,a) ;
1104
1105 a += 1.0 ;
1106 if( a == 0 ) result += b*log(e0/x0) ;
1107 else result += y0*(e0*pow(d,a-1) - x0)/a ;
1108 a += 1.0 ;
1109
1110 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0) ;
1111 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1112
1113 return result ;
1114
1115}
1116
1117/////////////////////////////////////////////////////////////////////////
1118//
1119//
1120
1121G4double G4PAIySection::GetStepEnergyLoss( G4double step )
1122{
1123 G4int iTransfer ;
1124 G4long numOfCollisions ;
1125 G4double loss = 0.0 ;
1126 G4double meanNumber, position ;
1127
1128 // G4cout<<" G4PAIySection::GetStepEnergyLoss "<<G4endl ;
1129
1130
1131
1132 meanNumber = fIntegralPAIySection[1]*step ;
1133 numOfCollisions = G4Poisson(meanNumber) ;
1134
1135 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1136
1137 while(numOfCollisions)
1138 {
1139 position = fIntegralPAIySection[1]*G4UniformRand() ;
1140
1141 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1142 {
1143 if( position >= fIntegralPAIySection[iTransfer] ) break ;
1144 }
1145 loss += fSplineEnergy[iTransfer] ;
1146 numOfCollisions-- ;
1147 }
1148 // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl ;
1149
1150 return loss ;
1151}
1152
1153/////////////////////////////////////////////////////////////////////////
1154//
1155//
1156
1157G4double G4PAIySection::GetStepCerenkovLoss( G4double step )
1158{
1159 G4int iTransfer ;
1160 G4long numOfCollisions ;
1161 G4double loss = 0.0 ;
1162 G4double meanNumber, position ;
1163
1164 // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl ;
1165
1166
1167
1168 meanNumber = fIntegralCerenkov[1]*step ;
1169 numOfCollisions = G4Poisson(meanNumber) ;
1170
1171 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1172
1173 while(numOfCollisions)
1174 {
1175 position = fIntegralCerenkov[1]*G4UniformRand() ;
1176
1177 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1178 {
1179 if( position >= fIntegralCerenkov[iTransfer] ) break ;
1180 }
1181 loss += fSplineEnergy[iTransfer] ;
1182 numOfCollisions-- ;
1183 }
1184 // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl ;
1185
1186 return loss ;
1187}
1188
1189/////////////////////////////////////////////////////////////////////////
1190//
1191//
1192
1193G4double G4PAIySection::GetStepPlasmonLoss( G4double step )
1194{
1195 G4int iTransfer ;
1196 G4long numOfCollisions ;
1197 G4double loss = 0.0 ;
1198 G4double meanNumber, position ;
1199
1200 // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl ;
1201
1202
1203
1204 meanNumber = fIntegralPlasmon[1]*step ;
1205 numOfCollisions = G4Poisson(meanNumber) ;
1206
1207 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1208
1209 while(numOfCollisions)
1210 {
1211 position = fIntegralPlasmon[1]*G4UniformRand() ;
1212
1213 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1214 {
1215 if( position >= fIntegralPlasmon[iTransfer] ) break ;
1216 }
1217 loss += fSplineEnergy[iTransfer] ;
1218 numOfCollisions-- ;
1219 }
1220 // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl ;
1221
1222 return loss ;
1223}
1224
1225
1226
1227/////////////////////////////////////////////////////////////////////////////
1228//
1229// Init array of Lorentz factors
1230//
1231
1232G4int G4PAIySection::fNumberOfGammas = 111 ;
1233
1234const G4double G4PAIySection::fLorentzFactor[112] = // fNumberOfGammas+1
1235{
12360.0,
12371.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
12381.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
12391.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
12401.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
12412.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
12423.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
12435.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
12448.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
12451.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
12462.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
12475.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
12481.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
12491.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
12503.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
12516.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
12521.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
12532.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
12544.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
12558.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
12561.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
12573.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
12585.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
12591.065799e+05
1260} ;
1261
1262///////////////////////////////////////////////////////////////////////
1263//
1264// The number of gamma for creation of spline (near ion-min , G ~ 4 )
1265//
1266
1267const
1268G4int G4PAIySection::fRefGammaNumber = 29 ;
1269
1270
1271//
1272// end of G4PAIySection implementation file
1273//
1274////////////////////////////////////////////////////////////////////////////
1275
Note: See TracBrowser for help on using the repository browser.