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