Geant4  10.02.p02
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 91726 2015-08-03 15:41:36Z 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 #include "G4PAIxSection.hh"
50 
51 #include "globals.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4ios.hh"
55 #include "G4Poisson.hh"
56 #include "G4Material.hh"
57 #include "G4MaterialCutsCouple.hh"
58 #include "G4SandiaTable.hh"
59 
60 using namespace std;
61 
62 /* ******************************************************************
63 
64 // Init array of Lorentz factors
65 
66 const G4double G4PAIxSection::fLorentzFactor[22] =
67 {
68  0.0 , 1.1 , 1.2 , 1.3 , 1.5 , 1.8 , 2.0 ,
69  2.5 , 3.0 , 4.0 , 7.0 , 10.0 , 20.0 , 40.0 ,
70  70.0 , 100.0 , 300.0 , 600.0 , 1000.0 , 3000.0 ,
71  10000.0 , 50000.0
72 };
73 
74 const G4int G4PAIxSection::
75 fRefGammaNumber = 29; // The number of gamma for creation of
76  // spline (9)
77 
78 ***************************************************************** */
79 
80 // Local class constants
81 
82 const G4double G4PAIxSection::fDelta = 0.005; // 0.005 energy shift from interval border
83 const G4double G4PAIxSection::fError = 0.005; // 0.005 error in lin-log approximation
84 
85 const G4int G4PAIxSection::fMaxSplineSize = 1000; // Max size of output spline
86  // arrays
88 //
89 // Constructor
90 //
91 
93 {
94  fSandia = 0;
95  fMatSandiaMatrix = 0;
96  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
97  fIntervalNumber = fSplineNumber = 0;
98  fVerbose = 0;
99 
100  fSplineEnergy = G4DataVector(fMaxSplineSize,0.0);
101  fRePartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
102  fImPartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
103  fIntegralTerm = G4DataVector(fMaxSplineSize,0.0);
104  fDifPAIxSection = G4DataVector(fMaxSplineSize,0.0);
105  fdNdxCerenkov = G4DataVector(fMaxSplineSize,0.0);
106  fdNdxPlasmon = G4DataVector(fMaxSplineSize,0.0);
107  fdNdxMM = G4DataVector(fMaxSplineSize,0.0);
108  fdNdxResonance = G4DataVector(fMaxSplineSize,0.0);
109  fIntegralPAIxSection = G4DataVector(fMaxSplineSize,0.0);
110  fIntegralPAIdEdx = G4DataVector(fMaxSplineSize,0.0);
111  fIntegralCerenkov = G4DataVector(fMaxSplineSize,0.0);
112  fIntegralPlasmon = G4DataVector(fMaxSplineSize,0.0);
113  fIntegralMM = G4DataVector(fMaxSplineSize,0.0);
114  fIntegralResonance = G4DataVector(fMaxSplineSize,0.0);
115 
116  fMaterialIndex = 0;
117 
118  for( G4int i = 0; i < 500; ++i )
119  {
120  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
121  }
122 }
123 
125 //
126 // Constructor
127 //
128 
130 {
131  fDensity = matCC->GetMaterial()->GetDensity();
132  G4int matIndex = matCC->GetMaterial()->GetIndex();
133  fMaterialIndex = matIndex;
134 
135  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
136  fSandia = (*theMaterialTable)[matIndex]->GetSandiaTable();
137 
138  fVerbose = 0;
139 
140  G4int i, j;
141  fMatSandiaMatrix = new G4OrderedTable();
142 
143  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
144  {
145  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
146  }
147  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
148  {
149  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
150 
151  for(j = 1; j < 5; j++)
152  {
153  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
154  }
155  }
156  ComputeLowEnergyCof();
157  // fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0;
158 }
159 
161 
163  G4double maxEnergyTransfer)
164 {
165  fSandia = 0;
166  fMatSandiaMatrix = 0;
167  fVerbose = 0;
168  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
169  G4int i, j;
170 
171  fMaterialIndex = materialIndex;
172  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
173  fElectronDensity = (*theMaterialTable)[materialIndex]->
174  GetElectronDensity();
175  fIntervalNumber = (*theMaterialTable)[materialIndex]->
176  GetSandiaTable()->GetMatNbOfIntervals();
177  fIntervalNumber--;
178  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
179 
180  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
181  fA1 = G4DataVector(fIntervalNumber+2,0.0);
182  fA2 = G4DataVector(fIntervalNumber+2,0.0);
183  fA3 = G4DataVector(fIntervalNumber+2,0.0);
184  fA4 = G4DataVector(fIntervalNumber+2,0.0);
185 
186  for(i = 1; i <= fIntervalNumber; i++ )
187  {
188  if(((*theMaterialTable)[materialIndex]->
189  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
190  i > fIntervalNumber )
191  {
192  fEnergyInterval[i] = maxEnergyTransfer;
193  fIntervalNumber = i;
194  break;
195  }
196  fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
197  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
198  fA1[i] = (*theMaterialTable)[materialIndex]->
199  GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
200  fA2[i] = (*theMaterialTable)[materialIndex]->
201  GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
202  fA3[i] = (*theMaterialTable)[materialIndex]->
203  GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
204  fA4[i] = (*theMaterialTable)[materialIndex]->
205  GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
206  // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
207  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
208  }
209  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
210  {
211  fIntervalNumber++;
212  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
213  }
214 
215  // Now checking, if two borders are too close together
216 
217  for(i=1;i<fIntervalNumber;i++)
218  {
219  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
220  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
221  {
222  continue;
223  }
224  else
225  {
226  for(j=i;j<fIntervalNumber;j++)
227  {
228  fEnergyInterval[j] = fEnergyInterval[j+1];
229  fA1[j] = fA1[j+1];
230  fA2[j] = fA2[j+1];
231  fA3[j] = fA3[j+1];
232  fA4[j] = fA4[j+1];
233  }
234  fIntervalNumber--;
235  i--;
236  }
237  }
238 
239 
240  /* *********************************
241 
242  fSplineEnergy = new G4double[fMaxSplineSize];
243  fRePartDielectricConst = new G4double[fMaxSplineSize];
244  fImPartDielectricConst = new G4double[fMaxSplineSize];
245  fIntegralTerm = new G4double[fMaxSplineSize];
246  fDifPAIxSection = new G4double[fMaxSplineSize];
247  fIntegralPAIxSection = new G4double[fMaxSplineSize];
248 
249  for(i=0;i<fMaxSplineSize;i++)
250  {
251  fSplineEnergy[i] = 0.0;
252  fRePartDielectricConst[i] = 0.0;
253  fImPartDielectricConst[i] = 0.0;
254  fIntegralTerm[i] = 0.0;
255  fDifPAIxSection[i] = 0.0;
256  fIntegralPAIxSection[i] = 0.0;
257  }
258  ************************************************** */
259  ComputeLowEnergyCof();
260  InitPAI(); // create arrays allocated above
261  /*
262  delete[] fEnergyInterval;
263  delete[] fA1;
264  delete[] fA2;
265  delete[] fA3;
266  delete[] fA4;
267  */
268 }
269 
271 //
272 // Constructor called from G4PAIPhotonModel !!!
273 
275  G4double maxEnergyTransfer,
276  G4double betaGammaSq,
277  G4double** photoAbsCof,
278  G4int intNumber )
279 {
280  fSandia = 0;
281  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
282  fIntervalNumber = fSplineNumber = 0;
283  fVerbose = 0;
284 
285  fSplineEnergy = G4DataVector(500,0.0);
286  fRePartDielectricConst = G4DataVector(500,0.0);
287  fImPartDielectricConst = G4DataVector(500,0.0);
288  fIntegralTerm = G4DataVector(500,0.0);
289  fDifPAIxSection = G4DataVector(500,0.0);
290  fdNdxCerenkov = G4DataVector(500,0.0);
291  fdNdxPlasmon = G4DataVector(500,0.0);
292  fdNdxMM = G4DataVector(500,0.0);
293  fdNdxResonance = G4DataVector(500,0.0);
294  fIntegralPAIxSection = G4DataVector(500,0.0);
295  fIntegralPAIdEdx = G4DataVector(500,0.0);
296  fIntegralCerenkov = G4DataVector(500,0.0);
297  fIntegralPlasmon = G4DataVector(500,0.0);
298  fIntegralMM = G4DataVector(500,0.0);
299  fIntegralResonance = G4DataVector(500,0.0);
300 
301  for( G4int i = 0; i < 500; ++i )
302  {
303  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
304  }
305 
306  fSandia = 0;
307  fMatSandiaMatrix = 0;
308  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
309  G4int i, j;
310 
311  fMaterialIndex = materialIndex;
312  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
313  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
314 
315  fIntervalNumber = intNumber;
316  fIntervalNumber--;
317  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
318 
319  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
320  fA1 = G4DataVector(fIntervalNumber+2,0.0);
321  fA2 = G4DataVector(fIntervalNumber+2,0.0);
322  fA3 = G4DataVector(fIntervalNumber+2,0.0);
323  fA4 = G4DataVector(fIntervalNumber+2,0.0);
324 
325 
326  /*
327  fEnergyInterval = new G4double[fIntervalNumber+2];
328  fA1 = new G4double[fIntervalNumber+2];
329  fA2 = new G4double[fIntervalNumber+2];
330  fA3 = new G4double[fIntervalNumber+2];
331  fA4 = new G4double[fIntervalNumber+2];
332  */
333  for( i = 1; i <= fIntervalNumber; i++ )
334  {
335  if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
336  i > fIntervalNumber )
337  {
338  fEnergyInterval[i] = maxEnergyTransfer;
339  fIntervalNumber = i;
340  break;
341  }
342  fEnergyInterval[i] = photoAbsCof[i-1][0];
343  fA1[i] = photoAbsCof[i-1][1];
344  fA2[i] = photoAbsCof[i-1][2];
345  fA3[i] = photoAbsCof[i-1][3];
346  fA4[i] = photoAbsCof[i-1][4];
347  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
348  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
349  }
350  // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
351 
352  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
353  {
354  fIntervalNumber++;
355  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
356  }
357  // G4cout<<"after check of max energy transfer"<<G4endl;
358 
359  for( i = 1; i <= fIntervalNumber; i++ )
360  {
361  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
362  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
363  }
364  // Now checking, if two borders are too close together
365 
366  for( i = 1; i < fIntervalNumber; i++ )
367  {
368  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
369  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
370  {
371  continue;
372  }
373  else
374  {
375  for(j=i;j<fIntervalNumber;j++)
376  {
377  fEnergyInterval[j] = fEnergyInterval[j+1];
378  fA1[j] = fA1[j+1];
379  fA2[j] = fA2[j+1];
380  fA3[j] = fA3[j+1];
381  fA4[j] = fA4[j+1];
382  }
383  fIntervalNumber--;
384  i--;
385  }
386  }
387  // G4cout<<"after check of close borders"<<G4endl;
388 
389  for( i = 1; i <= fIntervalNumber; i++ )
390  {
391  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
392  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
393  }
394 
395  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
396 
397  ComputeLowEnergyCof();
398  G4double betaGammaSqRef =
399  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
400 
401  NormShift(betaGammaSqRef);
402  SplainPAI(betaGammaSqRef);
403 
404  // Preparation of integral PAI cross section for input betaGammaSq
405 
406  for(i = 1; i <= fSplineNumber; i++)
407  {
408  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
409  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
410  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
411  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
412  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
413 
414  // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
415  // <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
416  }
417  IntegralCerenkov();
418  IntegralMM();
419  IntegralPlasmon();
420  IntegralResonance();
421  IntegralPAIxSection();
422  /*
423  delete[] fEnergyInterval;
424  delete[] fA1;
425  delete[] fA2;
426  delete[] fA3;
427  delete[] fA4;
428  */
429 }
430 
432 //
433 // Test Constructor with beta*gamma square value
434 
436  G4double maxEnergyTransfer,
437  G4double betaGammaSq )
438 {
439  fSandia = 0;
440  fMatSandiaMatrix = 0;
441  fVerbose = 0;
442  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
443 
444  G4int i, j, numberOfElements;
445 
446  fMaterialIndex = materialIndex;
447  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
448  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
449  numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements();
450 
451  G4int* thisMaterialZ = new G4int[numberOfElements];
452 
453  for( i = 0; i < numberOfElements; i++ )
454  {
455  thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
456  GetElement(i)->GetZ();
457  }
458  // fSandia = new G4SandiaTable(materialIndex);
459  fSandia = (*theMaterialTable)[materialIndex]->GetSandiaTable();
460  G4SandiaTable thisMaterialSandiaTable(materialIndex);
461  fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals(thisMaterialZ,
462  numberOfElements);
463  fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
464  ( thisMaterialZ ,
465  (*theMaterialTable)[materialIndex]->GetFractionVector() ,
466  numberOfElements,fIntervalNumber);
467 
468  fIntervalNumber--;
469 
470  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
471  fA1 = G4DataVector(fIntervalNumber+2,0.0);
472  fA2 = G4DataVector(fIntervalNumber+2,0.0);
473  fA3 = G4DataVector(fIntervalNumber+2,0.0);
474  fA4 = G4DataVector(fIntervalNumber+2,0.0);
475 
476  /*
477  fEnergyInterval = new G4double[fIntervalNumber+2];
478  fA1 = new G4double[fIntervalNumber+2];
479  fA2 = new G4double[fIntervalNumber+2];
480  fA3 = new G4double[fIntervalNumber+2];
481  fA4 = new G4double[fIntervalNumber+2];
482  */
483  for( i = 1; i <= fIntervalNumber; i++ )
484  {
485  if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
486  i > fIntervalNumber)
487  {
488  fEnergyInterval[i] = maxEnergyTransfer;
489  fIntervalNumber = i;
490  break;
491  }
492  fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0);
493  fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity;
494  fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity;
495  fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity;
496  fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity;
497 
498  }
499  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
500  {
501  fIntervalNumber++;
502  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
503  fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
504  fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
505  fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
506  fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
507  }
508  for(i=1;i<=fIntervalNumber;i++)
509  {
510  // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
511  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
512  }
513  // Now checking, if two borders are too close together
514 
515  for( i = 1; i < fIntervalNumber; i++ )
516  {
517  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
518  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
519  {
520  continue;
521  }
522  else
523  {
524  for( j = i; j < fIntervalNumber; j++ )
525  {
526  fEnergyInterval[j] = fEnergyInterval[j+1];
527  fA1[j] = fA1[j+1];
528  fA2[j] = fA2[j+1];
529  fA3[j] = fA3[j+1];
530  fA4[j] = fA4[j+1];
531  }
532  fIntervalNumber--;
533  i--;
534  }
535  }
536 
537  /* *********************************
538  fSplineEnergy = new G4double[fMaxSplineSize];
539  fRePartDielectricConst = new G4double[fMaxSplineSize];
540  fImPartDielectricConst = new G4double[fMaxSplineSize];
541  fIntegralTerm = new G4double[fMaxSplineSize];
542  fDifPAIxSection = new G4double[fMaxSplineSize];
543  fIntegralPAIxSection = new G4double[fMaxSplineSize];
544 
545  for(i=0;i<fMaxSplineSize;i++)
546  {
547  fSplineEnergy[i] = 0.0;
548  fRePartDielectricConst[i] = 0.0;
549  fImPartDielectricConst[i] = 0.0;
550  fIntegralTerm[i] = 0.0;
551  fDifPAIxSection[i] = 0.0;
552  fIntegralPAIxSection[i] = 0.0;
553  }
554  */
555 
556  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
557 
558  ComputeLowEnergyCof();
559  G4double betaGammaSqRef =
560  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
561 
562  NormShift(betaGammaSqRef);
563  SplainPAI(betaGammaSqRef);
564 
565  // Preparation of integral PAI cross section for input betaGammaSq
566 
567  for(i = 1; i <= fSplineNumber; i++)
568  {
569  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
570  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
571  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
572  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
573  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
574  }
575  IntegralPAIxSection();
576  IntegralCerenkov();
577  IntegralMM();
578  IntegralPlasmon();
579  IntegralResonance();
580 }
581 
583 //
584 // Destructor
585 
587 {
588  /* ************************
589  delete[] fSplineEnergy ;
590  delete[] fRePartDielectricConst;
591  delete[] fImPartDielectricConst;
592  delete[] fIntegralTerm ;
593  delete[] fDifPAIxSection ;
594  delete[] fIntegralPAIxSection ;
595  */
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  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1045  } // close 'while'
1046 
1047 } // end of SplainPAI
1048 
1049 
1051 //
1052 // Integration over electrons that could be considered
1053 // quasi-free at energy transfer of interest
1054 
1056  G4double x1,
1057  G4double x2 )
1058 {
1059  G4double c1, c2, c3;
1060  // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
1061  c1 = (x2 - x1)/x1/x2;
1062  c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2;
1063  c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
1064  // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
1065 
1066  return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
1067 
1068 } // end of RutherfordIntegral
1069 
1070 
1072 //
1073 // Imaginary part of dielectric constant
1074 // (G4int k - interval number, G4double en1 - energy point)
1075 
1077  G4double energy1 )
1078 {
1079  G4double energy2,energy3,energy4,result;
1080 
1081  energy2 = energy1*energy1;
1082  energy3 = energy2*energy1;
1083  energy4 = energy3*energy1;
1084 
1085  result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
1086  result *=hbarc/energy1;
1087 
1088  return result;
1089 
1090 } // end of ImPartDielectricConst
1091 
1093 //
1094 // Returns lambda of photon with energy1 in current material
1095 
1097 {
1098  G4int i;
1099  G4double energy2, energy3, energy4, result, lambda;
1100 
1101  energy2 = energy1*energy1;
1102  energy3 = energy2*energy1;
1103  energy4 = energy3*energy1;
1104 
1105  // G4double* SandiaCof = fSandia->GetSandiaCofForMaterialPAI(energy1);
1106  // result = SandiaCof[0]/energy1+SandiaCof[1]/energy2+SandiaCof[2]/energy3+SandiaCof[3]/energy4;
1107  // result *= fDensity;
1108 
1109  for( i = 1; i <= fIntervalNumber; i++ )
1110  {
1111  if( energy1 < fEnergyInterval[i]) break;
1112  }
1113  i--;
1114  if(i == 0) i = 1;
1115 
1116  result = fA1[i]/energy1+fA2[i]/energy2+fA3[i]/energy3+fA4[i]/energy4;
1117 
1118  if( result > DBL_MIN ) lambda = 1./result;
1119  else lambda = DBL_MAX;
1120 
1121  return lambda;
1122 }
1123 
1125 //
1126 // Return lambda of electron with energy1 in current material
1127 // parametrisation from NIM A554(2005)474-493
1128 
1130 {
1131  G4double range;
1132  /*
1133  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1134 
1135  G4double Z = (*theMaterialTable)[fMaterialIndex]->GetIonisation()->GetZeffective();
1136  G4double A = (*theMaterialTable)[fMaterialIndex]->GetA();
1137 
1138  energy /= keV; // energy in keV in parametrised formula
1139 
1140  if (energy < 10.)
1141  {
1142  range = 3.872e-3*A/Z;
1143  range *= pow(energy,1.492);
1144  }
1145  else
1146  {
1147  range = 6.97e-3*pow(energy,1.6);
1148  }
1149  */
1150  // Blum&Rolandi Particle Detection with Drift Chambers, p. 7
1151 
1152  G4double cofA = 5.37e-4*g/cm2/keV;
1153  G4double cofB = 0.9815;
1154  G4double cofC = 3.123e-3/keV;
1155  // energy /= keV;
1156 
1157  range = cofA*energy*( 1 - cofB/(1 + cofC*energy) );
1158 
1159  // range *= g/cm2;
1160  range /= fDensity;
1161 
1162  return range;
1163 }
1164 
1166 //
1167 // Real part of dielectric constant minus unit: epsilon_1 - 1
1168 // (G4double enb - energy point)
1169 //
1170 
1172 {
1173  G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
1174  c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
1175 
1176  x0 = enb;
1177  result = 0;
1178 
1179  for(G4int i=1;i<=fIntervalNumber-1;i++)
1180  {
1181  x1 = fEnergyInterval[i];
1182  x2 = fEnergyInterval[i+1];
1183  xx1 = x1 - x0;
1184  xx2 = x2 - x0;
1185  xx12 = xx2/xx1;
1186 
1187  if(xx12<0)
1188  {
1189  xx12 = -xx12;
1190  }
1191  xln1 = log(x2/x1);
1192  xln2 = log(xx12);
1193  xln3 = log((x2 + x0)/(x1 + x0));
1194  x02 = x0*x0;
1195  x03 = x02*x0;
1196  x04 = x03*x0;
1197  x05 = x04*x0;
1198  c1 = (x2 - x1)/x1/x2;
1199  c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2;
1200  c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
1201 
1202  result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
1203  result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
1204  result -= fA3[i]*c2/2/x02;
1205  result -= fA4[i]*c3/3/x02;
1206 
1207  cof1 = fA1[i]/x02 + fA3[i]/x04;
1208  cof2 = fA2[i]/x03 + fA4[i]/x05;
1209 
1210  result += 0.5*(cof1 +cof2)*xln2;
1211  result += 0.5*(cof1 - cof2)*xln3;
1212  }
1213  result *= 2*hbarc/pi;
1214 
1215  return result;
1216 
1217 } // end of RePartDielectricConst
1218 
1220 //
1221 // PAI differential cross-section in terms of
1222 // simplified Allison's equation
1223 //
1224 
1226  G4double betaGammaSq )
1227 {
1228  G4double cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
1229 
1230  G4double betaBohr = fine_structure_const;
1231  // G4double betaBohr2 = fine_structure_const*fine_structure_const;
1232  // G4double betaBohr3 = betaBohr*betaBohr2; // *4.0;
1233 
1234  G4double be2 = betaGammaSq/(1 + betaGammaSq);
1235  G4double beta = sqrt(be2);
1236  // G4double be3 = beta*be2;
1237 
1238  cof = 1.;
1239  x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
1240 
1241  if( betaGammaSq < 0.01 ) x2 = log(be2);
1242  else
1243  {
1244  x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1245  (1/betaGammaSq - fRePartDielectricConst[i]) +
1246  fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
1247  }
1248  if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
1249  {
1250  x6 = 0.;
1251  }
1252  else
1253  {
1254  x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
1255  x5 = -1 - fRePartDielectricConst[i] +
1256  be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1257  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1258 
1259  x7 = atan2(fImPartDielectricConst[i],x3);
1260  x6 = x5 * x7;
1261  }
1262  // if(fImPartDielectricConst[i] == 0) x6 = 0.;
1263 
1264  x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
1265 
1266  // if( x4 < 0.0 ) x4 = 0.0;
1267 
1268  x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1269  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1270 
1271  result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
1272 
1273  if( result < 1.0e-8 ) result = 1.0e-8;
1274 
1275  result *= fine_structure_const/be2/pi;
1276 
1277  // low energy correction
1278 
1279  G4double lowCof = fLowEnergyCof; // 6.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)?
1280 
1281  result *= (1 - exp(-beta/betaBohr/lowCof));
1282 
1283 
1284  // result *= (1 - exp(-be2/betaBohr2/lowCof));
1285 
1286  // result *= (1 - exp(-be3/betaBohr3/lowCof)); // ~ be for be<<betaBohr
1287 
1288  // result *= (1 - exp(-be4/betaBohr4/lowCof));
1289 
1290  if(fDensity >= 0.1)
1291  {
1292  result /= x8;
1293  }
1294  return result;
1295 
1296 } // end of DifPAIxSection
1297 
1299 //
1300 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
1301 
1303  G4double betaGammaSq )
1304 {
1305  G4double logarithm, x3, x5, argument, modul2, dNdxC;
1307 
1308  cofBetaBohr = 4.0;
1309  betaBohr2 = fine_structure_const*fine_structure_const;
1310  G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1311 
1312  be2 = betaGammaSq/(1 + betaGammaSq);
1313  G4double be4 = be2*be2;
1314 
1315  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1316  else
1317  {
1318  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1319  (1/betaGammaSq - fRePartDielectricConst[i]) +
1320  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1321  logarithm += log(1+1.0/betaGammaSq);
1322  }
1323 
1324  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1325  {
1326  argument = 0.0;
1327  }
1328  else
1329  {
1330  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1331  x5 = -1.0 - fRePartDielectricConst[i] +
1332  be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1333  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1334  if( x3 == 0.0 ) argument = 0.5*pi;
1335  else argument = atan2(fImPartDielectricConst[i],x3);
1336  argument *= x5 ;
1337  }
1338  dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
1339 
1340  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1341 
1342  dNdxC *= fine_structure_const/be2/pi;
1343 
1344  dNdxC *= (1-exp(-be4/betaBohr4));
1345 
1346  if(fDensity >= 0.1)
1347  {
1348  modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1349  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1350  dNdxC /= modul2;
1351  }
1352  return dNdxC;
1353 
1354 } // end of PAIdNdxCerenkov
1355 
1357 //
1358 // Calculation od dN/dx of collisions of MM with creation of Cerenkov pseudo-photons
1359 
1361  G4double betaGammaSq )
1362 {
1363  G4double logarithm, x3, x5, argument, dNdxC;
1365 
1366  cofBetaBohr = 4.0;
1367  betaBohr2 = fine_structure_const*fine_structure_const;
1368  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1369 
1370  be2 = betaGammaSq/(1 + betaGammaSq);
1371  be4 = be2*be2;
1372 
1373  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1374  else
1375  {
1376  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1377  (1/betaGammaSq - fRePartDielectricConst[i]) +
1378  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1379  logarithm += log(1+1.0/betaGammaSq);
1380  }
1381 
1382  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1383  {
1384  argument = 0.0;
1385  }
1386  else
1387  {
1388  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1389  x5 = be2*( 1.0 + fRePartDielectricConst[i] ) - 1.0;
1390  if( x3 == 0.0 ) argument = 0.5*pi;
1391  else argument = atan2(fImPartDielectricConst[i],x3);
1392  argument *= x5 ;
1393  }
1394  dNdxC = ( logarithm*fImPartDielectricConst[i]*be2 + argument )/hbarc;
1395 
1396  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1397 
1398  dNdxC *= fine_structure_const/be2/pi;
1399 
1400  dNdxC *= (1-exp(-be4/betaBohr4));
1401  return dNdxC;
1402 
1403 } // end of PAIdNdxMM
1404 
1406 //
1407 // Calculation od dN/dx of collisions with creation of longitudinal EM
1408 // excitations (plasmons, delta-electrons)
1409 
1411  G4double betaGammaSq )
1412 {
1413  G4double resonance, modul2, dNdxP, cof = 1.;
1414  G4double be2, betaBohr;
1415 
1416  betaBohr = fine_structure_const;
1417  be2 = betaGammaSq/(1 + betaGammaSq);
1418 
1419  G4double beta = sqrt(be2);
1420 
1421  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1422  resonance *= fImPartDielectricConst[i]/hbarc;
1423 
1424 
1425  dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
1426 
1427  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1428 
1429  dNdxP *= fine_structure_const/be2/pi;
1430 
1431  dNdxP *= (1 - exp(-beta/betaBohr/fLowEnergyCof));
1432 
1433  // dNdxP *= (1-exp(-be4/betaBohr4));
1434 
1435  if( fDensity >= 0.1 )
1436  {
1437  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1438  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1439  dNdxP /= modul2;
1440  }
1441  return dNdxP;
1442 
1443 } // end of PAIdNdxPlasmon
1444 
1446 //
1447 // Calculation od dN/dx of collisions with creation of longitudinal EM
1448 // resonance excitations (plasmons, delta-electrons)
1449 
1451  G4double betaGammaSq )
1452 {
1453  G4double resonance, modul2, dNdxP;
1455 
1456  cofBetaBohr = 4.0;
1457  betaBohr2 = fine_structure_const*fine_structure_const;
1458  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1459 
1460  be2 = betaGammaSq/(1 + betaGammaSq);
1461  be4 = be2*be2;
1462 
1463  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1464  resonance *= fImPartDielectricConst[i]/hbarc;
1465 
1466 
1467  dNdxP = resonance;
1468 
1469  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1470 
1471  dNdxP *= fine_structure_const/be2/pi;
1472  dNdxP *= (1-exp(-be4/betaBohr4));
1473 
1474  if( fDensity >= 0.1 )
1475  {
1476  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1477  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1478  dNdxP /= modul2;
1479  }
1480  return dNdxP;
1481 
1482 } // end of PAIdNdxResonance
1483 
1485 //
1486 // Calculation of the PAI integral cross-section
1487 // fIntegralPAIxSection[1] = specific primary ionisation, 1/cm
1488 // and fIntegralPAIxSection[0] = mean energy loss per cm in keV/cm
1489 
1491 {
1492  fIntegralPAIxSection[fSplineNumber] = 0;
1493  fIntegralPAIdEdx[fSplineNumber] = 0;
1494  fIntegralPAIxSection[0] = 0;
1495  G4int i, k = fIntervalNumber -1;
1496 
1497  for( i = fSplineNumber-1; i >= 1; i--)
1498  {
1499  if(fSplineEnergy[i] >= fEnergyInterval[k])
1500  {
1501  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i);
1502  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
1503  }
1504  else
1505  {
1506  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] +
1507  SumOverBorder(i+1,fEnergyInterval[k]);
1508  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
1509  SumOverBorderdEdx(i+1,fEnergyInterval[k]);
1510  k--;
1511  }
1512  if(fVerbose>0) G4cout<<"i = "<<i<<"; k = "<<k<<"; intPAIxsc[i] = "<<fIntegralPAIxSection[i]<<G4endl;
1513  }
1514 } // end of IntegralPAIxSection
1515 
1517 //
1518 // Calculation of the PAI Cerenkov integral cross-section
1519 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
1520 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
1521 
1523 {
1524  G4int i, k;
1525  fIntegralCerenkov[fSplineNumber] = 0;
1526  fIntegralCerenkov[0] = 0;
1527  k = fIntervalNumber -1;
1528 
1529  for( i = fSplineNumber-1; i >= 1; i-- )
1530  {
1531  if(fSplineEnergy[i] >= fEnergyInterval[k])
1532  {
1533  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
1534  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1535  }
1536  else
1537  {
1538  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
1539  SumOverBordCerenkov(i+1,fEnergyInterval[k]);
1540  k--;
1541  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1542  }
1543  }
1544 
1545 } // end of IntegralCerenkov
1546 
1548 //
1549 // Calculation of the PAI MM-Cerenkov integral cross-section
1550 // fIntegralMM[1] = specific MM-Cerenkov ionisation, 1/cm
1551 // and fIntegralMM[0] = mean MM-Cerenkov loss per cm in keV/cm
1552 
1554 {
1555  G4int i, k;
1556  fIntegralMM[fSplineNumber] = 0;
1557  fIntegralMM[0] = 0;
1558  k = fIntervalNumber -1;
1559 
1560  for( i = fSplineNumber-1; i >= 1; i-- )
1561  {
1562  if(fSplineEnergy[i] >= fEnergyInterval[k])
1563  {
1564  fIntegralMM[i] = fIntegralMM[i+1] + SumOverInterMM(i);
1565  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1566  }
1567  else
1568  {
1569  fIntegralMM[i] = fIntegralMM[i+1] +
1570  SumOverBordMM(i+1,fEnergyInterval[k]);
1571  k--;
1572  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1573  }
1574  }
1575 
1576 } // end of IntegralMM
1577 
1579 //
1580 // Calculation of the PAI Plasmon integral cross-section
1581 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
1582 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
1583 
1585 {
1586  fIntegralPlasmon[fSplineNumber] = 0;
1587  fIntegralPlasmon[0] = 0;
1588  G4int k = fIntervalNumber -1;
1589  for(G4int i=fSplineNumber-1;i>=1;i--)
1590  {
1591  if(fSplineEnergy[i] >= fEnergyInterval[k])
1592  {
1593  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
1594  }
1595  else
1596  {
1597  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
1598  SumOverBordPlasmon(i+1,fEnergyInterval[k]);
1599  k--;
1600  }
1601  }
1602 
1603 } // end of IntegralPlasmon
1604 
1606 //
1607 // Calculation of the PAI resonance integral cross-section
1608 // fIntegralResonance[1] = resonance primary ionisation, 1/cm
1609 // and fIntegralResonance[0] = mean resonance loss per cm in keV/cm
1610 
1612 {
1613  fIntegralResonance[fSplineNumber] = 0;
1614  fIntegralResonance[0] = 0;
1615  G4int k = fIntervalNumber -1;
1616  for(G4int i=fSplineNumber-1;i>=1;i--)
1617  {
1618  if(fSplineEnergy[i] >= fEnergyInterval[k])
1619  {
1620  fIntegralResonance[i] = fIntegralResonance[i+1] + SumOverInterResonance(i);
1621  }
1622  else
1623  {
1624  fIntegralResonance[i] = fIntegralResonance[i+1] +
1625  SumOverBordResonance(i+1,fEnergyInterval[k]);
1626  k--;
1627  }
1628  }
1629 
1630 } // end of IntegralResonance
1631 
1633 //
1634 // Calculation the PAI integral cross-section inside
1635 // of interval of continuous values of photo-ionisation
1636 // cross-section. Parameter 'i' is the number of interval.
1637 
1639 {
1640  G4double x0,x1,y0,yy1,a,b,c,result;
1641 
1642  x0 = fSplineEnergy[i];
1643  x1 = fSplineEnergy[i+1];
1644  if(fVerbose>0) G4cout<<"SumOverInterval i= " << i << " x0 = "<<x0<<"; x1 = "<<x1<<G4endl;
1645 
1646  if( x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1647 
1648  y0 = fDifPAIxSection[i];
1649  yy1 = fDifPAIxSection[i+1];
1650 
1651  if(fVerbose>0) G4cout<<"x0 = "<<x0<<"; x1 = "<<x1<<", y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1652 
1653  c = x1/x0;
1654  a = log10(yy1/y0)/log10(c);
1655 
1656  if(fVerbose>0) G4cout<<"SumOverInterval, a = "<<a<<"; c = "<<c<<G4endl;
1657 
1658  // b = log10(y0) - a*log10(x0);
1659  b = y0/pow(x0,a);
1660  a += 1.;
1661  if( std::fabs(a) < 1.e-6 )
1662  {
1663  result = b*log(x1/x0);
1664  }
1665  else
1666  {
1667  result = y0*(x1*pow(c,a-1) - x0)/a;
1668  }
1669  a += 1.;
1670  if( std::fabs(a) < 1.e-6 )
1671  {
1672  fIntegralPAIxSection[0] += b*log(x1/x0);
1673  }
1674  else
1675  {
1676  fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1677  }
1678  if(fVerbose>0) G4cout<<"SumOverInterval, result = "<<result<<G4endl;
1679  return result;
1680 
1681 } // end of SumOverInterval
1682 
1684 
1686 {
1687  G4double x0,x1,y0,yy1,a,b,c,result;
1688 
1689  x0 = fSplineEnergy[i];
1690  x1 = fSplineEnergy[i+1];
1691 
1692  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1693 
1694  y0 = fDifPAIxSection[i];
1695  yy1 = fDifPAIxSection[i+1];
1696  c = x1/x0;
1697  a = log10(yy1/y0)/log10(c);
1698  // b = log10(y0) - a*log10(x0);
1699  b = y0/pow(x0,a);
1700  a += 2;
1701  if(a == 0)
1702  {
1703  result = b*log(x1/x0);
1704  }
1705  else
1706  {
1707  result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1708  }
1709  return result;
1710 
1711 } // end of SumOverInterval
1712 
1714 //
1715 // Calculation the PAI Cerenkov integral cross-section inside
1716 // of interval of continuous values of photo-ionisation Cerenkov
1717 // cross-section. Parameter 'i' is the number of interval.
1718 
1720 {
1721  G4double x0,x1,y0,yy1,a,b,c,result;
1722 
1723  x0 = fSplineEnergy[i];
1724  x1 = fSplineEnergy[i+1];
1725 
1726  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1727 
1728  y0 = fdNdxCerenkov[i];
1729  yy1 = fdNdxCerenkov[i+1];
1730  // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1731  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1732 
1733  c = x1/x0;
1734  a = log10(yy1/y0)/log10(c);
1735  b = y0/pow(x0,a);
1736 
1737  a += 1.0;
1738  if(a == 0) result = b*log(c);
1739  else result = y0*(x1*pow(c,a-1) - x0)/a;
1740  a += 1.0;
1741 
1742  if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
1743  else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1744  // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1745  return result;
1746 
1747 } // end of SumOverInterCerenkov
1748 
1750 //
1751 // Calculation the PAI MM-Cerenkov integral cross-section inside
1752 // of interval of continuous values of photo-ionisation Cerenkov
1753 // cross-section. Parameter 'i' is the number of interval.
1754 
1756 {
1757  G4double x0,x1,y0,yy1,a,b,c,result;
1758 
1759  x0 = fSplineEnergy[i];
1760  x1 = fSplineEnergy[i+1];
1761 
1762  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1763 
1764  y0 = fdNdxMM[i];
1765  yy1 = fdNdxMM[i+1];
1766  //G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1767  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1768 
1769  c = x1/x0;
1770  //G4cout<<" c = "<<c<< " yy1/y0= " << yy1/y0 <<G4endl;
1771  a = log10(yy1/y0)/log10(c);
1772  if(a > 10.0) return 0.;
1773  b = y0/pow(x0,a);
1774 
1775  a += 1.0;
1776  if(a == 0) result = b*log(c);
1777  else result = y0*(x1*pow(c,a-1) - x0)/a;
1778  a += 1.0;
1779 
1780  if( a == 0 ) fIntegralMM[0] += b*log(c);
1781  else fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1782  //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1783  return result;
1784 
1785 } // end of SumOverInterMM
1786 
1788 //
1789 // Calculation the PAI Plasmon integral cross-section inside
1790 // of interval of continuous values of photo-ionisation Plasmon
1791 // cross-section. Parameter 'i' is the number of interval.
1792 
1794 {
1795  G4double x0,x1,y0,yy1,a,b,c,result;
1796 
1797  x0 = fSplineEnergy[i];
1798  x1 = fSplineEnergy[i+1];
1799 
1800  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1801 
1802  y0 = fdNdxPlasmon[i];
1803  yy1 = fdNdxPlasmon[i+1];
1804  c =x1/x0;
1805  a = log10(yy1/y0)/log10(c);
1806  if(a > 10.0) return 0.;
1807  // b = log10(y0) - a*log10(x0);
1808  b = y0/pow(x0,a);
1809 
1810  a += 1.0;
1811  if(a == 0) result = b*log(x1/x0);
1812  else result = y0*(x1*pow(c,a-1) - x0)/a;
1813  a += 1.0;
1814 
1815  if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
1816  else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1817 
1818  return result;
1819 
1820 } // end of SumOverInterPlasmon
1821 
1823 //
1824 // Calculation the PAI resonance integral cross-section inside
1825 // of interval of continuous values of photo-ionisation resonance
1826 // cross-section. Parameter 'i' is the number of interval.
1827 
1829 {
1830  G4double x0,x1,y0,yy1,a,b,c,result;
1831 
1832  x0 = fSplineEnergy[i];
1833  x1 = fSplineEnergy[i+1];
1834 
1835  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1836 
1837  y0 = fdNdxResonance[i];
1838  yy1 = fdNdxResonance[i+1];
1839  c =x1/x0;
1840  a = log10(yy1/y0)/log10(c);
1841  if(a > 10.0) return 0.;
1842  // b = log10(y0) - a*log10(x0);
1843  b = y0/pow(x0,a);
1844 
1845  a += 1.0;
1846  if(a == 0) result = b*log(x1/x0);
1847  else result = y0*(x1*pow(c,a-1) - x0)/a;
1848  a += 1.0;
1849 
1850  if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0);
1851  else fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1852 
1853  return result;
1854 
1855 } // end of SumOverInterResonance
1856 
1858 //
1859 // Integration of PAI cross-section for the case of
1860 // passing across border between intervals
1861 
1863  G4double en0 )
1864 {
1865  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1866 
1867  e0 = en0;
1868  x0 = fSplineEnergy[i];
1869  x1 = fSplineEnergy[i+1];
1870  y0 = fDifPAIxSection[i];
1871  yy1 = fDifPAIxSection[i+1];
1872 
1873  //c = x1/x0;
1874  d = e0/x0;
1875  a = log10(yy1/y0)/log10(x1/x0);
1876  if(a > 10.0) return 0.;
1877 
1878  if(fVerbose>0) G4cout<<"SumOverBorder, a = "<<a<<G4endl;
1879 
1880  // b0 = log10(y0) - a*log10(x0);
1881  b = y0/pow(x0,a); // pow(10.,b);
1882 
1883  a += 1.;
1884  if( std::fabs(a) < 1.e-6 )
1885  {
1886  result = b*log(x0/e0);
1887  }
1888  else
1889  {
1890  result = y0*(x0 - e0*pow(d,a-1))/a;
1891  }
1892  a += 1.;
1893  if( std::fabs(a) < 1.e-6 )
1894  {
1895  fIntegralPAIxSection[0] += b*log(x0/e0);
1896  }
1897  else
1898  {
1899  fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1900  }
1901  x0 = fSplineEnergy[i - 1];
1902  x1 = fSplineEnergy[i - 2];
1903  y0 = fDifPAIxSection[i - 1];
1904  yy1 = fDifPAIxSection[i - 2];
1905 
1906  //c = x1/x0;
1907  d = e0/x0;
1908  a = log10(yy1/y0)/log10(x1/x0);
1909  // b0 = log10(y0) - a*log10(x0);
1910  b = y0/pow(x0,a);
1911  a += 1.;
1912  if( std::fabs(a) < 1.e-6 )
1913  {
1914  result += b*log(e0/x0);
1915  }
1916  else
1917  {
1918  result += y0*(e0*pow(d,a-1) - x0)/a;
1919  }
1920  a += 1.;
1921  if( std::fabs(a) < 1.e-6 )
1922  {
1923  fIntegralPAIxSection[0] += b*log(e0/x0);
1924  }
1925  else
1926  {
1927  fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1928  }
1929  return result;
1930 
1931 }
1932 
1934 
1936  G4double en0 )
1937 {
1938  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1939 
1940  e0 = en0;
1941  x0 = fSplineEnergy[i];
1942  x1 = fSplineEnergy[i+1];
1943  y0 = fDifPAIxSection[i];
1944  yy1 = fDifPAIxSection[i+1];
1945 
1946  //c = x1/x0;
1947  d = e0/x0;
1948  a = log10(yy1/y0)/log10(x1/x0);
1949  if(a > 10.0) return 0.;
1950  // b0 = log10(y0) - a*log10(x0);
1951  b = y0/pow(x0,a); // pow(10.,b);
1952 
1953  a += 2;
1954  if(a == 0)
1955  {
1956  result = b*log(x0/e0);
1957  }
1958  else
1959  {
1960  result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1961  }
1962  x0 = fSplineEnergy[i - 1];
1963  x1 = fSplineEnergy[i - 2];
1964  y0 = fDifPAIxSection[i - 1];
1965  yy1 = fDifPAIxSection[i - 2];
1966 
1967  // c = x1/x0;
1968  d = e0/x0;
1969  a = log10(yy1/y0)/log10(x1/x0);
1970  // b0 = log10(y0) - a*log10(x0);
1971  b = y0/pow(x0,a);
1972  a += 2;
1973  if(a == 0)
1974  {
1975  result += b*log(e0/x0);
1976  }
1977  else
1978  {
1979  result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1980  }
1981  return result;
1982 
1983 }
1984 
1986 //
1987 // Integration of Cerenkov cross-section for the case of
1988 // passing across border between intervals
1989 
1991  G4double en0 )
1992 {
1993  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1994 
1995  e0 = en0;
1996  x0 = fSplineEnergy[i];
1997  x1 = fSplineEnergy[i+1];
1998  y0 = fdNdxCerenkov[i];
1999  yy1 = fdNdxCerenkov[i+1];
2000 
2001  // G4cout<<G4endl;
2002  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2003  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2004  c = x1/x0;
2005  d = e0/x0;
2006  a = log10(yy1/y0)/log10(c);
2007  if(a > 10.0) return 0.;
2008  // b0 = log10(y0) - a*log10(x0);
2009  b = y0/pow(x0,a); // pow(10.,b0);
2010 
2011  a += 1.0;
2012  if( a == 0 ) result = b*log(x0/e0);
2013  else result = y0*(x0 - e0*pow(d,a-1))/a;
2014  a += 1.0;
2015 
2016  if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
2017  else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2018 
2019 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2020 
2021  x0 = fSplineEnergy[i - 1];
2022  x1 = fSplineEnergy[i - 2];
2023  y0 = fdNdxCerenkov[i - 1];
2024  yy1 = fdNdxCerenkov[i - 2];
2025 
2026  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2027  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2028 
2029  c = x1/x0;
2030  d = e0/x0;
2031  a = log10(yy1/y0)/log10(x1/x0);
2032  // b0 = log10(y0) - a*log10(x0);
2033  b = y0/pow(x0,a); // pow(10.,b0);
2034 
2035  a += 1.0;
2036  if( a == 0 ) result += b*log(e0/x0);
2037  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2038  a += 1.0;
2039 
2040  if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
2041  else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2042 
2043  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2044  // <<b<<"; result = "<<result<<G4endl;
2045 
2046  return result;
2047 
2048 }
2049 
2051 //
2052 // Integration of MM-Cerenkov cross-section for the case of
2053 // passing across border between intervals
2054 
2056  G4double en0 )
2057 {
2058  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
2059 
2060  e0 = en0;
2061  x0 = fSplineEnergy[i];
2062  x1 = fSplineEnergy[i+1];
2063  y0 = fdNdxMM[i];
2064  yy1 = fdNdxMM[i+1];
2065 
2066  // G4cout<<G4endl;
2067  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2068  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2069  c = x1/x0;
2070  d = e0/x0;
2071  a = log10(yy1/y0)/log10(c);
2072  if(a > 10.0) return 0.;
2073  // b0 = log10(y0) - a*log10(x0);
2074  b = y0/pow(x0,a); // pow(10.,b0);
2075 
2076  a += 1.0;
2077  if( a == 0 ) result = b*log(x0/e0);
2078  else result = y0*(x0 - e0*pow(d,a-1))/a;
2079  a += 1.0;
2080 
2081  if( a == 0 ) fIntegralMM[0] += b*log(x0/e0);
2082  else fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2083 
2084 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2085 
2086  x0 = fSplineEnergy[i - 1];
2087  x1 = fSplineEnergy[i - 2];
2088  y0 = fdNdxMM[i - 1];
2089  yy1 = fdNdxMM[i - 2];
2090 
2091  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2092  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2093 
2094  c = x1/x0;
2095  d = e0/x0;
2096  a = log10(yy1/y0)/log10(x1/x0);
2097  // b0 = log10(y0) - a*log10(x0);
2098  b = y0/pow(x0,a); // pow(10.,b0);
2099 
2100  a += 1.0;
2101  if( a == 0 ) result += b*log(e0/x0);
2102  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2103  a += 1.0;
2104 
2105  if( a == 0 ) fIntegralMM[0] += b*log(e0/x0);
2106  else fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2107 
2108  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2109  // <<b<<"; result = "<<result<<G4endl;
2110 
2111  return result;
2112 
2113 }
2114 
2116 //
2117 // Integration of Plasmon cross-section for the case of
2118 // passing across border between intervals
2119 
2121  G4double en0 )
2122 {
2123  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2124 
2125  e0 = en0;
2126  x0 = fSplineEnergy[i];
2127  x1 = fSplineEnergy[i+1];
2128  y0 = fdNdxPlasmon[i];
2129  yy1 = fdNdxPlasmon[i+1];
2130 
2131  c = x1/x0;
2132  d = e0/x0;
2133  a = log10(yy1/y0)/log10(c);
2134  if(a > 10.0) return 0.;
2135  // b0 = log10(y0) - a*log10(x0);
2136  b = y0/pow(x0,a); //pow(10.,b);
2137 
2138  a += 1.0;
2139  if( a == 0 ) result = b*log(x0/e0);
2140  else result = y0*(x0 - e0*pow(d,a-1))/a;
2141  a += 1.0;
2142 
2143  if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
2144  else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2145 
2146  x0 = fSplineEnergy[i - 1];
2147  x1 = fSplineEnergy[i - 2];
2148  y0 = fdNdxPlasmon[i - 1];
2149  yy1 = fdNdxPlasmon[i - 2];
2150 
2151  c = x1/x0;
2152  d = e0/x0;
2153  a = log10(yy1/y0)/log10(c);
2154  // b0 = log10(y0) - a*log10(x0);
2155  b = y0/pow(x0,a);// pow(10.,b0);
2156 
2157  a += 1.0;
2158  if( a == 0 ) result += b*log(e0/x0);
2159  else result += y0*(e0*pow(d,a-1) - x0)/a;
2160  a += 1.0;
2161 
2162  if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
2163  else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2164 
2165  return result;
2166 
2167 }
2168 
2170 //
2171 // Integration of resonance cross-section for the case of
2172 // passing across border between intervals
2173 
2175  G4double en0 )
2176 {
2177  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2178 
2179  e0 = en0;
2180  x0 = fSplineEnergy[i];
2181  x1 = fSplineEnergy[i+1];
2182  y0 = fdNdxResonance[i];
2183  yy1 = fdNdxResonance[i+1];
2184 
2185  c = x1/x0;
2186  d = e0/x0;
2187  a = log10(yy1/y0)/log10(c);
2188  if(a > 10.0) return 0.;
2189  // b0 = log10(y0) - a*log10(x0);
2190  b = y0/pow(x0,a); //pow(10.,b);
2191 
2192  a += 1.0;
2193  if( a == 0 ) result = b*log(x0/e0);
2194  else result = y0*(x0 - e0*pow(d,a-1))/a;
2195  a += 1.0;
2196 
2197  if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0);
2198  else fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2199 
2200  x0 = fSplineEnergy[i - 1];
2201  x1 = fSplineEnergy[i - 2];
2202  y0 = fdNdxResonance[i - 1];
2203  yy1 = fdNdxResonance[i - 2];
2204 
2205  c = x1/x0;
2206  d = e0/x0;
2207  a = log10(yy1/y0)/log10(c);
2208  // b0 = log10(y0) - a*log10(x0);
2209  b = y0/pow(x0,a);// pow(10.,b0);
2210 
2211  a += 1.0;
2212  if( a == 0 ) result += b*log(e0/x0);
2213  else result += y0*(e0*pow(d,a-1) - x0)/a;
2214  a += 1.0;
2215 
2216  if( a == 0 ) fIntegralResonance[0] += b*log(e0/x0);
2217  else fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2218 
2219  return result;
2220 
2221 }
2222 
2224 //
2225 // Returns random PAI-total energy loss over step
2226 
2228 {
2229  G4long numOfCollisions;
2230  G4double meanNumber, loss = 0.0;
2231 
2232  // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl;
2233 
2234  meanNumber = fIntegralPAIxSection[1]*step;
2235  numOfCollisions = G4Poisson(meanNumber);
2236 
2237  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2238 
2239  while(numOfCollisions)
2240  {
2241  loss += GetEnergyTransfer();
2242  numOfCollisions--;
2243  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2244  }
2245  // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
2246 
2247  return loss;
2248 }
2249 
2251 //
2252 // Returns random PAI-total energy transfer in one collision
2253 
2255 {
2256  G4int iTransfer ;
2257 
2258  G4double energyTransfer, position;
2259 
2260  position = fIntegralPAIxSection[1]*G4UniformRand();
2261 
2262  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2263  {
2264  if( position >= fIntegralPAIxSection[iTransfer] ) break;
2265  }
2266  if(iTransfer > fSplineNumber) iTransfer--;
2267 
2268  energyTransfer = fSplineEnergy[iTransfer];
2269 
2270  if(iTransfer > 1)
2271  {
2272  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2273  }
2274  return energyTransfer;
2275 }
2276 
2278 //
2279 // Returns random Cerenkov energy loss over step
2280 
2282 {
2283  G4long numOfCollisions;
2284  G4double meanNumber, loss = 0.0;
2285 
2286  // G4cout<<" G4PAIxSection::GetStepCerenkovLoss "<<G4endl;
2287 
2288  meanNumber = fIntegralCerenkov[1]*step;
2289  numOfCollisions = G4Poisson(meanNumber);
2290 
2291  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2292 
2293  while(numOfCollisions)
2294  {
2295  loss += GetCerenkovEnergyTransfer();
2296  numOfCollisions--;
2297  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2298  }
2299  // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2300 
2301  return loss;
2302 }
2303 
2305 //
2306 // Returns random MM-Cerenkov energy loss over step
2307 
2309 {
2310  G4long numOfCollisions;
2311  G4double meanNumber, loss = 0.0;
2312 
2313  // G4cout<<" G4PAIxSection::GetStepMMLoss "<<G4endl;
2314 
2315  meanNumber = fIntegralMM[1]*step;
2316  numOfCollisions = G4Poisson(meanNumber);
2317 
2318  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2319 
2320  while(numOfCollisions)
2321  {
2322  loss += GetMMEnergyTransfer();
2323  numOfCollisions--;
2324  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2325  }
2326  // G4cout<<"PAI MM-Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2327 
2328  return loss;
2329 }
2330 
2332 //
2333 // Returns Cerenkov energy transfer in one collision
2334 
2336 {
2337  G4int iTransfer ;
2338 
2339  G4double energyTransfer, position;
2340 
2341  position = fIntegralCerenkov[1]*G4UniformRand();
2342 
2343  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2344  {
2345  if( position >= fIntegralCerenkov[iTransfer] ) break;
2346  }
2347  if(iTransfer > fSplineNumber) iTransfer--;
2348 
2349  energyTransfer = fSplineEnergy[iTransfer];
2350 
2351  if(iTransfer > 1)
2352  {
2353  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2354  }
2355  return energyTransfer;
2356 }
2357 
2359 //
2360 // Returns MM-Cerenkov energy transfer in one collision
2361 
2363 {
2364  G4int iTransfer ;
2365 
2366  G4double energyTransfer, position;
2367 
2368  position = fIntegralMM[1]*G4UniformRand();
2369 
2370  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2371  {
2372  if( position >= fIntegralMM[iTransfer] ) break;
2373  }
2374  if(iTransfer > fSplineNumber) iTransfer--;
2375 
2376  energyTransfer = fSplineEnergy[iTransfer];
2377 
2378  if(iTransfer > 1)
2379  {
2380  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2381  }
2382  return energyTransfer;
2383 }
2384 
2386 //
2387 // Returns random plasmon energy loss over step
2388 
2390 {
2391  G4long numOfCollisions;
2392  G4double meanNumber, loss = 0.0;
2393 
2394  // G4cout<<" G4PAIxSection::GetStepPlasmonLoss "<<G4endl;
2395 
2396  meanNumber = fIntegralPlasmon[1]*step;
2397  numOfCollisions = G4Poisson(meanNumber);
2398 
2399  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2400 
2401  while(numOfCollisions)
2402  {
2403  loss += GetPlasmonEnergyTransfer();
2404  numOfCollisions--;
2405  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2406  }
2407  // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
2408 
2409  return loss;
2410 }
2411 
2413 //
2414 // Returns plasmon energy transfer in one collision
2415 
2417 {
2418  G4int iTransfer ;
2419 
2420  G4double energyTransfer, position;
2421 
2422  position = fIntegralPlasmon[1]*G4UniformRand();
2423 
2424  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2425  {
2426  if( position >= fIntegralPlasmon[iTransfer] ) break;
2427  }
2428  if(iTransfer > fSplineNumber) iTransfer--;
2429 
2430  energyTransfer = fSplineEnergy[iTransfer];
2431 
2432  if(iTransfer > 1)
2433  {
2434  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2435  }
2436  return energyTransfer;
2437 }
2438 
2440 //
2441 // Returns random resonance energy loss over step
2442 
2444 {
2445  G4long numOfCollisions;
2446  G4double meanNumber, loss = 0.0;
2447 
2448  // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl;
2449 
2450  meanNumber = fIntegralResonance[1]*step;
2451  numOfCollisions = G4Poisson(meanNumber);
2452 
2453  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2454 
2455  while(numOfCollisions)
2456  {
2457  loss += GetResonanceEnergyTransfer();
2458  numOfCollisions--;
2459  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2460  }
2461  // G4cout<<"PAI resonance loss = "<<loss/keV<<" keV"<<G4endl;
2462 
2463  return loss;
2464 }
2465 
2466 
2468 //
2469 // Returns resonance energy transfer in one collision
2470 
2472 {
2473  G4int iTransfer ;
2474 
2475  G4double energyTransfer, position;
2476 
2477  position = fIntegralResonance[1]*G4UniformRand();
2478 
2479  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2480  {
2481  if( position >= fIntegralResonance[iTransfer] ) break;
2482  }
2483  if(iTransfer > fSplineNumber) iTransfer--;
2484 
2485  energyTransfer = fSplineEnergy[iTransfer];
2486 
2487  if(iTransfer > 1)
2488  {
2489  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2490  }
2491  return energyTransfer;
2492 }
2493 
2494 
2496 //
2497 // Returns Rutherford energy transfer in one collision
2498 
2500 {
2501  G4int iTransfer ;
2502 
2503  G4double energyTransfer, position;
2504 
2505  position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand();
2506 
2507  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2508  {
2509  if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break;
2510  }
2511  if(iTransfer > fSplineNumber) iTransfer--;
2512 
2513  energyTransfer = fSplineEnergy[iTransfer];
2514 
2515  if(iTransfer > 1)
2516  {
2517  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2518  }
2519  return energyTransfer;
2520 }
2521 
2523 //
2524 
2525 void G4PAIxSection::CallError(G4int i, const G4String& methodName) const
2526 {
2527  G4String head = "G4PAIxSection::" + methodName + "()";
2529  ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
2530  G4Exception(head,"pai001",FatalException,ed);
2531 }
2532 
2534 //
2535 // Init array of Lorentz factors
2536 //
2537 
2539 
2540 const G4double G4PAIxSection::fLorentzFactor[112] = // fNumberOfGammas+1
2541 {
2542 0.0,
2543 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
2544 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
2545 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
2546 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
2547 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
2548 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
2549 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
2550 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
2551 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
2552 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
2553 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
2554 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
2555 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
2556 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
2557 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
2558 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
2559 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
2560 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
2561 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
2562 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
2563 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
2564 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
2565 1.065799e+05
2566 };
2567 
2569 //
2570 // The number of gamma for creation of spline (near ion-min , G ~ 4 )
2571 //
2572 
2573 const
2575 
2576 
2577 //
2578 // end of G4PAIxSection implementation file
2579 //
2581 
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:119
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:262
G4double SumOverInterMM(G4int intervalNumber)
static const G4double betaBohr
G4double GetZ() const
Definition: G4Element.hh:131
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:589
std::vector< G4Material * > G4MaterialTable
static const G4double cofBetaBohr
G4double GetDensity() const
Definition: G4Material.hh:180
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 GetSandiaMatTablePAI(G4int, G4int) const
G4double SumOverInterPlasmon(G4int intervalNumber)
G4double a
Definition: TRTMaterials.hh:39
G4double SumOverInterval(G4int intervalNumber)
G4double RePartDielectricConst(G4double energy)
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:202
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 position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double GetElectronDensity() const
Definition: G4Material.hh:217
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:212
G4double GetStepEnergyLoss(G4double step)
static const double pi
Definition: G4SIunits.hh:74
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:180
static const G4double fDelta
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
G4double SumOverInterCerenkov(G4int intervalNumber)
const G4double x[NPOINTSGL]
#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:186
static const double keV
Definition: G4SIunits.hh:213
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)