Geant4  10.00.p01
G4PAIxSection.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 //
27 // $Id: G4PAIxSection.cc 79188 2014-02-20 09:22:48Z gcosmo $
28 // GEANT4 tag $Name: geant4-09-03-ref-06 $
29 //
30 //
31 // G4PAIxSection.cc -- class implementation file
32 //
33 // GEANT 4 class implementation file
34 //
35 // For information related to this code, please, contact
36 // the Geant4 Collaboration.
37 //
38 // R&D: Vladimir.Grichine@cern.ch
39 //
40 // History:
41 //
42 // 13.05.03 V. Grichine, bug fixed for maxEnergyTransfer > max interval energy
43 // 28.05.01 V.Ivanchenko minor changes to provide ANSI -wall compilation
44 // 17.05.01 V. Grichine, low energy extension down to 10*keV of proton
45 // 20.11.98 adapted to a new Material/SandiaTable interface, mma
46 // 11.06.97 V. Grichine, 1st version
47 //
48 
49 
50 
51 #include "G4PAIxSection.hh"
52 
53 #include "globals.hh"
54 #include "G4PhysicalConstants.hh"
55 #include "G4SystemOfUnits.hh"
56 #include "G4ios.hh"
57 #include "G4Poisson.hh"
58 #include "G4Material.hh"
59 #include "G4MaterialCutsCouple.hh"
60 #include "G4SandiaTable.hh"
61 
62 using namespace std;
63 
64 /* ******************************************************************
65 
66 // Init array of Lorentz factors
67 
68 const G4double G4PAIxSection::fLorentzFactor[22] =
69 {
70  0.0 , 1.1 , 1.2 , 1.3 , 1.5 , 1.8 , 2.0 ,
71  2.5 , 3.0 , 4.0 , 7.0 , 10.0 , 20.0 , 40.0 ,
72  70.0 , 100.0 , 300.0 , 600.0 , 1000.0 , 3000.0 ,
73  10000.0 , 50000.0
74 };
75 
76 const G4int G4PAIxSection::
77 fRefGammaNumber = 29; // The number of gamma for creation of
78  // spline (9)
79 
80 ***************************************************************** */
81 
82 // Local class constants
83 
84 const G4double G4PAIxSection::fDelta = 0.005; // 0.005 energy shift from interval border
85 const G4double G4PAIxSection::fError = 0.005; // 0.005 error in lin-log approximation
86 
87 const G4int G4PAIxSection::fMaxSplineSize = 1000; // Max size of output spline
88  // arrays
90 //
91 // Constructor
92 //
93 
95 {
96  fSandia = 0;
97  fMatSandiaMatrix = 0;
98  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
99  fIntervalNumber = fSplineNumber = 0;
100  fVerbose = 0;
101 
102  fSplineEnergy = G4DataVector(fMaxSplineSize,0.0);
103  fRePartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
104  fImPartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
105  fIntegralTerm = G4DataVector(fMaxSplineSize,0.0);
106  fDifPAIxSection = G4DataVector(fMaxSplineSize,0.0);
107  fdNdxCerenkov = G4DataVector(fMaxSplineSize,0.0);
108  fdNdxPlasmon = G4DataVector(fMaxSplineSize,0.0);
109  fdNdxMM = G4DataVector(fMaxSplineSize,0.0);
110  fdNdxResonance = G4DataVector(fMaxSplineSize,0.0);
111  fIntegralPAIxSection = G4DataVector(fMaxSplineSize,0.0);
112  fIntegralPAIdEdx = G4DataVector(fMaxSplineSize,0.0);
113  fIntegralCerenkov = G4DataVector(fMaxSplineSize,0.0);
114  fIntegralPlasmon = G4DataVector(fMaxSplineSize,0.0);
115  fIntegralMM = G4DataVector(fMaxSplineSize,0.0);
116  fIntegralResonance = G4DataVector(fMaxSplineSize,0.0);
117 
118  fMaterialIndex = 0;
119 
120  for( G4int i = 0; i < 500; ++i )
121  {
122  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
123  }
124 }
125 
127 //
128 // Constructor
129 //
130 
132 {
133  fDensity = matCC->GetMaterial()->GetDensity();
134  G4int matIndex = matCC->GetMaterial()->GetIndex();
135  fMaterialIndex = matIndex;
136  fSandia = new G4SandiaTable(matIndex);
137  fVerbose = 0;
138 
139  G4int i, j;
140  fMatSandiaMatrix = new G4OrderedTable();
141 
142  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
143  {
144  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
145  }
146  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
147  {
148  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
149 
150  for(j = 1; j < 5; j++)
151  {
152  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
153  }
154  }
155  ComputeLowEnergyCof();
156  // fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0;
157 }
158 
160 
162  G4double maxEnergyTransfer)
163 {
164  fSandia = 0;
165  fMatSandiaMatrix = 0;
166  fVerbose = 0;
167  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
168  G4int i, j;
169 
170  fMaterialIndex = materialIndex;
171  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
172  fElectronDensity = (*theMaterialTable)[materialIndex]->
173  GetElectronDensity();
174  fIntervalNumber = (*theMaterialTable)[materialIndex]->
175  GetSandiaTable()->GetMatNbOfIntervals();
176  fIntervalNumber--;
177  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
178 
179  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
180  fA1 = G4DataVector(fIntervalNumber+2,0.0);
181  fA2 = G4DataVector(fIntervalNumber+2,0.0);
182  fA3 = G4DataVector(fIntervalNumber+2,0.0);
183  fA4 = G4DataVector(fIntervalNumber+2,0.0);
184 
185  for(i = 1; i <= fIntervalNumber; i++ )
186  {
187  if(((*theMaterialTable)[materialIndex]->
188  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
189  i > fIntervalNumber )
190  {
191  fEnergyInterval[i] = maxEnergyTransfer;
192  fIntervalNumber = i;
193  break;
194  }
195  fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
196  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
197  fA1[i] = (*theMaterialTable)[materialIndex]->
198  GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
199  fA2[i] = (*theMaterialTable)[materialIndex]->
200  GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
201  fA3[i] = (*theMaterialTable)[materialIndex]->
202  GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
203  fA4[i] = (*theMaterialTable)[materialIndex]->
204  GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
205  // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
206  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
207  }
208  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
209  {
210  fIntervalNumber++;
211  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
212  }
213 
214  // Now checking, if two borders are too close together
215 
216  for(i=1;i<fIntervalNumber;i++)
217  {
218  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
219  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
220  {
221  continue;
222  }
223  else
224  {
225  for(j=i;j<fIntervalNumber;j++)
226  {
227  fEnergyInterval[j] = fEnergyInterval[j+1];
228  fA1[j] = fA1[j+1];
229  fA2[j] = fA2[j+1];
230  fA3[j] = fA3[j+1];
231  fA4[j] = fA4[j+1];
232  }
233  fIntervalNumber--;
234  i--;
235  }
236  }
237 
238 
239  /* *********************************
240 
241  fSplineEnergy = new G4double[fMaxSplineSize];
242  fRePartDielectricConst = new G4double[fMaxSplineSize];
243  fImPartDielectricConst = new G4double[fMaxSplineSize];
244  fIntegralTerm = new G4double[fMaxSplineSize];
245  fDifPAIxSection = new G4double[fMaxSplineSize];
246  fIntegralPAIxSection = new G4double[fMaxSplineSize];
247 
248  for(i=0;i<fMaxSplineSize;i++)
249  {
250  fSplineEnergy[i] = 0.0;
251  fRePartDielectricConst[i] = 0.0;
252  fImPartDielectricConst[i] = 0.0;
253  fIntegralTerm[i] = 0.0;
254  fDifPAIxSection[i] = 0.0;
255  fIntegralPAIxSection[i] = 0.0;
256  }
257  ************************************************** */
258  ComputeLowEnergyCof();
259  InitPAI(); // create arrays allocated above
260  /*
261  delete[] fEnergyInterval;
262  delete[] fA1;
263  delete[] fA2;
264  delete[] fA3;
265  delete[] fA4;
266  */
267 }
268 
270 //
271 // Constructor called from G4PAIPhotonModel !!!
272 
274  G4double maxEnergyTransfer,
275  G4double betaGammaSq,
276  G4double** photoAbsCof,
277  G4int intNumber )
278 {
279  fSandia = 0;
280  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
281  fIntervalNumber = fSplineNumber = 0;
282  fVerbose = 0;
283 
284  fSplineEnergy = G4DataVector(500,0.0);
285  fRePartDielectricConst = G4DataVector(500,0.0);
286  fImPartDielectricConst = G4DataVector(500,0.0);
287  fIntegralTerm = G4DataVector(500,0.0);
288  fDifPAIxSection = G4DataVector(500,0.0);
289  fdNdxCerenkov = G4DataVector(500,0.0);
290  fdNdxPlasmon = G4DataVector(500,0.0);
291  fdNdxMM = G4DataVector(500,0.0);
292  fdNdxResonance = G4DataVector(500,0.0);
293  fIntegralPAIxSection = G4DataVector(500,0.0);
294  fIntegralPAIdEdx = G4DataVector(500,0.0);
295  fIntegralCerenkov = G4DataVector(500,0.0);
296  fIntegralPlasmon = G4DataVector(500,0.0);
297  fIntegralMM = G4DataVector(500,0.0);
298  fIntegralResonance = G4DataVector(500,0.0);
299 
300  for( G4int i = 0; i < 500; ++i )
301  {
302  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
303  }
304 
305  fSandia = 0;
306  fMatSandiaMatrix = 0;
307  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
308  G4int i, j;
309 
310  fMaterialIndex = materialIndex;
311  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
312  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
313 
314  fIntervalNumber = intNumber;
315  fIntervalNumber--;
316  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
317 
318  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
319  fA1 = G4DataVector(fIntervalNumber+2,0.0);
320  fA2 = G4DataVector(fIntervalNumber+2,0.0);
321  fA3 = G4DataVector(fIntervalNumber+2,0.0);
322  fA4 = G4DataVector(fIntervalNumber+2,0.0);
323 
324 
325  /*
326  fEnergyInterval = new G4double[fIntervalNumber+2];
327  fA1 = new G4double[fIntervalNumber+2];
328  fA2 = new G4double[fIntervalNumber+2];
329  fA3 = new G4double[fIntervalNumber+2];
330  fA4 = new G4double[fIntervalNumber+2];
331  */
332  for( i = 1; i <= fIntervalNumber; i++ )
333  {
334  if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
335  i > fIntervalNumber )
336  {
337  fEnergyInterval[i] = maxEnergyTransfer;
338  fIntervalNumber = i;
339  break;
340  }
341  fEnergyInterval[i] = photoAbsCof[i-1][0];
342  fA1[i] = photoAbsCof[i-1][1];
343  fA2[i] = photoAbsCof[i-1][2];
344  fA3[i] = photoAbsCof[i-1][3];
345  fA4[i] = photoAbsCof[i-1][4];
346  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
347  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
348  }
349  // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
350 
351  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
352  {
353  fIntervalNumber++;
354  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
355  }
356  // G4cout<<"after check of max energy transfer"<<G4endl;
357 
358  for( i = 1; i <= fIntervalNumber; i++ )
359  {
360  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
361  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
362  }
363  // Now checking, if two borders are too close together
364 
365  for( i = 1; i < fIntervalNumber; i++ )
366  {
367  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
368  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
369  {
370  continue;
371  }
372  else
373  {
374  for(j=i;j<fIntervalNumber;j++)
375  {
376  fEnergyInterval[j] = fEnergyInterval[j+1];
377  fA1[j] = fA1[j+1];
378  fA2[j] = fA2[j+1];
379  fA3[j] = fA3[j+1];
380  fA4[j] = fA4[j+1];
381  }
382  fIntervalNumber--;
383  i--;
384  }
385  }
386  // G4cout<<"after check of close borders"<<G4endl;
387 
388  for( i = 1; i <= fIntervalNumber; i++ )
389  {
390  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
391  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
392  }
393 
394  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
395 
396  ComputeLowEnergyCof();
397  G4double betaGammaSqRef =
398  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
399 
400  NormShift(betaGammaSqRef);
401  SplainPAI(betaGammaSqRef);
402 
403  // Preparation of integral PAI cross section for input betaGammaSq
404 
405  for(i = 1; i <= fSplineNumber; i++)
406  {
407  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
408  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
409  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
410  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
411  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
412 
413  // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
414  // <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
415  }
416  IntegralCerenkov();
417  IntegralMM();
418  IntegralPlasmon();
419  IntegralResonance();
420  IntegralPAIxSection();
421  /*
422  delete[] fEnergyInterval;
423  delete[] fA1;
424  delete[] fA2;
425  delete[] fA3;
426  delete[] fA4;
427  */
428 }
429 
431 //
432 // Test Constructor with beta*gamma square value
433 
435  G4double maxEnergyTransfer,
436  G4double betaGammaSq )
437 {
438  fSandia = 0;
439  fMatSandiaMatrix = 0;
440  fVerbose = 0;
441  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
442 
443  G4int i, j, numberOfElements;
444 
445  fMaterialIndex = materialIndex;
446  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
447  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
448  numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements();
449 
450  G4int* thisMaterialZ = new G4int[numberOfElements];
451 
452  for( i = 0; i < numberOfElements; i++ )
453  {
454  thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
455  GetElement(i)->GetZ();
456  }
457  // fSandia = new G4SandiaTable(materialIndex);
458  fSandia = (*theMaterialTable)[materialIndex]->GetSandiaTable();
459  G4SandiaTable thisMaterialSandiaTable(materialIndex);
460  fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals(thisMaterialZ,
461  numberOfElements);
462  fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
463  ( thisMaterialZ ,
464  (*theMaterialTable)[materialIndex]->GetFractionVector() ,
465  numberOfElements,fIntervalNumber);
466 
467  fIntervalNumber--;
468 
469  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
470  fA1 = G4DataVector(fIntervalNumber+2,0.0);
471  fA2 = G4DataVector(fIntervalNumber+2,0.0);
472  fA3 = G4DataVector(fIntervalNumber+2,0.0);
473  fA4 = G4DataVector(fIntervalNumber+2,0.0);
474 
475  /*
476  fEnergyInterval = new G4double[fIntervalNumber+2];
477  fA1 = new G4double[fIntervalNumber+2];
478  fA2 = new G4double[fIntervalNumber+2];
479  fA3 = new G4double[fIntervalNumber+2];
480  fA4 = new G4double[fIntervalNumber+2];
481  */
482  for( i = 1; i <= fIntervalNumber; i++ )
483  {
484  if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
485  i > fIntervalNumber)
486  {
487  fEnergyInterval[i] = maxEnergyTransfer;
488  fIntervalNumber = i;
489  break;
490  }
491  fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0);
492  fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity;
493  fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity;
494  fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity;
495  fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity;
496 
497  }
498  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
499  {
500  fIntervalNumber++;
501  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
502  fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
503  fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
504  fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
505  fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
506  }
507  for(i=1;i<=fIntervalNumber;i++)
508  {
509  // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
510  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
511  }
512  // Now checking, if two borders are too close together
513 
514  for( i = 1; i < fIntervalNumber; i++ )
515  {
516  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
517  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
518  {
519  continue;
520  }
521  else
522  {
523  for( j = i; j < fIntervalNumber; j++ )
524  {
525  fEnergyInterval[j] = fEnergyInterval[j+1];
526  fA1[j] = fA1[j+1];
527  fA2[j] = fA2[j+1];
528  fA3[j] = fA3[j+1];
529  fA4[j] = fA4[j+1];
530  }
531  fIntervalNumber--;
532  i--;
533  }
534  }
535 
536  /* *********************************
537  fSplineEnergy = new G4double[fMaxSplineSize];
538  fRePartDielectricConst = new G4double[fMaxSplineSize];
539  fImPartDielectricConst = new G4double[fMaxSplineSize];
540  fIntegralTerm = new G4double[fMaxSplineSize];
541  fDifPAIxSection = new G4double[fMaxSplineSize];
542  fIntegralPAIxSection = new G4double[fMaxSplineSize];
543 
544  for(i=0;i<fMaxSplineSize;i++)
545  {
546  fSplineEnergy[i] = 0.0;
547  fRePartDielectricConst[i] = 0.0;
548  fImPartDielectricConst[i] = 0.0;
549  fIntegralTerm[i] = 0.0;
550  fDifPAIxSection[i] = 0.0;
551  fIntegralPAIxSection[i] = 0.0;
552  }
553  */
554 
555  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
556 
557  ComputeLowEnergyCof();
558  G4double betaGammaSqRef =
559  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
560 
561  NormShift(betaGammaSqRef);
562  SplainPAI(betaGammaSqRef);
563 
564  // Preparation of integral PAI cross section for input betaGammaSq
565 
566  for(i = 1; i <= fSplineNumber; i++)
567  {
568  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
569  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
570  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
571  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
572  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
573  }
574  IntegralPAIxSection();
575  IntegralCerenkov();
576  IntegralMM();
577  IntegralPlasmon();
578  IntegralResonance();
579 }
580 
582 //
583 // Destructor
584 
586 {
587  /* ************************
588  delete[] fSplineEnergy ;
589  delete[] fRePartDielectricConst;
590  delete[] fImPartDielectricConst;
591  delete[] fIntegralTerm ;
592  delete[] fDifPAIxSection ;
593  delete[] fIntegralPAIxSection ;
594  */
595  delete fSandia;
596  delete fMatSandiaMatrix;
597 }
598 
599 
600 
601 
603 //
604 // Constructor with beta*gamma square value called from G4PAIPhotModel/Data
605 
606 void G4PAIxSection::Initialize( const G4Material* material,
607  G4double maxEnergyTransfer,
608  G4double betaGammaSq,
609  G4SandiaTable* sandia)
610 {
611  if(fVerbose > 0)
612  {
613  G4cout<<G4endl;
614  G4cout<<"G4PAIxSection::Initialize(...,G4SandiaTable* sandia)"<<G4endl;
615  G4cout<<G4endl;
616  }
617  G4int i, j;
618 
619  fSandia = sandia;
620  fIntervalNumber = sandia->GetMaxInterval();
621  fDensity = material->GetDensity();
622  fElectronDensity = material->GetElectronDensity();
623 
624  // fIntervalNumber--;
625 
626  if( fVerbose > 0 )
627  {
628  G4cout<<"fDensity = "<<fDensity<<"\t"<<fElectronDensity<<"\t fIntervalNumber = "<<fIntervalNumber<<G4endl;
629  }
630  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
631  fA1 = G4DataVector(fIntervalNumber+2,0.0);
632  fA2 = G4DataVector(fIntervalNumber+2,0.0);
633  fA3 = G4DataVector(fIntervalNumber+2,0.0);
634  fA4 = G4DataVector(fIntervalNumber+2,0.0);
635 
636  for( i = 1; i <= fIntervalNumber; i++ )
637  {
638  if ( sandia->GetSandiaMatTablePAI(i-1,0) < 1.*eV && sandia->GetLowerI1() == false)
639  {
640  fIntervalNumber--;
641  continue;
642  }
643  if( ( sandia->GetSandiaMatTablePAI(i-1,0) >= maxEnergyTransfer ) || i >= fIntervalNumber )
644  {
645  fEnergyInterval[i] = maxEnergyTransfer;
646  fIntervalNumber = i;
647  break;
648  }
649  fEnergyInterval[i] = sandia->GetSandiaMatTablePAI(i-1,0);
650  fA1[i] = sandia->GetSandiaMatTablePAI(i-1,1);
651  fA2[i] = sandia->GetSandiaMatTablePAI(i-1,2);
652  fA3[i] = sandia->GetSandiaMatTablePAI(i-1,3);
653  fA4[i] = sandia->GetSandiaMatTablePAI(i-1,4);
654 
655  if( fVerbose > 0 )
656  {
657  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
658  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
659  }
660  }
661  if( fVerbose > 0 ) G4cout<<"last i = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
662 
663  if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer )
664  {
665  fIntervalNumber++;
666  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
667  }
668  if( fVerbose > 0 )
669  {
670  for( i = 1; i <= fIntervalNumber; i++ )
671  {
672  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
673  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
674  }
675  }
676  if( fVerbose > 0 ) G4cout<<"Now checking, if two borders are too close together"<<G4endl;
677 
678  for( i = 1; i < fIntervalNumber; i++ )
679  {
680  if( fEnergyInterval[i+1]-fEnergyInterval[i] >
681  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]) && fEnergyInterval[i] > 0.) continue;
682  else
683  {
684  if( fVerbose > 0 ) G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fEnergyInterval[i+1]/keV;
685 
686  for( j = i; j < fIntervalNumber; j++ )
687  {
688  fEnergyInterval[j] = fEnergyInterval[j+1];
689  fA1[j] = fA1[j+1];
690  fA2[j] = fA2[j+1];
691  fA3[j] = fA3[j+1];
692  fA4[j] = fA4[j+1];
693  }
694  fIntervalNumber--;
695  i--;
696  }
697  }
698  if( fVerbose > 0 )
699  {
700  for( i = 1; i <= fIntervalNumber; i++ )
701  {
702  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
703  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
704  }
705  }
706  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
707 
708  ComputeLowEnergyCof(material);
709 
710  G4double betaGammaSqRef =
711  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
712 
713  NormShift(betaGammaSqRef);
714  SplainPAI(betaGammaSqRef);
715 
716  // Preparation of integral PAI cross section for input betaGammaSq
717 
718  for( i = 1; i <= fSplineNumber; i++ )
719  {
720  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
721 
722 
723  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
724  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
725  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
726  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
727  }
728  IntegralPAIxSection();
729  IntegralCerenkov();
730  IntegralMM();
731  IntegralPlasmon();
732  IntegralResonance();
733 
734  for( i = 1; i <= fSplineNumber; i++ )
735  {
736  if(fVerbose>0) G4cout<<i<<"; w = "<<fSplineEnergy[i]/keV<<" keV; dN/dx_>w = "<<fIntegralPAIxSection[i]<<" 1/mm"<<G4endl;
737  }
738 }
739 
740 
742 //
743 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
744 //
745 
747 {
748  G4int i, numberOfElements = material->GetNumberOfElements();
749  G4double sumZ = 0., sumCof = 0.;
750 
751  static const G4double p0 = 1.20923e+00;
752  static const G4double p1 = 3.53256e-01;
753  static const G4double p2 = -1.45052e-03;
754 
755  G4double* thisMaterialZ = new G4double[numberOfElements];
756  G4double* thisMaterialCof = new G4double[numberOfElements];
757 
758  for( i = 0; i < numberOfElements; i++ )
759  {
760  thisMaterialZ[i] = material->GetElement(i)->GetZ();
761  sumZ += thisMaterialZ[i];
762  thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
763  }
764  for( i = 0; i < numberOfElements; i++ )
765  {
766  sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
767  }
768  fLowEnergyCof = sumCof;
769  delete [] thisMaterialZ;
770  delete [] thisMaterialCof;
771  // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
772 }
773 
775 //
776 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
777 //
778 
780 {
781  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
782  G4int i, numberOfElements = (*theMaterialTable)[fMaterialIndex]->GetNumberOfElements();
783  G4double sumZ = 0., sumCof = 0.;
784 
785  const G4double p0 = 1.20923e+00;
786  const G4double p1 = 3.53256e-01;
787  const G4double p2 = -1.45052e-03;
788 
789  G4double* thisMaterialZ = new G4double[numberOfElements];
790  G4double* thisMaterialCof = new G4double[numberOfElements];
791 
792  for( i = 0; i < numberOfElements; i++ )
793  {
794  thisMaterialZ[i] = (*theMaterialTable)[fMaterialIndex]->GetElement(i)->GetZ();
795  sumZ += thisMaterialZ[i];
796  thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
797  }
798  for( i = 0; i < numberOfElements; i++ )
799  {
800  sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
801  }
802  fLowEnergyCof = sumCof;
803  // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
804  delete [] thisMaterialZ;
805  delete [] thisMaterialCof;
806 }
807 
809 //
810 // General control function for class G4PAIxSection
811 //
812 
814 {
815  G4int i;
816  G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
817  fLorentzFactor[fRefGammaNumber] - 1;
818 
819  // Preparation of integral PAI cross section for reference gamma
820 
821  NormShift(betaGammaSq);
822  SplainPAI(betaGammaSq);
823 
824  IntegralPAIxSection();
825  IntegralCerenkov();
826  IntegralMM();
827  IntegralPlasmon();
828  IntegralResonance();
829 
830  for(i = 0; i<= fSplineNumber; i++)
831  {
832  fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i];
833  if(i != 0)
834  {
835  fPAItable[i][0] = fSplineEnergy[i];
836  }
837  }
838  fPAItable[0][0] = fSplineNumber;
839 
840  for(G4int j = 1; j < 112; j++) // for other gammas
841  {
842  if( j == fRefGammaNumber ) continue;
843 
844  betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
845 
846  for(i = 1; i <= fSplineNumber; i++)
847  {
848  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
849  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
850  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
851  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
852  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
853  }
854  IntegralPAIxSection();
855  IntegralCerenkov();
856  IntegralMM();
857  IntegralPlasmon();
858  IntegralResonance();
859 
860  for(i = 0; i <= fSplineNumber; i++)
861  {
862  fPAItable[i][j] = fIntegralPAIxSection[i];
863  }
864  }
865 
866 }
867 
869 //
870 // Shifting from borders to intervals Creation of first energy points
871 //
872 
874 {
875  G4int i, j;
876 
877  if(fVerbose>0) G4cout<<" G4PAIxSection::NormShift call "<<G4endl;
878 
879 
880  for( i = 1; i <= fIntervalNumber-1; i++ )
881  {
882  for( j = 1; j <= 2; j++ )
883  {
884  fSplineNumber = (i-1)*2 + j;
885 
886  if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
887  else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
888  if(fVerbose>0) G4cout<<"cn = "<<fSplineNumber<<"; "<<"w = "<<fSplineEnergy[fSplineNumber]/keV<<" keV"<<G4endl;
889  }
890  }
891  fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
892 
893  j = 1;
894 
895  for( i = 2; i <= fSplineNumber; i++ )
896  {
897  if( fSplineEnergy[i]<fEnergyInterval[j+1] )
898  {
899  fIntegralTerm[i] = fIntegralTerm[i-1] +
900  RutherfordIntegral(j,fSplineEnergy[i-1],
901  fSplineEnergy[i] );
902  }
903  else
904  {
905  G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
906  fEnergyInterval[j+1] );
907  j++;
908  fIntegralTerm[i] = fIntegralTerm[i-1] + x +
909  RutherfordIntegral(j,fEnergyInterval[j],
910  fSplineEnergy[i] );
911  }
912  if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV \t"<<fIntegralTerm[i]<<"\n"<<G4endl;
913  }
914  fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
915  fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
916 
917  // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl;
918 
919  // Calculation of PAI differrential cross-section (1/(keV*cm))
920  // in the energy points near borders of energy intervals
921 
922  for(G4int k = 1; k <= fIntervalNumber-1; k++ )
923  {
924  for( j = 1; j <= 2; j++ )
925  {
926  i = (k-1)*2 + j;
927  fImPartDielectricConst[i] = fNormalizationCof*
928  ImPartDielectricConst(k,fSplineEnergy[i]);
929  fRePartDielectricConst[i] = fNormalizationCof*
930  RePartDielectricConst(fSplineEnergy[i]);
931  fIntegralTerm[i] *= fNormalizationCof;
932 
933  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
934  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
935  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
936  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
937  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
938  if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV, xsc = "<<fDifPAIxSection[i]<<"\n"<<G4endl;
939  }
940  }
941 
942 } // end of NormShift
943 
945 //
946 // Creation of new energy points as geometrical mean of existing
947 // one, calculation PAI_cs for them, while the error of logarithmic
948 // linear approximation would be smaller than 'fError'
949 
951 {
952  G4int j, k = 1, i = 1;
953 
954  if(fVerbose>0) G4cout<<" G4PAIxSection::SplainPAI call "<<G4endl;
955 
956  while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
957  {
958  // if( std::abs(fSplineEnergy[i+1]-fEnergyInterval[k+1]) > (fSplineEnergy[i+1]+fEnergyInterval[k+1])*5.e-7 )
959  if( fSplineEnergy[i+1] > fEnergyInterval[k+1] )
960  {
961  k++; // Here next energy point is in next energy interval
962  i++;
963  if(fVerbose>0) G4cout<<" in if: i = "<<i<<"; k = "<<k<<G4endl;
964  continue;
965  }
966  if(fVerbose>0) G4cout<<" out if: i = "<<i<<"; k = "<<k<<G4endl;
967 
968  // Shifting of arrayes for inserting the geometrical
969  // average of 'i' and 'i+1' energy points to 'i+1' place
970  fSplineNumber++;
971 
972  for( j = fSplineNumber; j >= i+2; j-- )
973  {
974  fSplineEnergy[j] = fSplineEnergy[j-1];
975  fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
976  fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
977  fIntegralTerm[j] = fIntegralTerm[j-1];
978 
979  fDifPAIxSection[j] = fDifPAIxSection[j-1];
980  fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
981  fdNdxMM[j] = fdNdxMM[j-1];
982  fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
983  fdNdxResonance[j] = fdNdxResonance[j-1];
984  }
985  G4double x1 = fSplineEnergy[i];
986  G4double x2 = fSplineEnergy[i+1];
987  G4double yy1 = fDifPAIxSection[i];
988  G4double y2 = fDifPAIxSection[i+1];
989 
990  if(fVerbose>0) G4cout<<"Spline: x1 = "<<x1<<"; x2 = "<<x2<<", yy1 = "<<yy1<<"; y2 = "<<y2<<G4endl;
991 
992 
993  G4double en1 = sqrt(x1*x2);
994  // G4double en1 = 0.5*(x1 + x2);
995 
996 
997  fSplineEnergy[i+1] = en1;
998 
999  // Calculation of logarithmic linear approximation
1000  // in this (enr) energy point, which number is 'i+1' now
1001 
1002  G4double a = log10(y2/yy1)/log10(x2/x1);
1003  G4double b = log10(yy1) - a*log10(x1);
1004  G4double y = a*log10(en1) + b;
1005 
1006  y = pow(10.,y);
1007 
1008  // Calculation of the PAI dif. cross-section at this point
1009 
1010  fImPartDielectricConst[i+1] = fNormalizationCof*
1011  ImPartDielectricConst(k,fSplineEnergy[i+1]);
1012  fRePartDielectricConst[i+1] = fNormalizationCof*
1013  RePartDielectricConst(fSplineEnergy[i+1]);
1014  fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
1015  RutherfordIntegral(k,fSplineEnergy[i],
1016  fSplineEnergy[i+1]);
1017 
1018  fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
1019  fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
1020  fdNdxMM[i+1] = PAIdNdxMM(i+1,betaGammaSq);
1021  fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
1022  fdNdxResonance[i+1] = PAIdNdxResonance(i+1,betaGammaSq);
1023 
1024  // Condition for next division of this segment or to pass
1025 
1026  if(fVerbose>0) G4cout<<"Spline, a = "<<a<<"; b = "<<b<<"; new xsc = "<<y<<"; compxsc = "<<fDifPAIxSection[i+1]<<G4endl;
1027 
1028  // to higher energies
1029 
1030  G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
1031 
1032  G4double delta = 2.*(fSplineEnergy[i+1]-fSplineEnergy[i])/(fSplineEnergy[i+1]+fSplineEnergy[i]);
1033 
1034  if( x < 0 )
1035  {
1036  x = -x;
1037  }
1038  if( x > fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*fDelta )
1039  {
1040  continue; // next division
1041  }
1042  i += 2; // pass to next segment
1043 
1044  } // close 'while'
1045 
1046 } // end of SplainPAI
1047 
1048 
1050 //
1051 // Integration over electrons that could be considered
1052 // quasi-free at energy transfer of interest
1053 
1055  G4double x1,
1056  G4double x2 )
1057 {
1058  G4double c1, c2, c3;
1059  // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
1060  c1 = (x2 - x1)/x1/x2;
1061  c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2;
1062  c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
1063  // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
1064 
1065  return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
1066 
1067 } // end of RutherfordIntegral
1068 
1069 
1071 //
1072 // Imaginary part of dielectric constant
1073 // (G4int k - interval number, G4double en1 - energy point)
1074 
1076  G4double energy1 )
1077 {
1078  G4double energy2,energy3,energy4,result;
1079 
1080  energy2 = energy1*energy1;
1081  energy3 = energy2*energy1;
1082  energy4 = energy3*energy1;
1083 
1084  result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
1085  result *=hbarc/energy1;
1086 
1087  return result;
1088 
1089 } // end of ImPartDielectricConst
1090 
1092 //
1093 // Returns lambda of photon with energy1 in current material
1094 
1096 {
1097  G4int i;
1098  G4double energy2, energy3, energy4, result, lambda;
1099 
1100  energy2 = energy1*energy1;
1101  energy3 = energy2*energy1;
1102  energy4 = energy3*energy1;
1103 
1104  // G4double* SandiaCof = fSandia->GetSandiaCofForMaterialPAI(energy1);
1105  // result = SandiaCof[0]/energy1+SandiaCof[1]/energy2+SandiaCof[2]/energy3+SandiaCof[3]/energy4;
1106  // result *= fDensity;
1107 
1108  for( i = 1; i <= fIntervalNumber; i++ )
1109  {
1110  if( energy1 < fEnergyInterval[i]) break;
1111  }
1112  i--;
1113  if(i == 0) i = 1;
1114 
1115  result = fA1[i]/energy1+fA2[i]/energy2+fA3[i]/energy3+fA4[i]/energy4;
1116 
1117  if( result > DBL_MIN ) lambda = 1./result;
1118  else lambda = DBL_MAX;
1119 
1120  return lambda;
1121 }
1122 
1124 //
1125 // Return lambda of electron with energy1 in current material
1126 // parametrisation from NIM A554(2005)474-493
1127 
1129 {
1130  G4double range;
1131  /*
1132  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1133 
1134  G4double Z = (*theMaterialTable)[fMaterialIndex]->GetIonisation()->GetZeffective();
1135  G4double A = (*theMaterialTable)[fMaterialIndex]->GetA();
1136 
1137  energy /= keV; // energy in keV in parametrised formula
1138 
1139  if (energy < 10.)
1140  {
1141  range = 3.872e-3*A/Z;
1142  range *= pow(energy,1.492);
1143  }
1144  else
1145  {
1146  range = 6.97e-3*pow(energy,1.6);
1147  }
1148  */
1149  // Blum&Rolandi Particle Detection with Drift Chambers, p. 7
1150 
1151  G4double cofA = 5.37e-4*g/cm2/keV;
1152  G4double cofB = 0.9815;
1153  G4double cofC = 3.123e-3/keV;
1154  // energy /= keV;
1155 
1156  range = cofA*energy*( 1 - cofB/(1 + cofC*energy) );
1157 
1158  // range *= g/cm2;
1159  range /= fDensity;
1160 
1161  return range;
1162 }
1163 
1165 //
1166 // Real part of dielectric constant minus unit: epsilon_1 - 1
1167 // (G4double enb - energy point)
1168 //
1169 
1171 {
1172  G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
1173  c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
1174 
1175  x0 = enb;
1176  result = 0;
1177 
1178  for(G4int i=1;i<=fIntervalNumber-1;i++)
1179  {
1180  x1 = fEnergyInterval[i];
1181  x2 = fEnergyInterval[i+1];
1182  xx1 = x1 - x0;
1183  xx2 = x2 - x0;
1184  xx12 = xx2/xx1;
1185 
1186  if(xx12<0)
1187  {
1188  xx12 = -xx12;
1189  }
1190  xln1 = log(x2/x1);
1191  xln2 = log(xx12);
1192  xln3 = log((x2 + x0)/(x1 + x0));
1193  x02 = x0*x0;
1194  x03 = x02*x0;
1195  x04 = x03*x0;
1196  x05 = x04*x0;
1197  c1 = (x2 - x1)/x1/x2;
1198  c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2;
1199  c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
1200 
1201  result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
1202  result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
1203  result -= fA3[i]*c2/2/x02;
1204  result -= fA4[i]*c3/3/x02;
1205 
1206  cof1 = fA1[i]/x02 + fA3[i]/x04;
1207  cof2 = fA2[i]/x03 + fA4[i]/x05;
1208 
1209  result += 0.5*(cof1 +cof2)*xln2;
1210  result += 0.5*(cof1 - cof2)*xln3;
1211  }
1212  result *= 2*hbarc/pi;
1213 
1214  return result;
1215 
1216 } // end of RePartDielectricConst
1217 
1219 //
1220 // PAI differential cross-section in terms of
1221 // simplified Allison's equation
1222 //
1223 
1225  G4double betaGammaSq )
1226 {
1227  G4double cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
1228 
1229  G4double betaBohr = fine_structure_const;
1230  // G4double betaBohr2 = fine_structure_const*fine_structure_const;
1231  // G4double betaBohr3 = betaBohr*betaBohr2; // *4.0;
1232 
1233  G4double be2 = betaGammaSq/(1 + betaGammaSq);
1234  G4double beta = sqrt(be2);
1235  // G4double be3 = beta*be2;
1236 
1237  cof = 1.;
1238  x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
1239 
1240  if( betaGammaSq < 0.01 ) x2 = log(be2);
1241  else
1242  {
1243  x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1244  (1/betaGammaSq - fRePartDielectricConst[i]) +
1245  fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
1246  }
1247  if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
1248  {
1249  x6 = 0.;
1250  }
1251  else
1252  {
1253  x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
1254  x5 = -1 - fRePartDielectricConst[i] +
1255  be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1256  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1257 
1258  x7 = atan2(fImPartDielectricConst[i],x3);
1259  x6 = x5 * x7;
1260  }
1261  // if(fImPartDielectricConst[i] == 0) x6 = 0.;
1262 
1263  x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
1264 
1265  // if( x4 < 0.0 ) x4 = 0.0;
1266 
1267  x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1268  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1269 
1270  result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
1271 
1272  if( result < 1.0e-8 ) result = 1.0e-8;
1273 
1274  result *= fine_structure_const/be2/pi;
1275 
1276  // low energy correction
1277 
1278  G4double lowCof = fLowEnergyCof; // 6.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)?
1279 
1280  result *= (1 - exp(-beta/betaBohr/lowCof));
1281 
1282 
1283  // result *= (1 - exp(-be2/betaBohr2/lowCof));
1284 
1285  // result *= (1 - exp(-be3/betaBohr3/lowCof)); // ~ be for be<<betaBohr
1286 
1287  // result *= (1 - exp(-be4/betaBohr4/lowCof));
1288 
1289  if(fDensity >= 0.1)
1290  {
1291  result /= x8;
1292  }
1293  return result;
1294 
1295 } // end of DifPAIxSection
1296 
1298 //
1299 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
1300 
1302  G4double betaGammaSq )
1303 {
1304  G4double logarithm, x3, x5, argument, modul2, dNdxC;
1306 
1307  cofBetaBohr = 4.0;
1308  betaBohr2 = fine_structure_const*fine_structure_const;
1309  G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1310 
1311  be2 = betaGammaSq/(1 + betaGammaSq);
1312  G4double be4 = be2*be2;
1313 
1314  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1315  else
1316  {
1317  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1318  (1/betaGammaSq - fRePartDielectricConst[i]) +
1319  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1320  logarithm += log(1+1.0/betaGammaSq);
1321  }
1322 
1323  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1324  {
1325  argument = 0.0;
1326  }
1327  else
1328  {
1329  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1330  x5 = -1.0 - fRePartDielectricConst[i] +
1331  be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1332  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1333  if( x3 == 0.0 ) argument = 0.5*pi;
1334  else argument = atan2(fImPartDielectricConst[i],x3);
1335  argument *= x5 ;
1336  }
1337  dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
1338 
1339  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1340 
1341  dNdxC *= fine_structure_const/be2/pi;
1342 
1343  dNdxC *= (1-exp(-be4/betaBohr4));
1344 
1345  if(fDensity >= 0.1)
1346  {
1347  modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1348  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1349  dNdxC /= modul2;
1350  }
1351  return dNdxC;
1352 
1353 } // end of PAIdNdxCerenkov
1354 
1356 //
1357 // Calculation od dN/dx of collisions of MM with creation of Cerenkov pseudo-photons
1358 
1360  G4double betaGammaSq )
1361 {
1362  G4double logarithm, x3, x5, argument, dNdxC;
1364 
1365  cofBetaBohr = 4.0;
1366  betaBohr2 = fine_structure_const*fine_structure_const;
1367  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1368 
1369  be2 = betaGammaSq/(1 + betaGammaSq);
1370  be4 = be2*be2;
1371 
1372  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1373  else
1374  {
1375  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1376  (1/betaGammaSq - fRePartDielectricConst[i]) +
1377  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1378  logarithm += log(1+1.0/betaGammaSq);
1379  }
1380 
1381  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1382  {
1383  argument = 0.0;
1384  }
1385  else
1386  {
1387  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1388  x5 = be2*( 1.0 + fRePartDielectricConst[i] ) - 1.0;
1389  if( x3 == 0.0 ) argument = 0.5*pi;
1390  else argument = atan2(fImPartDielectricConst[i],x3);
1391  argument *= x5 ;
1392  }
1393  dNdxC = ( logarithm*fImPartDielectricConst[i]*be2 + argument )/hbarc;
1394 
1395  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1396 
1397  dNdxC *= fine_structure_const/be2/pi;
1398 
1399  dNdxC *= (1-exp(-be4/betaBohr4));
1400  return dNdxC;
1401 
1402 } // end of PAIdNdxMM
1403 
1405 //
1406 // Calculation od dN/dx of collisions with creation of longitudinal EM
1407 // excitations (plasmons, delta-electrons)
1408 
1410  G4double betaGammaSq )
1411 {
1412  G4double resonance, modul2, dNdxP, cof = 1.;
1413  G4double be2, betaBohr;
1414 
1415  betaBohr = fine_structure_const;
1416  be2 = betaGammaSq/(1 + betaGammaSq);
1417 
1418  G4double beta = sqrt(be2);
1419 
1420  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1421  resonance *= fImPartDielectricConst[i]/hbarc;
1422 
1423 
1424  dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
1425 
1426  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1427 
1428  dNdxP *= fine_structure_const/be2/pi;
1429 
1430  dNdxP *= (1 - exp(-beta/betaBohr/fLowEnergyCof));
1431 
1432  // dNdxP *= (1-exp(-be4/betaBohr4));
1433 
1434  if( fDensity >= 0.1 )
1435  {
1436  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1437  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1438  dNdxP /= modul2;
1439  }
1440  return dNdxP;
1441 
1442 } // end of PAIdNdxPlasmon
1443 
1445 //
1446 // Calculation od dN/dx of collisions with creation of longitudinal EM
1447 // resonance excitations (plasmons, delta-electrons)
1448 
1450  G4double betaGammaSq )
1451 {
1452  G4double resonance, modul2, dNdxP;
1454 
1455  cofBetaBohr = 4.0;
1456  betaBohr2 = fine_structure_const*fine_structure_const;
1457  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1458 
1459  be2 = betaGammaSq/(1 + betaGammaSq);
1460  be4 = be2*be2;
1461 
1462  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1463  resonance *= fImPartDielectricConst[i]/hbarc;
1464 
1465 
1466  dNdxP = resonance;
1467 
1468  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1469 
1470  dNdxP *= fine_structure_const/be2/pi;
1471  dNdxP *= (1-exp(-be4/betaBohr4));
1472 
1473  if( fDensity >= 0.1 )
1474  {
1475  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1476  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1477  dNdxP /= modul2;
1478  }
1479  return dNdxP;
1480 
1481 } // end of PAIdNdxResonance
1482 
1484 //
1485 // Calculation of the PAI integral cross-section
1486 // fIntegralPAIxSection[1] = specific primary ionisation, 1/cm
1487 // and fIntegralPAIxSection[0] = mean energy loss per cm in keV/cm
1488 
1490 {
1491  fIntegralPAIxSection[fSplineNumber] = 0;
1492  fIntegralPAIdEdx[fSplineNumber] = 0;
1493  fIntegralPAIxSection[0] = 0;
1494  G4int i, k = fIntervalNumber -1;
1495 
1496  for( i = fSplineNumber-1; i >= 1; i--)
1497  {
1498  if(fSplineEnergy[i] >= fEnergyInterval[k])
1499  {
1500  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i);
1501  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
1502  }
1503  else
1504  {
1505  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] +
1506  SumOverBorder(i+1,fEnergyInterval[k]);
1507  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
1508  SumOverBorderdEdx(i+1,fEnergyInterval[k]);
1509  k--;
1510  }
1511  if(fVerbose>0) G4cout<<"i = "<<i<<"; k = "<<k<<"; intPAIxsc[i] = "<<fIntegralPAIxSection[i]<<G4endl;
1512  }
1513 } // end of IntegralPAIxSection
1514 
1516 //
1517 // Calculation of the PAI Cerenkov integral cross-section
1518 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
1519 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
1520 
1522 {
1523  G4int i, k;
1524  fIntegralCerenkov[fSplineNumber] = 0;
1525  fIntegralCerenkov[0] = 0;
1526  k = fIntervalNumber -1;
1527 
1528  for( i = fSplineNumber-1; i >= 1; i-- )
1529  {
1530  if(fSplineEnergy[i] >= fEnergyInterval[k])
1531  {
1532  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
1533  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1534  }
1535  else
1536  {
1537  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
1538  SumOverBordCerenkov(i+1,fEnergyInterval[k]);
1539  k--;
1540  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1541  }
1542  }
1543 
1544 } // end of IntegralCerenkov
1545 
1547 //
1548 // Calculation of the PAI MM-Cerenkov integral cross-section
1549 // fIntegralMM[1] = specific MM-Cerenkov ionisation, 1/cm
1550 // and fIntegralMM[0] = mean MM-Cerenkov loss per cm in keV/cm
1551 
1553 {
1554  G4int i, k;
1555  fIntegralMM[fSplineNumber] = 0;
1556  fIntegralMM[0] = 0;
1557  k = fIntervalNumber -1;
1558 
1559  for( i = fSplineNumber-1; i >= 1; i-- )
1560  {
1561  if(fSplineEnergy[i] >= fEnergyInterval[k])
1562  {
1563  fIntegralMM[i] = fIntegralMM[i+1] + SumOverInterMM(i);
1564  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1565  }
1566  else
1567  {
1568  fIntegralMM[i] = fIntegralMM[i+1] +
1569  SumOverBordMM(i+1,fEnergyInterval[k]);
1570  k--;
1571  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1572  }
1573  }
1574 
1575 } // end of IntegralMM
1576 
1578 //
1579 // Calculation of the PAI Plasmon integral cross-section
1580 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
1581 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
1582 
1584 {
1585  fIntegralPlasmon[fSplineNumber] = 0;
1586  fIntegralPlasmon[0] = 0;
1587  G4int k = fIntervalNumber -1;
1588  for(G4int i=fSplineNumber-1;i>=1;i--)
1589  {
1590  if(fSplineEnergy[i] >= fEnergyInterval[k])
1591  {
1592  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
1593  }
1594  else
1595  {
1596  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
1597  SumOverBordPlasmon(i+1,fEnergyInterval[k]);
1598  k--;
1599  }
1600  }
1601 
1602 } // end of IntegralPlasmon
1603 
1605 //
1606 // Calculation of the PAI resonance integral cross-section
1607 // fIntegralResonance[1] = resonance primary ionisation, 1/cm
1608 // and fIntegralResonance[0] = mean resonance loss per cm in keV/cm
1609 
1611 {
1612  fIntegralResonance[fSplineNumber] = 0;
1613  fIntegralResonance[0] = 0;
1614  G4int k = fIntervalNumber -1;
1615  for(G4int i=fSplineNumber-1;i>=1;i--)
1616  {
1617  if(fSplineEnergy[i] >= fEnergyInterval[k])
1618  {
1619  fIntegralResonance[i] = fIntegralResonance[i+1] + SumOverInterResonance(i);
1620  }
1621  else
1622  {
1623  fIntegralResonance[i] = fIntegralResonance[i+1] +
1624  SumOverBordResonance(i+1,fEnergyInterval[k]);
1625  k--;
1626  }
1627  }
1628 
1629 } // end of IntegralResonance
1630 
1632 //
1633 // Calculation the PAI integral cross-section inside
1634 // of interval of continuous values of photo-ionisation
1635 // cross-section. Parameter 'i' is the number of interval.
1636 
1638 {
1639  G4double x0,x1,y0,yy1,a,b,c,result;
1640 
1641  x0 = fSplineEnergy[i];
1642  x1 = fSplineEnergy[i+1];
1643 
1644  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1645 
1646  y0 = fDifPAIxSection[i];
1647  yy1 = fDifPAIxSection[i+1];
1648 
1649  if(fVerbose>0) G4cout<<"x0 = "<<x0<<"; x1 = "<<x1<<", y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1650 
1651  c = x1/x0;
1652  a = log10(yy1/y0)/log10(c);
1653 
1654  if(fVerbose>0) G4cout<<"SumOverInterval, a = "<<a<<"; c = "<<c<<G4endl;
1655 
1656  // b = log10(y0) - a*log10(x0);
1657  b = y0/pow(x0,a);
1658  a += 1.;
1659  if( std::fabs(a) < 1.e-6 )
1660  {
1661  result = b*log(x1/x0);
1662  }
1663  else
1664  {
1665  result = y0*(x1*pow(c,a-1) - x0)/a;
1666  }
1667  a += 1.;
1668  if( std::fabs(a) < 1.e-6 )
1669  {
1670  fIntegralPAIxSection[0] += b*log(x1/x0);
1671  }
1672  else
1673  {
1674  fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1675  }
1676  return result;
1677 
1678 } // end of SumOverInterval
1679 
1681 
1683 {
1684  G4double x0,x1,y0,yy1,a,b,c,result;
1685 
1686  x0 = fSplineEnergy[i];
1687  x1 = fSplineEnergy[i+1];
1688 
1689  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1690 
1691  y0 = fDifPAIxSection[i];
1692  yy1 = fDifPAIxSection[i+1];
1693  c = x1/x0;
1694  a = log10(yy1/y0)/log10(c);
1695  // b = log10(y0) - a*log10(x0);
1696  b = y0/pow(x0,a);
1697  a += 2;
1698  if(a == 0)
1699  {
1700  result = b*log(x1/x0);
1701  }
1702  else
1703  {
1704  result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1705  }
1706  return result;
1707 
1708 } // end of SumOverInterval
1709 
1711 //
1712 // Calculation the PAI Cerenkov integral cross-section inside
1713 // of interval of continuous values of photo-ionisation Cerenkov
1714 // cross-section. Parameter 'i' is the number of interval.
1715 
1717 {
1718  G4double x0,x1,y0,yy1,a,b,c,result;
1719 
1720  x0 = fSplineEnergy[i];
1721  x1 = fSplineEnergy[i+1];
1722 
1723  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1724 
1725  y0 = fdNdxCerenkov[i];
1726  yy1 = fdNdxCerenkov[i+1];
1727  // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1728  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1729 
1730  c = x1/x0;
1731  a = log10(yy1/y0)/log10(c);
1732  b = y0/pow(x0,a);
1733 
1734  a += 1.0;
1735  if(a == 0) result = b*log(c);
1736  else result = y0*(x1*pow(c,a-1) - x0)/a;
1737  a += 1.0;
1738 
1739  if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
1740  else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1741  // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1742  return result;
1743 
1744 } // end of SumOverInterCerenkov
1745 
1747 //
1748 // Calculation the PAI MM-Cerenkov integral cross-section inside
1749 // of interval of continuous values of photo-ionisation Cerenkov
1750 // cross-section. Parameter 'i' is the number of interval.
1751 
1753 {
1754  G4double x0,x1,y0,yy1,a,b,c,result;
1755 
1756  x0 = fSplineEnergy[i];
1757  x1 = fSplineEnergy[i+1];
1758 
1759  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1760 
1761  y0 = fdNdxMM[i];
1762  yy1 = fdNdxMM[i+1];
1763  // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1764  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1765 
1766  c = x1/x0;
1767  a = log10(yy1/y0)/log10(c);
1768  b = y0/pow(x0,a);
1769 
1770  a += 1.0;
1771  if(a == 0) result = b*log(c);
1772  else result = y0*(x1*pow(c,a-1) - x0)/a;
1773  a += 1.0;
1774 
1775  if( a == 0 ) fIntegralMM[0] += b*log(x1/x0);
1776  else fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1777  // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1778  return result;
1779 
1780 } // end of SumOverInterMM
1781 
1783 //
1784 // Calculation the PAI Plasmon integral cross-section inside
1785 // of interval of continuous values of photo-ionisation Plasmon
1786 // cross-section. Parameter 'i' is the number of interval.
1787 
1789 {
1790  G4double x0,x1,y0,yy1,a,b,c,result;
1791 
1792  x0 = fSplineEnergy[i];
1793  x1 = fSplineEnergy[i+1];
1794 
1795  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1796 
1797  y0 = fdNdxPlasmon[i];
1798  yy1 = fdNdxPlasmon[i+1];
1799  c =x1/x0;
1800  a = log10(yy1/y0)/log10(c);
1801  // b = log10(y0) - a*log10(x0);
1802  b = y0/pow(x0,a);
1803 
1804  a += 1.0;
1805  if(a == 0) result = b*log(x1/x0);
1806  else result = y0*(x1*pow(c,a-1) - x0)/a;
1807  a += 1.0;
1808 
1809  if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
1810  else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1811 
1812  return result;
1813 
1814 } // end of SumOverInterPlasmon
1815 
1817 //
1818 // Calculation the PAI resonance integral cross-section inside
1819 // of interval of continuous values of photo-ionisation resonance
1820 // cross-section. Parameter 'i' is the number of interval.
1821 
1823 {
1824  G4double x0,x1,y0,yy1,a,b,c,result;
1825 
1826  x0 = fSplineEnergy[i];
1827  x1 = fSplineEnergy[i+1];
1828 
1829  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1830 
1831  y0 = fdNdxResonance[i];
1832  yy1 = fdNdxResonance[i+1];
1833  c =x1/x0;
1834  a = log10(yy1/y0)/log10(c);
1835  // b = log10(y0) - a*log10(x0);
1836  b = y0/pow(x0,a);
1837 
1838  a += 1.0;
1839  if(a == 0) result = b*log(x1/x0);
1840  else result = y0*(x1*pow(c,a-1) - x0)/a;
1841  a += 1.0;
1842 
1843  if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0);
1844  else fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1845 
1846  return result;
1847 
1848 } // end of SumOverInterResonance
1849 
1851 //
1852 // Integration of PAI cross-section for the case of
1853 // passing across border between intervals
1854 
1856  G4double en0 )
1857 {
1858  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1859 
1860  e0 = en0;
1861  x0 = fSplineEnergy[i];
1862  x1 = fSplineEnergy[i+1];
1863  y0 = fDifPAIxSection[i];
1864  yy1 = fDifPAIxSection[i+1];
1865 
1866  //c = x1/x0;
1867  d = e0/x0;
1868  a = log10(yy1/y0)/log10(x1/x0);
1869 
1870  if(fVerbose>0) G4cout<<"SumOverBorder, a = "<<a<<G4endl;
1871 
1872  // b0 = log10(y0) - a*log10(x0);
1873  b = y0/pow(x0,a); // pow(10.,b);
1874 
1875  a += 1.;
1876  if( std::fabs(a) < 1.e-6 )
1877  {
1878  result = b*log(x0/e0);
1879  }
1880  else
1881  {
1882  result = y0*(x0 - e0*pow(d,a-1))/a;
1883  }
1884  a += 1.;
1885  if( std::fabs(a) < 1.e-6 )
1886  {
1887  fIntegralPAIxSection[0] += b*log(x0/e0);
1888  }
1889  else
1890  {
1891  fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1892  }
1893  x0 = fSplineEnergy[i - 1];
1894  x1 = fSplineEnergy[i - 2];
1895  y0 = fDifPAIxSection[i - 1];
1896  yy1 = fDifPAIxSection[i - 2];
1897 
1898  //c = x1/x0;
1899  d = e0/x0;
1900  a = log10(yy1/y0)/log10(x1/x0);
1901  // b0 = log10(y0) - a*log10(x0);
1902  b = y0/pow(x0,a);
1903  a += 1.;
1904  if( std::fabs(a) < 1.e-6 )
1905  {
1906  result += b*log(e0/x0);
1907  }
1908  else
1909  {
1910  result += y0*(e0*pow(d,a-1) - x0)/a;
1911  }
1912  a += 1.;
1913  if( std::fabs(a) < 1.e-6 )
1914  {
1915  fIntegralPAIxSection[0] += b*log(e0/x0);
1916  }
1917  else
1918  {
1919  fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1920  }
1921  return result;
1922 
1923 }
1924 
1926 
1928  G4double en0 )
1929 {
1930  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1931 
1932  e0 = en0;
1933  x0 = fSplineEnergy[i];
1934  x1 = fSplineEnergy[i+1];
1935  y0 = fDifPAIxSection[i];
1936  yy1 = fDifPAIxSection[i+1];
1937 
1938  //c = x1/x0;
1939  d = e0/x0;
1940  a = log10(yy1/y0)/log10(x1/x0);
1941  // b0 = log10(y0) - a*log10(x0);
1942  b = y0/pow(x0,a); // pow(10.,b);
1943 
1944  a += 2;
1945  if(a == 0)
1946  {
1947  result = b*log(x0/e0);
1948  }
1949  else
1950  {
1951  result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1952  }
1953  x0 = fSplineEnergy[i - 1];
1954  x1 = fSplineEnergy[i - 2];
1955  y0 = fDifPAIxSection[i - 1];
1956  yy1 = fDifPAIxSection[i - 2];
1957 
1958  // c = x1/x0;
1959  d = e0/x0;
1960  a = log10(yy1/y0)/log10(x1/x0);
1961  // b0 = log10(y0) - a*log10(x0);
1962  b = y0/pow(x0,a);
1963  a += 2;
1964  if(a == 0)
1965  {
1966  result += b*log(e0/x0);
1967  }
1968  else
1969  {
1970  result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1971  }
1972  return result;
1973 
1974 }
1975 
1977 //
1978 // Integration of Cerenkov cross-section for the case of
1979 // passing across border between intervals
1980 
1982  G4double en0 )
1983 {
1984  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1985 
1986  e0 = en0;
1987  x0 = fSplineEnergy[i];
1988  x1 = fSplineEnergy[i+1];
1989  y0 = fdNdxCerenkov[i];
1990  yy1 = fdNdxCerenkov[i+1];
1991 
1992  // G4cout<<G4endl;
1993  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1994  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1995  c = x1/x0;
1996  d = e0/x0;
1997  a = log10(yy1/y0)/log10(c);
1998  // b0 = log10(y0) - a*log10(x0);
1999  b = y0/pow(x0,a); // pow(10.,b0);
2000 
2001  a += 1.0;
2002  if( a == 0 ) result = b*log(x0/e0);
2003  else result = y0*(x0 - e0*pow(d,a-1))/a;
2004  a += 1.0;
2005 
2006  if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
2007  else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2008 
2009 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2010 
2011  x0 = fSplineEnergy[i - 1];
2012  x1 = fSplineEnergy[i - 2];
2013  y0 = fdNdxCerenkov[i - 1];
2014  yy1 = fdNdxCerenkov[i - 2];
2015 
2016  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2017  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2018 
2019  c = x1/x0;
2020  d = e0/x0;
2021  a = log10(yy1/y0)/log10(x1/x0);
2022  // b0 = log10(y0) - a*log10(x0);
2023  b = y0/pow(x0,a); // pow(10.,b0);
2024 
2025  a += 1.0;
2026  if( a == 0 ) result += b*log(e0/x0);
2027  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2028  a += 1.0;
2029 
2030  if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
2031  else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2032 
2033  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2034  // <<b<<"; result = "<<result<<G4endl;
2035 
2036  return result;
2037 
2038 }
2039 
2041 //
2042 // Integration of MM-Cerenkov cross-section for the case of
2043 // passing across border between intervals
2044 
2046  G4double en0 )
2047 {
2048  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
2049 
2050  e0 = en0;
2051  x0 = fSplineEnergy[i];
2052  x1 = fSplineEnergy[i+1];
2053  y0 = fdNdxMM[i];
2054  yy1 = fdNdxMM[i+1];
2055 
2056  // G4cout<<G4endl;
2057  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2058  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2059  c = x1/x0;
2060  d = e0/x0;
2061  a = log10(yy1/y0)/log10(c);
2062  // b0 = log10(y0) - a*log10(x0);
2063  b = y0/pow(x0,a); // pow(10.,b0);
2064 
2065  a += 1.0;
2066  if( a == 0 ) result = b*log(x0/e0);
2067  else result = y0*(x0 - e0*pow(d,a-1))/a;
2068  a += 1.0;
2069 
2070  if( a == 0 ) fIntegralMM[0] += b*log(x0/e0);
2071  else fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2072 
2073 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2074 
2075  x0 = fSplineEnergy[i - 1];
2076  x1 = fSplineEnergy[i - 2];
2077  y0 = fdNdxMM[i - 1];
2078  yy1 = fdNdxMM[i - 2];
2079 
2080  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2081  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2082 
2083  c = x1/x0;
2084  d = e0/x0;
2085  a = log10(yy1/y0)/log10(x1/x0);
2086  // b0 = log10(y0) - a*log10(x0);
2087  b = y0/pow(x0,a); // pow(10.,b0);
2088 
2089  a += 1.0;
2090  if( a == 0 ) result += b*log(e0/x0);
2091  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2092  a += 1.0;
2093 
2094  if( a == 0 ) fIntegralMM[0] += b*log(e0/x0);
2095  else fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2096 
2097  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2098  // <<b<<"; result = "<<result<<G4endl;
2099 
2100  return result;
2101 
2102 }
2103 
2105 //
2106 // Integration of Plasmon cross-section for the case of
2107 // passing across border between intervals
2108 
2110  G4double en0 )
2111 {
2112  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2113 
2114  e0 = en0;
2115  x0 = fSplineEnergy[i];
2116  x1 = fSplineEnergy[i+1];
2117  y0 = fdNdxPlasmon[i];
2118  yy1 = fdNdxPlasmon[i+1];
2119 
2120  c = x1/x0;
2121  d = e0/x0;
2122  a = log10(yy1/y0)/log10(c);
2123  // b0 = log10(y0) - a*log10(x0);
2124  b = y0/pow(x0,a); //pow(10.,b);
2125 
2126  a += 1.0;
2127  if( a == 0 ) result = b*log(x0/e0);
2128  else result = y0*(x0 - e0*pow(d,a-1))/a;
2129  a += 1.0;
2130 
2131  if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
2132  else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2133 
2134  x0 = fSplineEnergy[i - 1];
2135  x1 = fSplineEnergy[i - 2];
2136  y0 = fdNdxPlasmon[i - 1];
2137  yy1 = fdNdxPlasmon[i - 2];
2138 
2139  c = x1/x0;
2140  d = e0/x0;
2141  a = log10(yy1/y0)/log10(c);
2142  // b0 = log10(y0) - a*log10(x0);
2143  b = y0/pow(x0,a);// pow(10.,b0);
2144 
2145  a += 1.0;
2146  if( a == 0 ) result += b*log(e0/x0);
2147  else result += y0*(e0*pow(d,a-1) - x0)/a;
2148  a += 1.0;
2149 
2150  if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
2151  else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2152 
2153  return result;
2154 
2155 }
2156 
2158 //
2159 // Integration of resonance cross-section for the case of
2160 // passing across border between intervals
2161 
2163  G4double en0 )
2164 {
2165  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2166 
2167  e0 = en0;
2168  x0 = fSplineEnergy[i];
2169  x1 = fSplineEnergy[i+1];
2170  y0 = fdNdxResonance[i];
2171  yy1 = fdNdxResonance[i+1];
2172 
2173  c = x1/x0;
2174  d = e0/x0;
2175  a = log10(yy1/y0)/log10(c);
2176  // b0 = log10(y0) - a*log10(x0);
2177  b = y0/pow(x0,a); //pow(10.,b);
2178 
2179  a += 1.0;
2180  if( a == 0 ) result = b*log(x0/e0);
2181  else result = y0*(x0 - e0*pow(d,a-1))/a;
2182  a += 1.0;
2183 
2184  if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0);
2185  else fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2186 
2187  x0 = fSplineEnergy[i - 1];
2188  x1 = fSplineEnergy[i - 2];
2189  y0 = fdNdxResonance[i - 1];
2190  yy1 = fdNdxResonance[i - 2];
2191 
2192  c = x1/x0;
2193  d = e0/x0;
2194  a = log10(yy1/y0)/log10(c);
2195  // b0 = log10(y0) - a*log10(x0);
2196  b = y0/pow(x0,a);// pow(10.,b0);
2197 
2198  a += 1.0;
2199  if( a == 0 ) result += b*log(e0/x0);
2200  else result += y0*(e0*pow(d,a-1) - x0)/a;
2201  a += 1.0;
2202 
2203  if( a == 0 ) fIntegralResonance[0] += b*log(e0/x0);
2204  else fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2205 
2206  return result;
2207 
2208 }
2209 
2211 //
2212 // Returns random PAI-total energy loss over step
2213 
2215 {
2216  G4long numOfCollisions;
2217  G4double meanNumber, loss = 0.0;
2218 
2219  // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl;
2220 
2221  meanNumber = fIntegralPAIxSection[1]*step;
2222  numOfCollisions = G4Poisson(meanNumber);
2223 
2224  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2225 
2226  while(numOfCollisions)
2227  {
2228  loss += GetEnergyTransfer();
2229  numOfCollisions--;
2230  }
2231  // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
2232 
2233  return loss;
2234 }
2235 
2237 //
2238 // Returns random PAI-total energy transfer in one collision
2239 
2241 {
2242  G4int iTransfer ;
2243 
2244  G4double energyTransfer, position;
2245 
2246  position = fIntegralPAIxSection[1]*G4UniformRand();
2247 
2248  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2249  {
2250  if( position >= fIntegralPAIxSection[iTransfer] ) break;
2251  }
2252  if(iTransfer > fSplineNumber) iTransfer--;
2253 
2254  energyTransfer = fSplineEnergy[iTransfer];
2255 
2256  if(iTransfer > 1)
2257  {
2258  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2259  }
2260  return energyTransfer;
2261 }
2262 
2264 //
2265 // Returns random Cerenkov energy loss over step
2266 
2268 {
2269  G4long numOfCollisions;
2270  G4double meanNumber, loss = 0.0;
2271 
2272  // G4cout<<" G4PAIxSection::GetStepCerenkovLoss "<<G4endl;
2273 
2274  meanNumber = fIntegralCerenkov[1]*step;
2275  numOfCollisions = G4Poisson(meanNumber);
2276 
2277  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2278 
2279  while(numOfCollisions)
2280  {
2281  loss += GetCerenkovEnergyTransfer();
2282  numOfCollisions--;
2283  }
2284  // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2285 
2286  return loss;
2287 }
2288 
2290 //
2291 // Returns random MM-Cerenkov energy loss over step
2292 
2294 {
2295  G4long numOfCollisions;
2296  G4double meanNumber, loss = 0.0;
2297 
2298  // G4cout<<" G4PAIxSection::GetStepMMLoss "<<G4endl;
2299 
2300  meanNumber = fIntegralMM[1]*step;
2301  numOfCollisions = G4Poisson(meanNumber);
2302 
2303  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2304 
2305  while(numOfCollisions)
2306  {
2307  loss += GetMMEnergyTransfer();
2308  numOfCollisions--;
2309  }
2310  // G4cout<<"PAI MM-Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2311 
2312  return loss;
2313 }
2314 
2316 //
2317 // Returns Cerenkov energy transfer in one collision
2318 
2320 {
2321  G4int iTransfer ;
2322 
2323  G4double energyTransfer, position;
2324 
2325  position = fIntegralCerenkov[1]*G4UniformRand();
2326 
2327  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2328  {
2329  if( position >= fIntegralCerenkov[iTransfer] ) break;
2330  }
2331  if(iTransfer > fSplineNumber) iTransfer--;
2332 
2333  energyTransfer = fSplineEnergy[iTransfer];
2334 
2335  if(iTransfer > 1)
2336  {
2337  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2338  }
2339  return energyTransfer;
2340 }
2341 
2343 //
2344 // Returns MM-Cerenkov energy transfer in one collision
2345 
2347 {
2348  G4int iTransfer ;
2349 
2350  G4double energyTransfer, position;
2351 
2352  position = fIntegralMM[1]*G4UniformRand();
2353 
2354  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2355  {
2356  if( position >= fIntegralMM[iTransfer] ) break;
2357  }
2358  if(iTransfer > fSplineNumber) iTransfer--;
2359 
2360  energyTransfer = fSplineEnergy[iTransfer];
2361 
2362  if(iTransfer > 1)
2363  {
2364  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2365  }
2366  return energyTransfer;
2367 }
2368 
2370 //
2371 // Returns random plasmon energy loss over step
2372 
2374 {
2375  G4long numOfCollisions;
2376  G4double meanNumber, loss = 0.0;
2377 
2378  // G4cout<<" G4PAIxSection::GetStepPlasmonLoss "<<G4endl;
2379 
2380  meanNumber = fIntegralPlasmon[1]*step;
2381  numOfCollisions = G4Poisson(meanNumber);
2382 
2383  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2384 
2385  while(numOfCollisions)
2386  {
2387  loss += GetPlasmonEnergyTransfer();
2388  numOfCollisions--;
2389  }
2390  // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
2391 
2392  return loss;
2393 }
2394 
2396 //
2397 // Returns plasmon energy transfer in one collision
2398 
2400 {
2401  G4int iTransfer ;
2402 
2403  G4double energyTransfer, position;
2404 
2405  position = fIntegralPlasmon[1]*G4UniformRand();
2406 
2407  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2408  {
2409  if( position >= fIntegralPlasmon[iTransfer] ) break;
2410  }
2411  if(iTransfer > fSplineNumber) iTransfer--;
2412 
2413  energyTransfer = fSplineEnergy[iTransfer];
2414 
2415  if(iTransfer > 1)
2416  {
2417  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2418  }
2419  return energyTransfer;
2420 }
2421 
2423 //
2424 // Returns random resonance energy loss over step
2425 
2427 {
2428  G4long numOfCollisions;
2429  G4double meanNumber, loss = 0.0;
2430 
2431  // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl;
2432 
2433  meanNumber = fIntegralResonance[1]*step;
2434  numOfCollisions = G4Poisson(meanNumber);
2435 
2436  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2437 
2438  while(numOfCollisions)
2439  {
2440  loss += GetResonanceEnergyTransfer();
2441  numOfCollisions--;
2442  }
2443  // G4cout<<"PAI resonance loss = "<<loss/keV<<" keV"<<G4endl;
2444 
2445  return loss;
2446 }
2447 
2448 
2450 //
2451 // Returns resonance energy transfer in one collision
2452 
2454 {
2455  G4int iTransfer ;
2456 
2457  G4double energyTransfer, position;
2458 
2459  position = fIntegralResonance[1]*G4UniformRand();
2460 
2461  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2462  {
2463  if( position >= fIntegralResonance[iTransfer] ) break;
2464  }
2465  if(iTransfer > fSplineNumber) iTransfer--;
2466 
2467  energyTransfer = fSplineEnergy[iTransfer];
2468 
2469  if(iTransfer > 1)
2470  {
2471  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2472  }
2473  return energyTransfer;
2474 }
2475 
2476 
2478 //
2479 // Returns Rutherford energy transfer in one collision
2480 
2482 {
2483  G4int iTransfer ;
2484 
2485  G4double energyTransfer, position;
2486 
2487  position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand();
2488 
2489  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2490  {
2491  if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break;
2492  }
2493  if(iTransfer > fSplineNumber) iTransfer--;
2494 
2495  energyTransfer = fSplineEnergy[iTransfer];
2496 
2497  if(iTransfer > 1)
2498  {
2499  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2500  }
2501  return energyTransfer;
2502 }
2503 
2505 //
2506 
2507 void G4PAIxSection::CallError(G4int i, const G4String& methodName) const
2508 {
2509  G4String head = "G4PAIxSection::" + methodName + "()";
2511  ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
2512  G4Exception(head,"pai001",FatalException,ed);
2513 }
2514 
2516 //
2517 // Init array of Lorentz factors
2518 //
2519 
2521 
2522 const G4double G4PAIxSection::fLorentzFactor[112] = // fNumberOfGammas+1
2523 {
2524 0.0,
2525 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
2526 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
2527 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
2528 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
2529 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
2530 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
2531 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
2532 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
2533 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
2534 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
2535 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
2536 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
2537 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
2538 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
2539 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
2540 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
2541 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
2542 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
2543 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
2544 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
2545 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
2546 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
2547 1.065799e+05
2548 };
2549 
2551 //
2552 // The number of gamma for creation of spline (near ion-min , G ~ 4 )
2553 //
2554 
2555 const
2557 
2558 
2559 //
2560 // end of G4PAIxSection implementation file
2561 //
2563 
void IntegralPAIxSection()
void SplainPAI(G4double betaGammaSq)
G4double GetRutherfordEnergyTransfer()
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
void IntegralCerenkov()
static c2_factory< G4double > c2
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4double GetPlasmonEnergyTransfer()
static const double cm2
Definition: G4SIunits.hh:107
void IntegralPlasmon()
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void IntegralResonance()
G4double GetStepResonanceLoss(G4double step)
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
G4double SumOverBordPlasmon(G4int intervalNumber, G4double energy)
size_t GetIndex() const
Definition: G4Material.hh:260
G4double SumOverInterMM(G4int intervalNumber)
static const G4double betaBohr
const G4double pi
G4double GetZ() const
Definition: G4Element.hh:131
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:564
std::vector< G4Material * > G4MaterialTable
static const G4double cofBetaBohr
G4double GetDensity() const
Definition: G4Material.hh:178
long G4long
Definition: G4Types.hh:80
void NormShift(G4double betaGammaSq)
G4int GetMaxInterval() const
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
G4double SumOverBordResonance(G4int intervalNumber, G4double energy)
static G4int fNumberOfGammas
G4double SumOverInterPlasmon(G4int intervalNumber)
G4double GetSandiaMatTablePAI(G4int, G4int)
G4double a
Definition: TRTMaterials.hh:39
G4double SumOverInterval(G4int intervalNumber)
G4double RePartDielectricConst(G4double energy)
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:200
void ComputeLowEnergyCof()
int G4int
Definition: G4Types.hh:78
void CallError(G4int i, const G4String &methodName) const
static const G4int fMaxSplineSize
G4int SandiaMixing(G4int Z[], const G4double *fractionW, G4int el, G4int mi)
static const G4double betaBohr4
G4double GetResonanceEnergyTransfer()
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4double GetElectronDensity() const
Definition: G4Material.hh:215
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
G4bool GetLowerI1()
G4double GetPhotonRange(G4double energy)
G4double SumOverBordMM(G4int intervalNumber, G4double energy)
static const G4double betaBohr2
static const G4double c3
G4double GetCerenkovEnergyTransfer()
static const G4int fRefGammaNumber
static const G4double fError
static const G4double c1
G4double SumOverInterResonance(G4int intervalNumber)
G4int SandiaIntervals(G4int Z[], G4int el)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetStepPlasmonLoss(G4double step)
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
static const double eV
Definition: G4SIunits.hh:194
G4double GetStepEnergyLoss(G4double step)
G4double SumOverBordCerenkov(G4int intervalNumber, G4double energy)
G4double GetMMEnergyTransfer()
G4double SumOverBorder(G4int intervalNumber, G4double energy)
G4double energy(const ThreeVector &p, const G4double m)
static const double g
Definition: G4SIunits.hh:162
static const G4double fDelta
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
int position
Definition: filter.cc:7
G4double SumOverInterCerenkov(G4int intervalNumber)
#define DBL_MIN
Definition: templates.hh:75
G4double GetPhotoAbsorpCof(G4int i, G4int j) const
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
static const double keV
Definition: G4SIunits.hh:195
double G4double
Definition: G4Types.hh:76
G4double GetEnergyTransfer()
G4double GetStepMMLoss(G4double step)
G4double SumOverIntervaldEdx(G4int intervalNumber)
#define DBL_MAX
Definition: templates.hh:83
G4double GetStepCerenkovLoss(G4double step)
G4double SumOverBorderdEdx(G4int intervalNumber, G4double energy)
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)
static const G4double fLorentzFactor[112]
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
const G4Material * GetMaterial() const
G4double GetElectronRange(G4double energy)