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