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

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

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