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