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

Last change on this file since 1340 was 1340, checked in by garnier, 14 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.