Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4PAIySection.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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