ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIySection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PAIySection.cc
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 
71 //
72 // Constructor
73 //
74 
76 {
77  fSandia = 0;
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;
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 
106 //
107 //
108 
110 {
111  return fLorentzFactor[j];
112 }
113 
115 //
116 // Constructor with beta*gamma square value called from G4PAIModel
117 
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 
241 //
242 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
243 //
244 
246 {
247  G4int i, numberOfElements = 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 
274 //
275 // General control function for class G4PAIySection
276 //
277 
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 
325 //
326 // Shifting from borders to intervals Creation of first energy points
327 //
328 
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 =
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 
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 
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 
485 //
486 // Integration over electrons that could be considered
487 // quasi-free at energy transfer of interest
488 
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 
507 //
508 // Imaginary part of dielectric constant
509 // (G4int k - interval number, G4double en1 - energy point)
510 
512  G4double energy1 )
513 {
514  G4double energy2,energy3,energy4,result;
515 
516  energy2 = energy1*energy1;
517  energy3 = energy2*energy1;
518  energy4 = energy3*energy1;
519 
520  result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
521  result *=hbarc/energy1;
522 
523  return result;
524 
525 } // end of ImPartDielectricConst
526 
527 
529 //
530 // Real part of dielectric constant minus unit: epsilon_1 - 1
531 // (G4double enb - energy point)
532 //
533 
535 {
536  G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
537  c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
538 
539  x0 = enb;
540  result = 0;
541 
542  for(G4int i=1;i<=fIntervalNumber-1;i++)
543  {
544  x1 = fEnergyInterval[i];
545  x2 = fEnergyInterval[i+1];
546  xx1 = x1 - x0;
547  xx2 = x2 - x0;
548  xx12 = xx2/xx1;
549 
550  if(xx12<0)
551  {
552  xx12 = -xx12;
553  }
554  xln1 = log(x2/x1);
555  xln2 = log(xx12);
556  xln3 = log((x2 + x0)/(x1 + x0));
557  x02 = x0*x0;
558  x03 = x02*x0;
559  x04 = x03*x0;
560  x05 = x04*x0;
561  G4double x12 = x1*x2;
562  c1 = (x2 - x1)/x12;
563  c2 = (x2 - x1)*(x2 +x1)/(x12*x12);
564  c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/(x12*x12*x12);
565 
566  result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
567  result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
568  result -= fA3[i]*c2/2/x02;
569  result -= fA4[i]*c3/3/x02;
570 
571  cof1 = fA1[i]/x02 + fA3[i]/x04;
572  cof2 = fA2[i]/x03 + fA4[i]/x05;
573 
574  result += 0.5*(cof1 +cof2)*xln2;
575  result += 0.5*(cof1 - cof2)*xln3;
576  }
577  result *= 2*hbarc/pi;
578 
579  return result;
580 
581 } // end of RePartDielectricConst
582 
584 //
585 // PAI differential cross-section in terms of
586 // simplified Allison's equation
587 //
588 
590  G4double betaGammaSq )
591 {
592  G4double beta, be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
593  //G4double beta, be4;
594  //G4double be4;
595  // G4double betaBohr2 = fine_structure_const*fine_structure_const;
596  // G4double betaBohr4 = betaBohr2*betaBohr2*4.0;
597  be2 = betaGammaSq/(1 + betaGammaSq);
598  //be4 = be2*be2;
599  beta = sqrt(be2);
600  cof = 1;
601  x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
602 
603  if( betaGammaSq < 0.01 ) x2 = log(be2);
604  else
605  {
606  x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
607  (1/betaGammaSq - fRePartDielectricConst[i]) +
608  fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
609  }
610  if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
611  {
612  x6=0;
613  }
614  else
615  {
616  x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
617  x5 = -1 - fRePartDielectricConst[i] +
618  be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
619  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
620 
621  x7 = atan2(fImPartDielectricConst[i],x3);
622  x6 = x5 * x7;
623  }
624  // if(fImPartDielectricConst[i] == 0) x6 = 0;
625 
626  x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
627  // if( x4 < 0.0 ) x4 = 0.0;
628  x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
629  fImPartDielectricConst[i]*fImPartDielectricConst[i];
630 
631  result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
632  result = std::max(result, 1.0e-8);
633  result *= fine_structure_const/be2/pi;
634  // low energy correction
635 
636  G4double lowCof = fLowEnergyCof; // 6.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)?
637 
638  result *= (1 - exp(-beta/betaBohr/lowCof));
639 
640  // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr));
641  // result *= (1-exp(-be2/betaBohr2));
642  // result *= (1-exp(-be4/betaBohr4));
643  // if(fDensity >= 0.1)
644  if(x8 > 0.)
645  {
646  result /= x8;
647  }
648  return result;
649 
650 } // end of DifPAIySection
651 
653 //
654 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
655 
657  G4double betaGammaSq )
658 {
659  G4double logarithm, x3, x5, argument, modul2, dNdxC;
660  G4double be2, be4;
661 
662  //G4double cof = 1.0;
663 
664  be2 = betaGammaSq/(1 + betaGammaSq);
665  be4 = be2*be2;
666 
667  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
668  else
669  {
670  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
671  (1/betaGammaSq - fRePartDielectricConst[i]) +
672  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
673  logarithm += log(1+1.0/betaGammaSq);
674  }
675 
676  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
677  {
678  argument = 0.0;
679  }
680  else
681  {
682  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
683  x5 = -1.0 - fRePartDielectricConst[i] +
684  be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
685  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
686  if( x3 == 0.0 ) argument = 0.5*pi;
687  else argument = atan2(fImPartDielectricConst[i],x3);
688  argument *= x5 ;
689  }
690  dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
691 
692  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
693 
694  dNdxC *= fine_structure_const/be2/pi;
695 
696  dNdxC *= (1-exp(-be4/betaBohr4));
697 
698  // if(fDensity >= 0.1)
699  // {
700  modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
701  fImPartDielectricConst[i]*fImPartDielectricConst[i];
702  if(modul2 > 0.)
703  {
704  dNdxC /= modul2;
705  }
706  return dNdxC;
707 
708 } // end of PAIdNdxCerenkov
709 
711 //
712 // Calculation od dN/dx of collisions with creation of longitudinal EM
713 // excitations (plasmons, delta-electrons)
714 
716  G4double betaGammaSq )
717 {
718  G4double cof, resonance, modul2, dNdxP;
719  G4double be2, be4;
720 
721  cof = 1;
722 
723  be2 = betaGammaSq/(1 + betaGammaSq);
724  be4 = be2*be2;
725 
726  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
727  resonance *= fImPartDielectricConst[i]/hbarc;
728 
729  dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
730 
731  dNdxP = std::max(dNdxP, 1.0e-8);
732 
733  dNdxP *= fine_structure_const/be2/pi;
734  dNdxP *= (1-exp(-be4/betaBohr4));
735 
736 // if( fDensity >= 0.1 )
737 // {
738  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
739  fImPartDielectricConst[i]*fImPartDielectricConst[i];
740  if(modul2 > 0.)
741  {
742  dNdxP /= modul2;
743  }
744  return dNdxP;
745 
746 } // end of PAIdNdxPlasmon
747 
749 //
750 // Calculation of the PAI integral cross-section
751 // fIntegralPAIySection[1] = specific primary ionisation, 1/cm
752 // and fIntegralPAIySection[0] = mean energy loss per cm in keV/cm
753 
755 {
756  fIntegralPAIySection[fSplineNumber] = 0;
757  fIntegralPAIdEdx[fSplineNumber] = 0;
758  fIntegralPAIySection[0] = 0;
759  G4int k = fIntervalNumber -1;
760 
761  for(G4int i = fSplineNumber-1; i >= 1; i--)
762  {
763  if(fSplineEnergy[i] >= fEnergyInterval[k])
764  {
765  fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + SumOverInterval(i);
766  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
767  }
768  else
769  {
770  fIntegralPAIySection[i] = fIntegralPAIySection[i+1] +
771  SumOverBorder(i+1,fEnergyInterval[k]);
772  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
773  SumOverBorderdEdx(i+1,fEnergyInterval[k]);
774  k--;
775  }
776  }
777 } // end of IntegralPAIySection
778 
780 //
781 // Calculation of the PAI Cerenkov integral cross-section
782 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
783 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
784 
786 {
787  G4int i, k;
788  fIntegralCerenkov[fSplineNumber] = 0;
789  fIntegralCerenkov[0] = 0;
790  k = fIntervalNumber -1;
791 
792  for( i = fSplineNumber-1; i >= 1; i-- )
793  {
794  if(fSplineEnergy[i] >= fEnergyInterval[k])
795  {
796  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
797  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
798  }
799  else
800  {
801  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
802  SumOverBordCerenkov(i+1,fEnergyInterval[k]);
803  k--;
804  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
805  }
806  }
807 
808 } // end of IntegralCerenkov
809 
811 //
812 // Calculation of the PAI Plasmon integral cross-section
813 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
814 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
815 
817 {
818  fIntegralPlasmon[fSplineNumber] = 0;
819  fIntegralPlasmon[0] = 0;
820  G4int k = fIntervalNumber -1;
821  for(G4int i=fSplineNumber-1;i>=1;i--)
822  {
823  if(fSplineEnergy[i] >= fEnergyInterval[k])
824  {
825  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
826  }
827  else
828  {
829  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
830  SumOverBordPlasmon(i+1,fEnergyInterval[k]);
831  k--;
832  }
833  }
834 
835 } // end of IntegralPlasmon
836 
838 //
839 // Calculation the PAI integral cross-section inside
840 // of interval of continuous values of photo-ionisation
841 // cross-section. Parameter 'i' is the number of interval.
842 
844 {
845  G4double x0,x1,y0,yy1,a,b,c,result;
846 
847  x0 = fSplineEnergy[i];
848  x1 = fSplineEnergy[i+1];
849 
850  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
851 
852  y0 = fDifPAIySection[i];
853  yy1 = fDifPAIySection[i+1];
854  //G4cout << "## x0= " << x0 << " x1= " << x1 << G4endl;
855  c = x1/x0;
856  //G4cout << "c= " << c << " y0= " << y0 << " yy1= " << yy1 << G4endl;
857  a = log10(yy1/y0)/log10(c);
858  //G4cout << "a= " << a << G4endl;
859  // b = log10(y0) - a*log10(x0);
860  b = y0/pow(x0,a);
861  a += 1;
862  if(a == 0)
863  {
864  result = b*log(x1/x0);
865  }
866  else
867  {
868  result = y0*(x1*pow(c,a-1) - x0)/a;
869  }
870  a++;
871  if(a == 0)
872  {
873  fIntegralPAIySection[0] += b*log(x1/x0);
874  }
875  else
876  {
877  fIntegralPAIySection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
878  }
879  return result;
880 
881 } // end of SumOverInterval
882 
884 
886 {
887  G4double x0,x1,y0,yy1,a,b,c,result;
888 
889  x0 = fSplineEnergy[i];
890  x1 = fSplineEnergy[i+1];
891 
892  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
893 
894  y0 = fDifPAIySection[i];
895  yy1 = fDifPAIySection[i+1];
896  c = x1/x0;
897  a = log10(yy1/y0)/log10(c);
898  // b = log10(y0) - a*log10(x0);
899  b = y0/pow(x0,a);
900  a += 2;
901  if(a == 0)
902  {
903  result = b*log(x1/x0);
904  }
905  else
906  {
907  result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
908  }
909  return result;
910 
911 } // end of SumOverInterval
912 
914 //
915 // Calculation the PAI Cerenkov integral cross-section inside
916 // of interval of continuous values of photo-ionisation Cerenkov
917 // cross-section. Parameter 'i' is the number of interval.
918 
920 {
921  G4double x0,x1,y0,yy1,a,c,result;
922 
923  x0 = fSplineEnergy[i];
924  x1 = fSplineEnergy[i+1];
925 
926  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
927 
928  y0 = fdNdxCerenkov[i];
929  yy1 = fdNdxCerenkov[i+1];
930  // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
931  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
932 
933  c = x1/x0;
934  a = log10(yy1/y0)/log10(c);
935  G4double b = 0.0;
936  if(a < 20.) b = y0/pow(x0,a);
937 
938  a += 1.0;
939  if(a == 0) result = b*log(c);
940  else result = y0*(x1*pow(c,a-1) - x0)/a;
941  a += 1.0;
942 
943  if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
944  else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
945  // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
946  return result;
947 
948 } // end of SumOverInterCerenkov
949 
951 //
952 // Calculation the PAI Plasmon integral cross-section inside
953 // of interval of continuous values of photo-ionisation Plasmon
954 // cross-section. Parameter 'i' is the number of interval.
955 
957 {
958  G4double x0,x1,y0,yy1,a,c,result;
959 
960  x0 = fSplineEnergy[i];
961  x1 = fSplineEnergy[i+1];
962 
963  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
964 
965  y0 = fdNdxPlasmon[i];
966  yy1 = fdNdxPlasmon[i+1];
967  c = x1/x0;
968  a = log10(yy1/y0)/log10(c);
969 
970  G4double b = 0.0;
971  if(a < 20.) b = y0/pow(x0,a);
972 
973  a += 1.0;
974  if(a == 0) result = b*log(x1/x0);
975  else result = y0*(x1*pow(c,a-1) - x0)/a;
976  a += 1.0;
977 
978  if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
979  else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
980 
981  return result;
982 
983 } // end of SumOverInterPlasmon
984 
986 //
987 // Integration of PAI cross-section for the case of
988 // passing across border between intervals
989 
991  G4double en0 )
992 {
993  G4double x0,x1,y0,yy1,a,/*c,*/d,e0,result;
994 
995  e0 = en0;
996  x0 = fSplineEnergy[i];
997  x1 = fSplineEnergy[i+1];
998  y0 = fDifPAIySection[i];
999  yy1 = fDifPAIySection[i+1];
1000 
1001  //c = x1/x0;
1002  d = e0/x0;
1003  a = log10(yy1/y0)/log10(x1/x0);
1004 
1005  G4double b = 0.0;
1006  if(a < 20.) b = y0/pow(x0,a);
1007 
1008  a += 1;
1009  if(a == 0)
1010  {
1011  result = b*log(x0/e0);
1012  }
1013  else
1014  {
1015  result = y0*(x0 - e0*pow(d,a-1))/a;
1016  }
1017  a++;
1018  if(a == 0)
1019  {
1020  fIntegralPAIySection[0] += b*log(x0/e0);
1021  }
1022  else
1023  {
1024  fIntegralPAIySection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1025  }
1026  x0 = fSplineEnergy[i - 1];
1027  x1 = fSplineEnergy[i - 2];
1028  y0 = fDifPAIySection[i - 1];
1029  yy1 = fDifPAIySection[i - 2];
1030 
1031  //c = x1/x0;
1032  d = e0/x0;
1033  a = log10(yy1/y0)/log10(x1/x0);
1034  // b0 = log10(y0) - a*log10(x0);
1035  b = y0/pow(x0,a);
1036  a += 1;
1037  if(a == 0)
1038  {
1039  result += b*log(e0/x0);
1040  }
1041  else
1042  {
1043  result += y0*(e0*pow(d,a-1) - x0)/a;
1044  }
1045  a++;
1046  if(a == 0)
1047  {
1048  fIntegralPAIySection[0] += b*log(e0/x0);
1049  }
1050  else
1051  {
1052  fIntegralPAIySection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1053  }
1054  return result;
1055 
1056 }
1057 
1059 
1061  G4double en0 )
1062 {
1063  G4double x0,x1,y0,yy1,a,/*c,*/d,e0,result;
1064 
1065  e0 = en0;
1066  x0 = fSplineEnergy[i];
1067  x1 = fSplineEnergy[i+1];
1068  y0 = fDifPAIySection[i];
1069  yy1 = fDifPAIySection[i+1];
1070 
1071  //c = x1/x0;
1072  d = e0/x0;
1073  a = log10(yy1/y0)/log10(x1/x0);
1074 
1075  G4double b = 0.0;
1076  if(a < 20.) b = y0/pow(x0,a);
1077 
1078  a += 2;
1079  if(a == 0)
1080  {
1081  result = b*log(x0/e0);
1082  }
1083  else
1084  {
1085  result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1086  }
1087  x0 = fSplineEnergy[i - 1];
1088  x1 = fSplineEnergy[i - 2];
1089  y0 = fDifPAIySection[i - 1];
1090  yy1 = fDifPAIySection[i - 2];
1091 
1092  //c = x1/x0;
1093  d = e0/x0;
1094  a = log10(yy1/y0)/log10(x1/x0);
1095 
1096  if(a < 20.) b = y0/pow(x0,a);
1097 
1098  a += 2;
1099  if(a == 0)
1100  {
1101  result += b*log(e0/x0);
1102  }
1103  else
1104  {
1105  result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1106  }
1107  return result;
1108 
1109 }
1110 
1112 //
1113 // Integration of Cerenkov cross-section for the case of
1114 // passing across border between intervals
1115 
1117  G4double en0 )
1118 {
1119  G4double x0,x1,y0,yy1,a,e0,c,d,result;
1120 
1121  e0 = en0;
1122  x0 = fSplineEnergy[i];
1123  x1 = fSplineEnergy[i+1];
1124  y0 = fdNdxCerenkov[i];
1125  yy1 = fdNdxCerenkov[i+1];
1126 
1127  // G4cout<<G4endl;
1128  //G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1129  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1130  c = x1/x0;
1131  d = e0/x0;
1132  a = log10(yy1/y0)/log10(c);
1133 
1134  G4double b = 0.0;
1135  if(a < 20.) b = y0/pow(x0,a);
1136 
1137  a += 1.0;
1138  if( a == 0 ) result = b*log(x0/e0);
1139  else result = y0*(x0 - e0*pow(d,a-1))/a;
1140  a += 1.0;
1141 
1142  if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
1143  else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1144 
1145  //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1146 
1147  x0 = fSplineEnergy[i - 1];
1148  x1 = fSplineEnergy[i - 2];
1149  y0 = fdNdxCerenkov[i - 1];
1150  yy1 = fdNdxCerenkov[i - 2];
1151 
1152  //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1153  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1154 
1155  c = x1/x0;
1156  d = e0/x0;
1157  a = log10(yy1/y0)/log10(x1/x0);
1158 
1159  // G4cout << "a= " << a << G4endl;
1160  if(a > 20.0) b = 0.0;
1161  else b = y0/pow(x0,a);
1162 
1163  //G4cout << "b= " << b << G4endl;
1164 
1165  a += 1.0;
1166  if( a == 0 ) result += b*log(e0/x0);
1167  else result += y0*(e0*pow(d,a-1) - x0 )/a;
1168  a += 1.0;
1169  //G4cout << "result= " << result << G4endl;
1170 
1171  if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
1172  else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1173 
1174  //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1175 
1176  return result;
1177 
1178 }
1179 
1181 //
1182 // Integration of Plasmon cross-section for the case of
1183 // passing across border between intervals
1184 
1186  G4double en0 )
1187 {
1188  G4double x0,x1,y0,yy1,a,c,d,e0,result;
1189 
1190  e0 = en0;
1191  x0 = fSplineEnergy[i];
1192  x1 = fSplineEnergy[i+1];
1193  y0 = fdNdxPlasmon[i];
1194  yy1 = fdNdxPlasmon[i+1];
1195 
1196  c = x1/x0;
1197  d = e0/x0;
1198  a = log10(yy1/y0)/log10(c);
1199 
1200  G4double b = 0.0;
1201  if(a < 20.) b = y0/pow(x0,a);
1202 
1203  a += 1.0;
1204  if( a == 0 ) result = b*log(x0/e0);
1205  else result = y0*(x0 - e0*pow(d,a-1))/a;
1206  a += 1.0;
1207 
1208  if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
1209  else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1210 
1211  x0 = fSplineEnergy[i - 1];
1212  x1 = fSplineEnergy[i - 2];
1213  y0 = fdNdxPlasmon[i - 1];
1214  yy1 = fdNdxPlasmon[i - 2];
1215 
1216  c = x1/x0;
1217  d = e0/x0;
1218  a = log10(yy1/y0)/log10(c);
1219 
1220  if(a < 20.) b = y0/pow(x0,a);
1221 
1222  a += 1.0;
1223  if( a == 0 ) result += b*log(e0/x0);
1224  else result += y0*(e0*pow(d,a-1) - x0)/a;
1225  a += 1.0;
1226 
1227  if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
1228  else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1229 
1230  return result;
1231 
1232 }
1233 
1235 //
1236 //
1237 
1239 {
1240  G4int iTransfer ;
1241  G4long numOfCollisions;
1242  G4double loss = 0.0;
1243  G4double meanNumber, position;
1244 
1245  // G4cout<<" G4PAIySection::GetStepEnergyLoss "<<G4endl;
1246 
1247 
1248 
1249  meanNumber = fIntegralPAIySection[1]*step;
1250  numOfCollisions = G4Poisson(meanNumber);
1251 
1252  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1253 
1254  while(numOfCollisions)
1255  {
1256  position = fIntegralPAIySection[1]*G4UniformRand();
1257 
1258  for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1259  {
1260  if( position >= fIntegralPAIySection[iTransfer] ) break;
1261  }
1262  loss += fSplineEnergy[iTransfer] ;
1263  numOfCollisions--;
1264  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1265  }
1266  // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
1267 
1268  return loss;
1269 }
1270 
1272 //
1273 //
1274 
1276 {
1277  G4int iTransfer ;
1278  G4long numOfCollisions;
1279  G4double loss = 0.0;
1280  G4double meanNumber, position;
1281 
1282  // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl;
1283 
1284 
1285 
1286  meanNumber = fIntegralCerenkov[1]*step;
1287  numOfCollisions = G4Poisson(meanNumber);
1288 
1289  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1290 
1291  while(numOfCollisions)
1292  {
1293  position = fIntegralCerenkov[1]*G4UniformRand();
1294 
1295  for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1296  {
1297  if( position >= fIntegralCerenkov[iTransfer] ) break;
1298  }
1299  loss += fSplineEnergy[iTransfer] ;
1300  numOfCollisions--;
1301  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1302  }
1303  // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
1304 
1305  return loss;
1306 }
1307 
1309 //
1310 //
1311 
1313 {
1314  G4int iTransfer ;
1315  G4long numOfCollisions;
1316  G4double loss = 0.0;
1317  G4double meanNumber, position;
1318 
1319  // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl;
1320 
1321 
1322 
1323  meanNumber = fIntegralPlasmon[1]*step;
1324  numOfCollisions = G4Poisson(meanNumber);
1325 
1326  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1327 
1328  while(numOfCollisions)
1329  {
1330  position = fIntegralPlasmon[1]*G4UniformRand();
1331 
1332  for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1333  {
1334  if( position >= fIntegralPlasmon[iTransfer] ) break;
1335  }
1336  loss += fSplineEnergy[iTransfer] ;
1337  numOfCollisions--;
1338  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1339  }
1340  // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
1341 
1342  return loss;
1343 }
1344 
1346 //
1347 
1348 void G4PAIySection::CallError(G4int i, const G4String& methodName) const
1349 {
1350  G4String head = "G4PAIySection::" + methodName + "()";
1352  ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
1353  G4Exception(head,"pai001",FatalException,ed);
1354 }
1355 
1357 //
1358 // Init array of Lorentz factors
1359 //
1360 
1362 
1363 const G4double G4PAIySection::fLorentzFactor[112] = // fNumberOfGammas+1
1364 {
1365 0.0,
1366 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
1367 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
1368 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
1369 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
1370 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
1371 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
1372 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
1373 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
1374 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
1375 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
1376 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
1377 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
1378 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
1379 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
1380 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
1381 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
1382 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
1383 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
1384 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
1385 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
1386 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
1387 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
1388 1.065799e+05
1389 };
1390 
1392 //
1393 // The number of gamma for creation of spline (near ion-min , G ~ 4 )
1394 //
1395 
1397 
1398 //
1399 // end of G4PAIySection implementation file
1400 //
1402