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