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

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

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