Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4PAIxSection.cc
이 파일의 문서화 페이지로 가기
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4PAIxSection.cc 108737 2018-03-02 13:49:56Z 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 
92 G4PAIxSection::G4PAIxSection()
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 
129 G4PAIxSection::G4PAIxSection(G4MaterialCutsCouple* matCC)
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 
162 G4PAIxSection::G4PAIxSection(G4int materialIndex,
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 
274 G4PAIxSection::G4PAIxSection( G4int materialIndex,
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 
435 G4PAIxSection::G4PAIxSection( G4int materialIndex,
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 
586 G4PAIxSection::~G4PAIxSection()
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 G4double G4PAIxSection::GetLorentzFactor(G4int j) const
600 {
601  return fLorentzFactor[j];
602 }
603 
605 //
606 // Constructor with beta*gamma square value called from G4PAIPhotModel/Data
607 
608 void G4PAIxSection::Initialize( const G4Material* material,
609  G4double maxEnergyTransfer,
610  G4double betaGammaSq,
611  G4SandiaTable* sandia)
612 {
613  if(fVerbose > 0)
614  {
615  G4cout<<G4endl;
616  G4cout<<"G4PAIxSection::Initialize(...,G4SandiaTable* sandia)"<<G4endl;
617  G4cout<<G4endl;
618  }
619  G4int i, j;
620 
621  fSandia = sandia;
622  fIntervalNumber = sandia->GetMaxInterval();
623  fDensity = material->GetDensity();
624  fElectronDensity = material->GetElectronDensity();
625 
626  // fIntervalNumber--;
627 
628  if( fVerbose > 0 )
629  {
630  G4cout<<"fDensity = "<<fDensity<<"\t"<<fElectronDensity<<"\t fIntervalNumber = "<<fIntervalNumber<<G4endl;
631  }
632  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
633  fA1 = G4DataVector(fIntervalNumber+2,0.0);
634  fA2 = G4DataVector(fIntervalNumber+2,0.0);
635  fA3 = G4DataVector(fIntervalNumber+2,0.0);
636  fA4 = G4DataVector(fIntervalNumber+2,0.0);
637 
638  for( i = 1; i <= fIntervalNumber; i++ )
639  {
640  if ( sandia->GetSandiaMatTablePAI(i-1,0) < 1.*eV && sandia->GetLowerI1() == false)
641  {
642  fIntervalNumber--;
643  continue;
644  }
645  if( ( sandia->GetSandiaMatTablePAI(i-1,0) >= maxEnergyTransfer ) || i >= fIntervalNumber )
646  {
647  fEnergyInterval[i] = maxEnergyTransfer;
648  fIntervalNumber = i;
649  break;
650  }
651  fEnergyInterval[i] = sandia->GetSandiaMatTablePAI(i-1,0);
652  fA1[i] = sandia->GetSandiaMatTablePAI(i-1,1);
653  fA2[i] = sandia->GetSandiaMatTablePAI(i-1,2);
654  fA3[i] = sandia->GetSandiaMatTablePAI(i-1,3);
655  fA4[i] = sandia->GetSandiaMatTablePAI(i-1,4);
656 
657  if( fVerbose > 0 )
658  {
659  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
660  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
661  }
662  }
663  if( fVerbose > 0 ) G4cout<<"last i = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
664 
665  if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer )
666  {
667  fIntervalNumber++;
668  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
669  }
670  if( fVerbose > 0 )
671  {
672  for( i = 1; i <= fIntervalNumber; i++ )
673  {
674  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
675  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
676  }
677  }
678  if( fVerbose > 0 ) G4cout<<"Now checking, if two borders are too close together"<<G4endl;
679 
680  for( i = 1; i < fIntervalNumber; i++ )
681  {
682  if( fEnergyInterval[i+1]-fEnergyInterval[i] >
683  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]) && fEnergyInterval[i] > 0.) continue;
684  else
685  {
686  if( fVerbose > 0 ) G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fEnergyInterval[i+1]/keV;
687 
688  for( j = i; j < fIntervalNumber; j++ )
689  {
690  fEnergyInterval[j] = fEnergyInterval[j+1];
691  fA1[j] = fA1[j+1];
692  fA2[j] = fA2[j+1];
693  fA3[j] = fA3[j+1];
694  fA4[j] = fA4[j+1];
695  }
696  fIntervalNumber--;
697  i--;
698  }
699  }
700  if( fVerbose > 0 )
701  {
702  for( i = 1; i <= fIntervalNumber; i++ )
703  {
704  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
705  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
706  }
707  }
708  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
709 
710  ComputeLowEnergyCof(material);
711 
712  G4double betaGammaSqRef =
713  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
714 
715  NormShift(betaGammaSqRef);
716  SplainPAI(betaGammaSqRef);
717 
718  // Preparation of integral PAI cross section for input betaGammaSq
719 
720  for( i = 1; i <= fSplineNumber; i++ )
721  {
722  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
723 
724 
725  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
726  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
727  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
728  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
729  }
730  IntegralPAIxSection();
731  IntegralCerenkov();
732  IntegralMM();
733  IntegralPlasmon();
734  IntegralResonance();
735 
736  for( i = 1; i <= fSplineNumber; i++ )
737  {
738  if(fVerbose>0) G4cout<<i<<"; w = "<<fSplineEnergy[i]/keV<<" keV; dN/dx_>w = "<<fIntegralPAIxSection[i]<<" 1/mm"<<G4endl;
739  }
740 }
741 
742 
744 //
745 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
746 //
747 
748 void G4PAIxSection::ComputeLowEnergyCof(const G4Material* material)
749 {
750  G4int i, numberOfElements = material->GetNumberOfElements();
751  G4double sumZ = 0., sumCof = 0.;
752 
753  static const G4double p0 = 1.20923e+00;
754  static const G4double p1 = 3.53256e-01;
755  static const G4double p2 = -1.45052e-03;
756 
757  G4double* thisMaterialZ = new G4double[numberOfElements];
758  G4double* thisMaterialCof = new G4double[numberOfElements];
759 
760  for( i = 0; i < numberOfElements; i++ )
761  {
762  thisMaterialZ[i] = material->GetElement(i)->GetZ();
763  sumZ += thisMaterialZ[i];
764  thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
765  }
766  for( i = 0; i < numberOfElements; i++ )
767  {
768  sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
769  }
770  fLowEnergyCof = sumCof;
771  delete [] thisMaterialZ;
772  delete [] thisMaterialCof;
773  // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
774 }
775 
777 //
778 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
779 //
780 
781 void G4PAIxSection::ComputeLowEnergyCof()
782 {
783  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
784  G4int i, numberOfElements = (*theMaterialTable)[fMaterialIndex]->GetNumberOfElements();
785  G4double sumZ = 0., sumCof = 0.;
786 
787  const G4double p0 = 1.20923e+00;
788  const G4double p1 = 3.53256e-01;
789  const G4double p2 = -1.45052e-03;
790 
791  G4double* thisMaterialZ = new G4double[numberOfElements];
792  G4double* thisMaterialCof = new G4double[numberOfElements];
793 
794  for( i = 0; i < numberOfElements; i++ )
795  {
796  thisMaterialZ[i] = (*theMaterialTable)[fMaterialIndex]->GetElement(i)->GetZ();
797  sumZ += thisMaterialZ[i];
798  thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
799  }
800  for( i = 0; i < numberOfElements; i++ )
801  {
802  sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
803  }
804  fLowEnergyCof = sumCof;
805  // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
806  delete [] thisMaterialZ;
807  delete [] thisMaterialCof;
808 }
809 
811 //
812 // General control function for class G4PAIxSection
813 //
814 
815 void G4PAIxSection::InitPAI()
816 {
817  G4int i;
818  G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
819  fLorentzFactor[fRefGammaNumber] - 1;
820 
821  // Preparation of integral PAI cross section for reference gamma
822 
823  NormShift(betaGammaSq);
824  SplainPAI(betaGammaSq);
825 
826  IntegralPAIxSection();
827  IntegralCerenkov();
828  IntegralMM();
829  IntegralPlasmon();
830  IntegralResonance();
831 
832  for(i = 0; i<= fSplineNumber; i++)
833  {
834  fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i];
835  if(i != 0)
836  {
837  fPAItable[i][0] = fSplineEnergy[i];
838  }
839  }
840  fPAItable[0][0] = fSplineNumber;
841 
842  for(G4int j = 1; j < 112; j++) // for other gammas
843  {
844  if( j == fRefGammaNumber ) continue;
845 
846  betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
847 
848  for(i = 1; i <= fSplineNumber; i++)
849  {
850  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
851  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
852  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
853  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
854  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
855  }
856  IntegralPAIxSection();
857  IntegralCerenkov();
858  IntegralMM();
859  IntegralPlasmon();
860  IntegralResonance();
861 
862  for(i = 0; i <= fSplineNumber; i++)
863  {
864  fPAItable[i][j] = fIntegralPAIxSection[i];
865  }
866  }
867 
868 }
869 
871 //
872 // Shifting from borders to intervals Creation of first energy points
873 //
874 
875 void G4PAIxSection::NormShift(G4double betaGammaSq)
876 {
877  G4int i, j;
878 
879  if(fVerbose>0) G4cout<<" G4PAIxSection::NormShift call "<<G4endl;
880 
881 
882  for( i = 1; i <= fIntervalNumber-1; i++ )
883  {
884  for( j = 1; j <= 2; j++ )
885  {
886  fSplineNumber = (i-1)*2 + j;
887 
888  if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
889  else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
890  if(fVerbose>0) G4cout<<"cn = "<<fSplineNumber<<"; "<<"w = "<<fSplineEnergy[fSplineNumber]/keV<<" keV"<<G4endl;
891  }
892  }
893  fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
894 
895  j = 1;
896 
897  for( i = 2; i <= fSplineNumber; i++ )
898  {
899  if( fSplineEnergy[i]<fEnergyInterval[j+1] )
900  {
901  fIntegralTerm[i] = fIntegralTerm[i-1] +
902  RutherfordIntegral(j,fSplineEnergy[i-1],
903  fSplineEnergy[i] );
904  }
905  else
906  {
907  G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
908  fEnergyInterval[j+1] );
909  j++;
910  fIntegralTerm[i] = fIntegralTerm[i-1] + x +
911  RutherfordIntegral(j,fEnergyInterval[j],
912  fSplineEnergy[i] );
913  }
914  if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV \t"<<fIntegralTerm[i]<<"\n"<<G4endl;
915  }
916  fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
917  fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
918 
919  // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl;
920 
921  // Calculation of PAI differrential cross-section (1/(keV*cm))
922  // in the energy points near borders of energy intervals
923 
924  for(G4int k = 1; k <= fIntervalNumber-1; k++ )
925  {
926  for( j = 1; j <= 2; j++ )
927  {
928  i = (k-1)*2 + j;
929  fImPartDielectricConst[i] = fNormalizationCof*
930  ImPartDielectricConst(k,fSplineEnergy[i]);
931  fRePartDielectricConst[i] = fNormalizationCof*
932  RePartDielectricConst(fSplineEnergy[i]);
933  fIntegralTerm[i] *= fNormalizationCof;
934 
935  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
936  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
937  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
938  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
939  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
940  if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV, xsc = "<<fDifPAIxSection[i]<<"\n"<<G4endl;
941  }
942  }
943 
944 } // end of NormShift
945 
947 //
948 // Creation of new energy points as geometrical mean of existing
949 // one, calculation PAI_cs for them, while the error of logarithmic
950 // linear approximation would be smaller than 'fError'
951 
952 void G4PAIxSection::SplainPAI(G4double betaGammaSq)
953 {
954  G4int j, k = 1, i = 1;
955 
956  if(fVerbose>0) G4cout<<" G4PAIxSection::SplainPAI call "<<G4endl;
957 
958  while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
959  {
960  // if( std::abs(fSplineEnergy[i+1]-fEnergyInterval[k+1]) > (fSplineEnergy[i+1]+fEnergyInterval[k+1])*5.e-7 )
961  if( fSplineEnergy[i+1] > fEnergyInterval[k+1] )
962  {
963  k++; // Here next energy point is in next energy interval
964  i++;
965  if(fVerbose>0) G4cout<<" in if: i = "<<i<<"; k = "<<k<<G4endl;
966  continue;
967  }
968  if(fVerbose>0) G4cout<<" out if: i = "<<i<<"; k = "<<k<<G4endl;
969 
970  // Shifting of arrayes for inserting the geometrical
971  // average of 'i' and 'i+1' energy points to 'i+1' place
972  fSplineNumber++;
973 
974  for( j = fSplineNumber; j >= i+2; j-- )
975  {
976  fSplineEnergy[j] = fSplineEnergy[j-1];
977  fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
978  fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
979  fIntegralTerm[j] = fIntegralTerm[j-1];
980 
981  fDifPAIxSection[j] = fDifPAIxSection[j-1];
982  fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
983  fdNdxMM[j] = fdNdxMM[j-1];
984  fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
985  fdNdxResonance[j] = fdNdxResonance[j-1];
986  }
987  G4double x1 = fSplineEnergy[i];
988  G4double x2 = fSplineEnergy[i+1];
989  G4double yy1 = fDifPAIxSection[i];
990  G4double y2 = fDifPAIxSection[i+1];
991 
992  if(fVerbose>0) G4cout<<"Spline: x1 = "<<x1<<"; x2 = "<<x2<<", yy1 = "<<yy1<<"; y2 = "<<y2<<G4endl;
993 
994 
995  G4double en1 = sqrt(x1*x2);
996  // G4double en1 = 0.5*(x1 + x2);
997 
998 
999  fSplineEnergy[i+1] = en1;
1000 
1001  // Calculation of logarithmic linear approximation
1002  // in this (enr) energy point, which number is 'i+1' now
1003 
1004  G4double a = log10(y2/yy1)/log10(x2/x1);
1005  G4double b = log10(yy1) - a*log10(x1);
1006  G4double y = a*log10(en1) + b;
1007 
1008  y = pow(10.,y);
1009 
1010  // Calculation of the PAI dif. cross-section at this point
1011 
1012  fImPartDielectricConst[i+1] = fNormalizationCof*
1013  ImPartDielectricConst(k,fSplineEnergy[i+1]);
1014  fRePartDielectricConst[i+1] = fNormalizationCof*
1015  RePartDielectricConst(fSplineEnergy[i+1]);
1016  fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
1017  RutherfordIntegral(k,fSplineEnergy[i],
1018  fSplineEnergy[i+1]);
1019 
1020  fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
1021  fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
1022  fdNdxMM[i+1] = PAIdNdxMM(i+1,betaGammaSq);
1023  fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
1024  fdNdxResonance[i+1] = PAIdNdxResonance(i+1,betaGammaSq);
1025 
1026  // Condition for next division of this segment or to pass
1027 
1028  if(fVerbose>0) G4cout<<"Spline, a = "<<a<<"; b = "<<b<<"; new xsc = "<<y<<"; compxsc = "<<fDifPAIxSection[i+1]<<G4endl;
1029 
1030  // to higher energies
1031 
1032  G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
1033 
1034  G4double delta = 2.*(fSplineEnergy[i+1]-fSplineEnergy[i])/(fSplineEnergy[i+1]+fSplineEnergy[i]);
1035 
1036  if( x < 0 )
1037  {
1038  x = -x;
1039  }
1040  if( x > fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*fDelta )
1041  {
1042  continue; // next division
1043  }
1044  i += 2; // pass to next segment
1045 
1046  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1047  } // close 'while'
1048 
1049 } // end of SplainPAI
1050 
1051 
1053 //
1054 // Integration over electrons that could be considered
1055 // quasi-free at energy transfer of interest
1056 
1057 G4double G4PAIxSection::RutherfordIntegral( G4int k,
1058  G4double x1,
1059  G4double x2 )
1060 {
1061  G4double c1, c2, c3;
1062  // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
1063  c1 = (x2 - x1)/x1/x2;
1064  c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2;
1065  c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
1066  // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
1067 
1068  return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
1069 
1070 } // end of RutherfordIntegral
1071 
1072 
1074 //
1075 // Imaginary part of dielectric constant
1076 // (G4int k - interval number, G4double en1 - energy point)
1077 
1078 G4double G4PAIxSection::ImPartDielectricConst( G4int k ,
1079  G4double energy1 )
1080 {
1081  G4double energy2,energy3,energy4,result;
1082 
1083  energy2 = energy1*energy1;
1084  energy3 = energy2*energy1;
1085  energy4 = energy3*energy1;
1086 
1087  result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
1088  result *=hbarc/energy1;
1089 
1090  return result;
1091 
1092 } // end of ImPartDielectricConst
1093 
1095 //
1096 // Returns lambda of photon with energy1 in current material
1097 
1098 G4double G4PAIxSection::GetPhotonRange( G4double energy1 )
1099 {
1100  G4int i;
1101  G4double energy2, energy3, energy4, result, lambda;
1102 
1103  energy2 = energy1*energy1;
1104  energy3 = energy2*energy1;
1105  energy4 = energy3*energy1;
1106 
1107  // G4double* SandiaCof = fSandia->GetSandiaCofForMaterialPAI(energy1);
1108  // result = SandiaCof[0]/energy1+SandiaCof[1]/energy2+SandiaCof[2]/energy3+SandiaCof[3]/energy4;
1109  // result *= fDensity;
1110 
1111  for( i = 1; i <= fIntervalNumber; i++ )
1112  {
1113  if( energy1 < fEnergyInterval[i]) break;
1114  }
1115  i--;
1116  if(i == 0) i = 1;
1117 
1118  result = fA1[i]/energy1+fA2[i]/energy2+fA3[i]/energy3+fA4[i]/energy4;
1119 
1120  if( result > DBL_MIN ) lambda = 1./result;
1121  else lambda = DBL_MAX;
1122 
1123  return lambda;
1124 }
1125 
1127 //
1128 // Return lambda of electron with energy1 in current material
1129 // parametrisation from NIM A554(2005)474-493
1130 
1131 G4double G4PAIxSection::GetElectronRange( G4double energy )
1132 {
1133  G4double range;
1134  /*
1135  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1136 
1137  G4double Z = (*theMaterialTable)[fMaterialIndex]->GetIonisation()->GetZeffective();
1138  G4double A = (*theMaterialTable)[fMaterialIndex]->GetA();
1139 
1140  energy /= keV; // energy in keV in parametrised formula
1141 
1142  if (energy < 10.)
1143  {
1144  range = 3.872e-3*A/Z;
1145  range *= pow(energy,1.492);
1146  }
1147  else
1148  {
1149  range = 6.97e-3*pow(energy,1.6);
1150  }
1151  */
1152  // Blum&Rolandi Particle Detection with Drift Chambers, p. 7
1153 
1154  G4double cofA = 5.37e-4*g/cm2/keV;
1155  G4double cofB = 0.9815;
1156  G4double cofC = 3.123e-3/keV;
1157  // energy /= keV;
1158 
1159  range = cofA*energy*( 1 - cofB/(1 + cofC*energy) );
1160 
1161  // range *= g/cm2;
1162  range /= fDensity;
1163 
1164  return range;
1165 }
1166 
1168 //
1169 // Real part of dielectric constant minus unit: epsilon_1 - 1
1170 // (G4double enb - energy point)
1171 //
1172 
1173 G4double G4PAIxSection::RePartDielectricConst(G4double enb)
1174 {
1175  G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
1176  c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
1177 
1178  x0 = enb;
1179  result = 0;
1180 
1181  for(G4int i=1;i<=fIntervalNumber-1;i++)
1182  {
1183  x1 = fEnergyInterval[i];
1184  x2 = fEnergyInterval[i+1];
1185  xx1 = x1 - x0;
1186  xx2 = x2 - x0;
1187  xx12 = xx2/xx1;
1188 
1189  if(xx12<0)
1190  {
1191  xx12 = -xx12;
1192  }
1193  xln1 = log(x2/x1);
1194  xln2 = log(xx12);
1195  xln3 = log((x2 + x0)/(x1 + x0));
1196  x02 = x0*x0;
1197  x03 = x02*x0;
1198  x04 = x03*x0;
1199  x05 = x04*x0;
1200  c1 = (x2 - x1)/x1/x2;
1201  c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2;
1202  c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
1203 
1204  result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
1205  result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
1206  result -= fA3[i]*c2/2/x02;
1207  result -= fA4[i]*c3/3/x02;
1208 
1209  cof1 = fA1[i]/x02 + fA3[i]/x04;
1210  cof2 = fA2[i]/x03 + fA4[i]/x05;
1211 
1212  result += 0.5*(cof1 +cof2)*xln2;
1213  result += 0.5*(cof1 - cof2)*xln3;
1214  }
1215  result *= 2*hbarc/pi;
1216 
1217  return result;
1218 
1219 } // end of RePartDielectricConst
1220 
1222 //
1223 // PAI differential cross-section in terms of
1224 // simplified Allison's equation
1225 //
1226 
1227 G4double G4PAIxSection::DifPAIxSection( G4int i ,
1228  G4double betaGammaSq )
1229 {
1230  G4double cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
1231 
1232  G4double betaBohr = fine_structure_const;
1233  // G4double betaBohr2 = fine_structure_const*fine_structure_const;
1234  // G4double betaBohr3 = betaBohr*betaBohr2; // *4.0;
1235 
1236  G4double be2 = betaGammaSq/(1 + betaGammaSq);
1237  G4double beta = sqrt(be2);
1238  // G4double be3 = beta*be2;
1239 
1240  cof = 1.;
1241  x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
1242 
1243  if( betaGammaSq < 0.01 ) x2 = log(be2);
1244  else
1245  {
1246  x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1247  (1/betaGammaSq - fRePartDielectricConst[i]) +
1248  fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
1249  }
1250  if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
1251  {
1252  x6 = 0.;
1253  }
1254  else
1255  {
1256  x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
1257  x5 = -1 - fRePartDielectricConst[i] +
1258  be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1259  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1260 
1261  x7 = atan2(fImPartDielectricConst[i],x3);
1262  x6 = x5 * x7;
1263  }
1264  // if(fImPartDielectricConst[i] == 0) x6 = 0.;
1265 
1266  x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
1267 
1268  // if( x4 < 0.0 ) x4 = 0.0;
1269 
1270  x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1271  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1272 
1273  result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
1274 
1275  if( result < 1.0e-8 ) result = 1.0e-8;
1276 
1277  result *= fine_structure_const/be2/pi;
1278 
1279  // low energy correction
1280 
1281  G4double lowCof = fLowEnergyCof; // 6.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)?
1282 
1283  result *= (1 - exp(-beta/betaBohr/lowCof));
1284 
1285 
1286  // result *= (1 - exp(-be2/betaBohr2/lowCof));
1287 
1288  // result *= (1 - exp(-be3/betaBohr3/lowCof)); // ~ be for be<<betaBohr
1289 
1290  // result *= (1 - exp(-be4/betaBohr4/lowCof));
1291 
1292  if(fDensity >= 0.1)
1293  {
1294  result /= x8;
1295  }
1296  return result;
1297 
1298 } // end of DifPAIxSection
1299 
1301 //
1302 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
1303 
1304 G4double G4PAIxSection::PAIdNdxCerenkov( G4int i ,
1305  G4double betaGammaSq )
1306 {
1307  G4double logarithm, x3, x5, argument, modul2, dNdxC;
1308  G4double be2, betaBohr2, cofBetaBohr;
1309 
1310  cofBetaBohr = 4.0;
1311  betaBohr2 = fine_structure_const*fine_structure_const;
1312  G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1313 
1314  be2 = betaGammaSq/(1 + betaGammaSq);
1315  G4double be4 = be2*be2;
1316 
1317  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1318  else
1319  {
1320  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1321  (1/betaGammaSq - fRePartDielectricConst[i]) +
1322  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1323  logarithm += log(1+1.0/betaGammaSq);
1324  }
1325 
1326  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1327  {
1328  argument = 0.0;
1329  }
1330  else
1331  {
1332  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1333  x5 = -1.0 - fRePartDielectricConst[i] +
1334  be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1335  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1336  if( x3 == 0.0 ) argument = 0.5*pi;
1337  else argument = atan2(fImPartDielectricConst[i],x3);
1338  argument *= x5 ;
1339  }
1340  dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
1341 
1342  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1343 
1344  dNdxC *= fine_structure_const/be2/pi;
1345 
1346  dNdxC *= (1-exp(-be4/betaBohr4));
1347 
1348  if(fDensity >= 0.1)
1349  {
1350  modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1351  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1352  dNdxC /= modul2;
1353  }
1354  return dNdxC;
1355 
1356 } // end of PAIdNdxCerenkov
1357 
1359 //
1360 // Calculation od dN/dx of collisions of MM with creation of Cerenkov pseudo-photons
1361 
1362 G4double G4PAIxSection::PAIdNdxMM( G4int i ,
1363  G4double betaGammaSq )
1364 {
1365  G4double logarithm, x3, x5, argument, dNdxC;
1366  G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
1367 
1368  cofBetaBohr = 4.0;
1369  betaBohr2 = fine_structure_const*fine_structure_const;
1370  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1371 
1372  be2 = betaGammaSq/(1 + betaGammaSq);
1373  be4 = be2*be2;
1374 
1375  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1376  else
1377  {
1378  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1379  (1/betaGammaSq - fRePartDielectricConst[i]) +
1380  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1381  logarithm += log(1+1.0/betaGammaSq);
1382  }
1383 
1384  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1385  {
1386  argument = 0.0;
1387  }
1388  else
1389  {
1390  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1391  x5 = be2*( 1.0 + fRePartDielectricConst[i] ) - 1.0;
1392  if( x3 == 0.0 ) argument = 0.5*pi;
1393  else argument = atan2(fImPartDielectricConst[i],x3);
1394  argument *= x5 ;
1395  }
1396  dNdxC = ( logarithm*fImPartDielectricConst[i]*be2 + argument )/hbarc;
1397 
1398  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1399 
1400  dNdxC *= fine_structure_const/be2/pi;
1401 
1402  dNdxC *= (1-exp(-be4/betaBohr4));
1403  return dNdxC;
1404 
1405 } // end of PAIdNdxMM
1406 
1408 //
1409 // Calculation od dN/dx of collisions with creation of longitudinal EM
1410 // excitations (plasmons, delta-electrons)
1411 
1412 G4double G4PAIxSection::PAIdNdxPlasmon( G4int i ,
1413  G4double betaGammaSq )
1414 {
1415  G4double resonance, modul2, dNdxP, cof = 1.;
1416  G4double be2, betaBohr;
1417 
1418  betaBohr = fine_structure_const;
1419  be2 = betaGammaSq/(1 + betaGammaSq);
1420 
1421  G4double beta = sqrt(be2);
1422 
1423  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1424  resonance *= fImPartDielectricConst[i]/hbarc;
1425 
1426 
1427  dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
1428 
1429  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1430 
1431  dNdxP *= fine_structure_const/be2/pi;
1432 
1433  dNdxP *= (1 - exp(-beta/betaBohr/fLowEnergyCof));
1434 
1435  // dNdxP *= (1-exp(-be4/betaBohr4));
1436 
1437  if( fDensity >= 0.1 )
1438  {
1439  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1440  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1441  dNdxP /= modul2;
1442  }
1443  return dNdxP;
1444 
1445 } // end of PAIdNdxPlasmon
1446 
1448 //
1449 // Calculation od dN/dx of collisions with creation of longitudinal EM
1450 // resonance excitations (plasmons, delta-electrons)
1451 
1452 G4double G4PAIxSection::PAIdNdxResonance( G4int i ,
1453  G4double betaGammaSq )
1454 {
1455  G4double resonance, modul2, dNdxP;
1456  G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
1457 
1458  cofBetaBohr = 4.0;
1459  betaBohr2 = fine_structure_const*fine_structure_const;
1460  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1461 
1462  be2 = betaGammaSq/(1 + betaGammaSq);
1463  be4 = be2*be2;
1464 
1465  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1466  resonance *= fImPartDielectricConst[i]/hbarc;
1467 
1468 
1469  dNdxP = resonance;
1470 
1471  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1472 
1473  dNdxP *= fine_structure_const/be2/pi;
1474  dNdxP *= (1-exp(-be4/betaBohr4));
1475 
1476  if( fDensity >= 0.1 )
1477  {
1478  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1479  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1480  dNdxP /= modul2;
1481  }
1482  return dNdxP;
1483 
1484 } // end of PAIdNdxResonance
1485 
1487 //
1488 // Calculation of the PAI integral cross-section
1489 // fIntegralPAIxSection[1] = specific primary ionisation, 1/cm
1490 // and fIntegralPAIxSection[0] = mean energy loss per cm in keV/cm
1491 
1492 void G4PAIxSection::IntegralPAIxSection()
1493 {
1494  fIntegralPAIxSection[fSplineNumber] = 0;
1495  fIntegralPAIdEdx[fSplineNumber] = 0;
1496  fIntegralPAIxSection[0] = 0;
1497  G4int i, k = fIntervalNumber -1;
1498 
1499  for( i = fSplineNumber-1; i >= 1; i--)
1500  {
1501  if(fSplineEnergy[i] >= fEnergyInterval[k])
1502  {
1503  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i);
1504  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
1505  }
1506  else
1507  {
1508  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] +
1509  SumOverBorder(i+1,fEnergyInterval[k]);
1510  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
1511  SumOverBorderdEdx(i+1,fEnergyInterval[k]);
1512  k--;
1513  }
1514  if(fVerbose>0) G4cout<<"i = "<<i<<"; k = "<<k<<"; intPAIxsc[i] = "<<fIntegralPAIxSection[i]<<G4endl;
1515  }
1516 } // end of IntegralPAIxSection
1517 
1519 //
1520 // Calculation of the PAI Cerenkov integral cross-section
1521 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
1522 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
1523 
1524 void G4PAIxSection::IntegralCerenkov()
1525 {
1526  G4int i, k;
1527  fIntegralCerenkov[fSplineNumber] = 0;
1528  fIntegralCerenkov[0] = 0;
1529  k = fIntervalNumber -1;
1530 
1531  for( i = fSplineNumber-1; i >= 1; i-- )
1532  {
1533  if(fSplineEnergy[i] >= fEnergyInterval[k])
1534  {
1535  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
1536  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1537  }
1538  else
1539  {
1540  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
1541  SumOverBordCerenkov(i+1,fEnergyInterval[k]);
1542  k--;
1543  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1544  }
1545  }
1546 
1547 } // end of IntegralCerenkov
1548 
1550 //
1551 // Calculation of the PAI MM-Cerenkov integral cross-section
1552 // fIntegralMM[1] = specific MM-Cerenkov ionisation, 1/cm
1553 // and fIntegralMM[0] = mean MM-Cerenkov loss per cm in keV/cm
1554 
1555 void G4PAIxSection::IntegralMM()
1556 {
1557  G4int i, k;
1558  fIntegralMM[fSplineNumber] = 0;
1559  fIntegralMM[0] = 0;
1560  k = fIntervalNumber -1;
1561 
1562  for( i = fSplineNumber-1; i >= 1; i-- )
1563  {
1564  if(fSplineEnergy[i] >= fEnergyInterval[k])
1565  {
1566  fIntegralMM[i] = fIntegralMM[i+1] + SumOverInterMM(i);
1567  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1568  }
1569  else
1570  {
1571  fIntegralMM[i] = fIntegralMM[i+1] +
1572  SumOverBordMM(i+1,fEnergyInterval[k]);
1573  k--;
1574  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1575  }
1576  }
1577 
1578 } // end of IntegralMM
1579 
1581 //
1582 // Calculation of the PAI Plasmon integral cross-section
1583 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
1584 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
1585 
1586 void G4PAIxSection::IntegralPlasmon()
1587 {
1588  fIntegralPlasmon[fSplineNumber] = 0;
1589  fIntegralPlasmon[0] = 0;
1590  G4int k = fIntervalNumber -1;
1591  for(G4int i=fSplineNumber-1;i>=1;i--)
1592  {
1593  if(fSplineEnergy[i] >= fEnergyInterval[k])
1594  {
1595  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
1596  }
1597  else
1598  {
1599  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
1600  SumOverBordPlasmon(i+1,fEnergyInterval[k]);
1601  k--;
1602  }
1603  }
1604 
1605 } // end of IntegralPlasmon
1606 
1608 //
1609 // Calculation of the PAI resonance integral cross-section
1610 // fIntegralResonance[1] = resonance primary ionisation, 1/cm
1611 // and fIntegralResonance[0] = mean resonance loss per cm in keV/cm
1612 
1613 void G4PAIxSection::IntegralResonance()
1614 {
1615  fIntegralResonance[fSplineNumber] = 0;
1616  fIntegralResonance[0] = 0;
1617  G4int k = fIntervalNumber -1;
1618  for(G4int i=fSplineNumber-1;i>=1;i--)
1619  {
1620  if(fSplineEnergy[i] >= fEnergyInterval[k])
1621  {
1622  fIntegralResonance[i] = fIntegralResonance[i+1] + SumOverInterResonance(i);
1623  }
1624  else
1625  {
1626  fIntegralResonance[i] = fIntegralResonance[i+1] +
1627  SumOverBordResonance(i+1,fEnergyInterval[k]);
1628  k--;
1629  }
1630  }
1631 
1632 } // end of IntegralResonance
1633 
1635 //
1636 // Calculation the PAI integral cross-section inside
1637 // of interval of continuous values of photo-ionisation
1638 // cross-section. Parameter 'i' is the number of interval.
1639 
1640 G4double G4PAIxSection::SumOverInterval( G4int i )
1641 {
1642  G4double x0,x1,y0,yy1,a,b,c,result;
1643 
1644  x0 = fSplineEnergy[i];
1645  x1 = fSplineEnergy[i+1];
1646  if(fVerbose>0) G4cout<<"SumOverInterval i= " << i << " x0 = "<<x0<<"; x1 = "<<x1<<G4endl;
1647 
1648  if( x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1649 
1650  y0 = fDifPAIxSection[i];
1651  yy1 = fDifPAIxSection[i+1];
1652 
1653  if(fVerbose>0) G4cout<<"x0 = "<<x0<<"; x1 = "<<x1<<", y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1654 
1655  c = x1/x0;
1656  a = log10(yy1/y0)/log10(c);
1657 
1658  if(fVerbose>0) G4cout<<"SumOverInterval, a = "<<a<<"; c = "<<c<<G4endl;
1659 
1660  // b = log10(y0) - a*log10(x0);
1661  b = y0/pow(x0,a);
1662  a += 1.;
1663  if( std::abs(a) < 1.e-6 )
1664  {
1665  result = b*log(x1/x0);
1666  }
1667  else
1668  {
1669  result = y0*(x1*pow(c,a-1) - x0)/a;
1670  }
1671  a += 1.;
1672  if( std::abs(a) < 1.e-6 )
1673  {
1674  fIntegralPAIxSection[0] += b*log(x1/x0);
1675  }
1676  else
1677  {
1678  fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1679  }
1680  if(fVerbose>0) G4cout<<"SumOverInterval, result = "<<result<<G4endl;
1681  return result;
1682 
1683 } // end of SumOverInterval
1684 
1686 
1687 G4double G4PAIxSection::SumOverIntervaldEdx( G4int i )
1688 {
1689  G4double x0,x1,y0,yy1,a,b,c,result;
1690 
1691  x0 = fSplineEnergy[i];
1692  x1 = fSplineEnergy[i+1];
1693 
1694  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1695 
1696  y0 = fDifPAIxSection[i];
1697  yy1 = fDifPAIxSection[i+1];
1698  c = x1/x0;
1699  a = log10(yy1/y0)/log10(c);
1700  // b = log10(y0) - a*log10(x0);
1701  b = y0/pow(x0,a);
1702  a += 2;
1703  if(a == 0)
1704  {
1705  result = b*log(x1/x0);
1706  }
1707  else
1708  {
1709  result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1710  }
1711  return result;
1712 
1713 } // end of SumOverInterval
1714 
1716 //
1717 // Calculation the PAI Cerenkov integral cross-section inside
1718 // of interval of continuous values of photo-ionisation Cerenkov
1719 // cross-section. Parameter 'i' is the number of interval.
1720 
1721 G4double G4PAIxSection::SumOverInterCerenkov( G4int i )
1722 {
1723  G4double x0,x1,y0,yy1,a,b,c,result;
1724 
1725  x0 = fSplineEnergy[i];
1726  x1 = fSplineEnergy[i+1];
1727 
1728  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1729 
1730  y0 = fdNdxCerenkov[i];
1731  yy1 = fdNdxCerenkov[i+1];
1732  // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1733  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1734 
1735  c = x1/x0;
1736  a = log10(yy1/y0)/log10(c);
1737  b = y0/pow(x0,a);
1738 
1739  a += 1.0;
1740  if(a == 0) result = b*log(c);
1741  else result = y0*(x1*pow(c,a-1) - x0)/a;
1742  a += 1.0;
1743 
1744  if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
1745  else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1746  // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1747  return result;
1748 
1749 } // end of SumOverInterCerenkov
1750 
1752 //
1753 // Calculation the PAI MM-Cerenkov integral cross-section inside
1754 // of interval of continuous values of photo-ionisation Cerenkov
1755 // cross-section. Parameter 'i' is the number of interval.
1756 
1757 G4double G4PAIxSection::SumOverInterMM( G4int i )
1758 {
1759  G4double x0,x1,y0,yy1,a,b,c,result;
1760 
1761  x0 = fSplineEnergy[i];
1762  x1 = fSplineEnergy[i+1];
1763 
1764  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1765 
1766  y0 = fdNdxMM[i];
1767  yy1 = fdNdxMM[i+1];
1768  //G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1769  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1770 
1771  c = x1/x0;
1772  //G4cout<<" c = "<<c<< " yy1/y0= " << yy1/y0 <<G4endl;
1773  a = log10(yy1/y0)/log10(c);
1774  if(a > 10.0) return 0.;
1775  b = y0/pow(x0,a);
1776 
1777  a += 1.0;
1778  if(a == 0) result = b*log(c);
1779  else result = y0*(x1*pow(c,a-1) - x0)/a;
1780  a += 1.0;
1781 
1782  if( a == 0 ) fIntegralMM[0] += b*log(c);
1783  else fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1784  //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1785  return result;
1786 
1787 } // end of SumOverInterMM
1788 
1790 //
1791 // Calculation the PAI Plasmon integral cross-section inside
1792 // of interval of continuous values of photo-ionisation Plasmon
1793 // cross-section. Parameter 'i' is the number of interval.
1794 
1795 G4double G4PAIxSection::SumOverInterPlasmon( G4int i )
1796 {
1797  G4double x0,x1,y0,yy1,a,b,c,result;
1798 
1799  x0 = fSplineEnergy[i];
1800  x1 = fSplineEnergy[i+1];
1801 
1802  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1803 
1804  y0 = fdNdxPlasmon[i];
1805  yy1 = fdNdxPlasmon[i+1];
1806  c =x1/x0;
1807  a = log10(yy1/y0)/log10(c);
1808  if(a > 10.0) return 0.;
1809  // b = log10(y0) - a*log10(x0);
1810  b = y0/pow(x0,a);
1811 
1812  a += 1.0;
1813  if(a == 0) result = b*log(x1/x0);
1814  else result = y0*(x1*pow(c,a-1) - x0)/a;
1815  a += 1.0;
1816 
1817  if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
1818  else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1819 
1820  return result;
1821 
1822 } // end of SumOverInterPlasmon
1823 
1825 //
1826 // Calculation the PAI resonance integral cross-section inside
1827 // of interval of continuous values of photo-ionisation resonance
1828 // cross-section. Parameter 'i' is the number of interval.
1829 
1830 G4double G4PAIxSection::SumOverInterResonance( G4int i )
1831 {
1832  G4double x0,x1,y0,yy1,a,b,c,result;
1833 
1834  x0 = fSplineEnergy[i];
1835  x1 = fSplineEnergy[i+1];
1836 
1837  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1838 
1839  y0 = fdNdxResonance[i];
1840  yy1 = fdNdxResonance[i+1];
1841  c =x1/x0;
1842  a = log10(yy1/y0)/log10(c);
1843  if(a > 10.0) return 0.;
1844  // b = log10(y0) - a*log10(x0);
1845  b = y0/pow(x0,a);
1846 
1847  a += 1.0;
1848  if(a == 0) result = b*log(x1/x0);
1849  else result = y0*(x1*pow(c,a-1) - x0)/a;
1850  a += 1.0;
1851 
1852  if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0);
1853  else fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1854 
1855  return result;
1856 
1857 } // end of SumOverInterResonance
1858 
1860 //
1861 // Integration of PAI cross-section for the case of
1862 // passing across border between intervals
1863 
1864 G4double G4PAIxSection::SumOverBorder( G4int i ,
1865  G4double en0 )
1866 {
1867  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1868 
1869  e0 = en0;
1870  x0 = fSplineEnergy[i];
1871  x1 = fSplineEnergy[i+1];
1872  y0 = fDifPAIxSection[i];
1873  yy1 = fDifPAIxSection[i+1];
1874 
1875  //c = x1/x0;
1876  d = e0/x0;
1877  a = log10(yy1/y0)/log10(x1/x0);
1878  if(a > 10.0) return 0.;
1879 
1880  if(fVerbose>0) G4cout<<"SumOverBorder, a = "<<a<<G4endl;
1881 
1882  // b0 = log10(y0) - a*log10(x0);
1883  b = y0/pow(x0,a); // pow(10.,b);
1884 
1885  a += 1.;
1886  if( std::abs(a) < 1.e-6 )
1887  {
1888  result = b*log(x0/e0);
1889  }
1890  else
1891  {
1892  result = y0*(x0 - e0*pow(d,a-1))/a;
1893  }
1894  a += 1.;
1895  if( std::abs(a) < 1.e-6 )
1896  {
1897  fIntegralPAIxSection[0] += b*log(x0/e0);
1898  }
1899  else
1900  {
1901  fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1902  }
1903  x0 = fSplineEnergy[i - 1];
1904  x1 = fSplineEnergy[i - 2];
1905  y0 = fDifPAIxSection[i - 1];
1906  yy1 = fDifPAIxSection[i - 2];
1907 
1908  //c = x1/x0;
1909  d = e0/x0;
1910  a = log10(yy1/y0)/log10(x1/x0);
1911  // b0 = log10(y0) - a*log10(x0);
1912  b = y0/pow(x0,a);
1913  a += 1.;
1914  if( std::abs(a) < 1.e-6 )
1915  {
1916  result += b*log(e0/x0);
1917  }
1918  else
1919  {
1920  result += y0*(e0*pow(d,a-1) - x0)/a;
1921  }
1922  a += 1.;
1923  if( std::abs(a) < 1.e-6 )
1924  {
1925  fIntegralPAIxSection[0] += b*log(e0/x0);
1926  }
1927  else
1928  {
1929  fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1930  }
1931  return result;
1932 
1933 }
1934 
1936 
1937 G4double G4PAIxSection::SumOverBorderdEdx( G4int i ,
1938  G4double en0 )
1939 {
1940  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1941 
1942  e0 = en0;
1943  x0 = fSplineEnergy[i];
1944  x1 = fSplineEnergy[i+1];
1945  y0 = fDifPAIxSection[i];
1946  yy1 = fDifPAIxSection[i+1];
1947 
1948  //c = x1/x0;
1949  d = e0/x0;
1950  a = log10(yy1/y0)/log10(x1/x0);
1951  if(a > 10.0) return 0.;
1952  // b0 = log10(y0) - a*log10(x0);
1953  b = y0/pow(x0,a); // pow(10.,b);
1954 
1955  a += 2;
1956  if(a == 0)
1957  {
1958  result = b*log(x0/e0);
1959  }
1960  else
1961  {
1962  result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1963  }
1964  x0 = fSplineEnergy[i - 1];
1965  x1 = fSplineEnergy[i - 2];
1966  y0 = fDifPAIxSection[i - 1];
1967  yy1 = fDifPAIxSection[i - 2];
1968 
1969  // c = x1/x0;
1970  d = e0/x0;
1971  a = log10(yy1/y0)/log10(x1/x0);
1972  // b0 = log10(y0) - a*log10(x0);
1973  b = y0/pow(x0,a);
1974  a += 2;
1975  if(a == 0)
1976  {
1977  result += b*log(e0/x0);
1978  }
1979  else
1980  {
1981  result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1982  }
1983  return result;
1984 
1985 }
1986 
1988 //
1989 // Integration of Cerenkov cross-section for the case of
1990 // passing across border between intervals
1991 
1992 G4double G4PAIxSection::SumOverBordCerenkov( G4int i ,
1993  G4double en0 )
1994 {
1995  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1996 
1997  e0 = en0;
1998  x0 = fSplineEnergy[i];
1999  x1 = fSplineEnergy[i+1];
2000  y0 = fdNdxCerenkov[i];
2001  yy1 = fdNdxCerenkov[i+1];
2002 
2003  // G4cout<<G4endl;
2004  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2005  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2006  c = x1/x0;
2007  d = e0/x0;
2008  a = log10(yy1/y0)/log10(c);
2009  if(a > 10.0) return 0.;
2010  // b0 = log10(y0) - a*log10(x0);
2011  b = y0/pow(x0,a); // pow(10.,b0);
2012 
2013  a += 1.0;
2014  if( a == 0 ) result = b*log(x0/e0);
2015  else result = y0*(x0 - e0*pow(d,a-1))/a;
2016  a += 1.0;
2017 
2018  if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
2019  else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2020 
2021 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2022 
2023  x0 = fSplineEnergy[i - 1];
2024  x1 = fSplineEnergy[i - 2];
2025  y0 = fdNdxCerenkov[i - 1];
2026  yy1 = fdNdxCerenkov[i - 2];
2027 
2028  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2029  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2030 
2031  c = x1/x0;
2032  d = e0/x0;
2033  a = log10(yy1/y0)/log10(x1/x0);
2034  // b0 = log10(y0) - a*log10(x0);
2035  b = y0/pow(x0,a); // pow(10.,b0);
2036 
2037  a += 1.0;
2038  if( a == 0 ) result += b*log(e0/x0);
2039  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2040  a += 1.0;
2041 
2042  if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
2043  else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2044 
2045  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2046  // <<b<<"; result = "<<result<<G4endl;
2047 
2048  return result;
2049 
2050 }
2051 
2053 //
2054 // Integration of MM-Cerenkov cross-section for the case of
2055 // passing across border between intervals
2056 
2057 G4double G4PAIxSection::SumOverBordMM( G4int i ,
2058  G4double en0 )
2059 {
2060  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
2061 
2062  e0 = en0;
2063  x0 = fSplineEnergy[i];
2064  x1 = fSplineEnergy[i+1];
2065  y0 = fdNdxMM[i];
2066  yy1 = fdNdxMM[i+1];
2067 
2068  // G4cout<<G4endl;
2069  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2070  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2071  c = x1/x0;
2072  d = e0/x0;
2073  a = log10(yy1/y0)/log10(c);
2074  if(a > 10.0) return 0.;
2075  // b0 = log10(y0) - a*log10(x0);
2076  b = y0/pow(x0,a); // pow(10.,b0);
2077 
2078  a += 1.0;
2079  if( a == 0 ) result = b*log(x0/e0);
2080  else result = y0*(x0 - e0*pow(d,a-1))/a;
2081  a += 1.0;
2082 
2083  if( a == 0 ) fIntegralMM[0] += b*log(x0/e0);
2084  else fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2085 
2086 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2087 
2088  x0 = fSplineEnergy[i - 1];
2089  x1 = fSplineEnergy[i - 2];
2090  y0 = fdNdxMM[i - 1];
2091  yy1 = fdNdxMM[i - 2];
2092 
2093  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2094  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2095 
2096  c = x1/x0;
2097  d = e0/x0;
2098  a = log10(yy1/y0)/log10(x1/x0);
2099  // b0 = log10(y0) - a*log10(x0);
2100  b = y0/pow(x0,a); // pow(10.,b0);
2101 
2102  a += 1.0;
2103  if( a == 0 ) result += b*log(e0/x0);
2104  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2105  a += 1.0;
2106 
2107  if( a == 0 ) fIntegralMM[0] += b*log(e0/x0);
2108  else fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2109 
2110  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2111  // <<b<<"; result = "<<result<<G4endl;
2112 
2113  return result;
2114 
2115 }
2116 
2118 //
2119 // Integration of Plasmon cross-section for the case of
2120 // passing across border between intervals
2121 
2122 G4double G4PAIxSection::SumOverBordPlasmon( G4int i ,
2123  G4double en0 )
2124 {
2125  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2126 
2127  e0 = en0;
2128  x0 = fSplineEnergy[i];
2129  x1 = fSplineEnergy[i+1];
2130  y0 = fdNdxPlasmon[i];
2131  yy1 = fdNdxPlasmon[i+1];
2132 
2133  c = x1/x0;
2134  d = e0/x0;
2135  a = log10(yy1/y0)/log10(c);
2136  if(a > 10.0) return 0.;
2137  // b0 = log10(y0) - a*log10(x0);
2138  b = y0/pow(x0,a); //pow(10.,b);
2139 
2140  a += 1.0;
2141  if( a == 0 ) result = b*log(x0/e0);
2142  else result = y0*(x0 - e0*pow(d,a-1))/a;
2143  a += 1.0;
2144 
2145  if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
2146  else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2147 
2148  x0 = fSplineEnergy[i - 1];
2149  x1 = fSplineEnergy[i - 2];
2150  y0 = fdNdxPlasmon[i - 1];
2151  yy1 = fdNdxPlasmon[i - 2];
2152 
2153  c = x1/x0;
2154  d = e0/x0;
2155  a = log10(yy1/y0)/log10(c);
2156  // b0 = log10(y0) - a*log10(x0);
2157  b = y0/pow(x0,a);// pow(10.,b0);
2158 
2159  a += 1.0;
2160  if( a == 0 ) result += b*log(e0/x0);
2161  else result += y0*(e0*pow(d,a-1) - x0)/a;
2162  a += 1.0;
2163 
2164  if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
2165  else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2166 
2167  return result;
2168 
2169 }
2170 
2172 //
2173 // Integration of resonance cross-section for the case of
2174 // passing across border between intervals
2175 
2176 G4double G4PAIxSection::SumOverBordResonance( G4int i ,
2177  G4double en0 )
2178 {
2179  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2180 
2181  e0 = en0;
2182  x0 = fSplineEnergy[i];
2183  x1 = fSplineEnergy[i+1];
2184  y0 = fdNdxResonance[i];
2185  yy1 = fdNdxResonance[i+1];
2186 
2187  c = x1/x0;
2188  d = e0/x0;
2189  a = log10(yy1/y0)/log10(c);
2190  if(a > 10.0) return 0.;
2191  // b0 = log10(y0) - a*log10(x0);
2192  b = y0/pow(x0,a); //pow(10.,b);
2193 
2194  a += 1.0;
2195  if( a == 0 ) result = b*log(x0/e0);
2196  else result = y0*(x0 - e0*pow(d,a-1))/a;
2197  a += 1.0;
2198 
2199  if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0);
2200  else fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2201 
2202  x0 = fSplineEnergy[i - 1];
2203  x1 = fSplineEnergy[i - 2];
2204  y0 = fdNdxResonance[i - 1];
2205  yy1 = fdNdxResonance[i - 2];
2206 
2207  c = x1/x0;
2208  d = e0/x0;
2209  a = log10(yy1/y0)/log10(c);
2210  // b0 = log10(y0) - a*log10(x0);
2211  b = y0/pow(x0,a);// pow(10.,b0);
2212 
2213  a += 1.0;
2214  if( a == 0 ) result += b*log(e0/x0);
2215  else result += y0*(e0*pow(d,a-1) - x0)/a;
2216  a += 1.0;
2217 
2218  if( a == 0 ) fIntegralResonance[0] += b*log(e0/x0);
2219  else fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2220 
2221  return result;
2222 
2223 }
2224 
2226 //
2227 // Returns random PAI-total energy loss over step
2228 
2229 G4double G4PAIxSection::GetStepEnergyLoss( G4double step )
2230 {
2231  G4long numOfCollisions;
2232  G4double meanNumber, loss = 0.0;
2233 
2234  // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl;
2235 
2236  meanNumber = fIntegralPAIxSection[1]*step;
2237  numOfCollisions = G4Poisson(meanNumber);
2238 
2239  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2240 
2241  while(numOfCollisions)
2242  {
2243  loss += GetEnergyTransfer();
2244  numOfCollisions--;
2245  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2246  }
2247  // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
2248 
2249  return loss;
2250 }
2251 
2253 //
2254 // Returns random PAI-total energy transfer in one collision
2255 
2256 G4double G4PAIxSection::GetEnergyTransfer()
2257 {
2258  G4int iTransfer ;
2259 
2260  G4double energyTransfer, position;
2261 
2262  position = fIntegralPAIxSection[1]*G4UniformRand();
2263 
2264  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2265  {
2266  if( position >= fIntegralPAIxSection[iTransfer] ) break;
2267  }
2268  if(iTransfer > fSplineNumber) iTransfer--;
2269 
2270  energyTransfer = fSplineEnergy[iTransfer];
2271 
2272  if(iTransfer > 1)
2273  {
2274  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2275  }
2276  return energyTransfer;
2277 }
2278 
2280 //
2281 // Returns random Cerenkov energy loss over step
2282 
2283 G4double G4PAIxSection::GetStepCerenkovLoss( G4double step )
2284 {
2285  G4long numOfCollisions;
2286  G4double meanNumber, loss = 0.0;
2287 
2288  // G4cout<<" G4PAIxSection::GetStepCerenkovLoss "<<G4endl;
2289 
2290  meanNumber = fIntegralCerenkov[1]*step;
2291  numOfCollisions = G4Poisson(meanNumber);
2292 
2293  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2294 
2295  while(numOfCollisions)
2296  {
2297  loss += GetCerenkovEnergyTransfer();
2298  numOfCollisions--;
2299  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2300  }
2301  // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2302 
2303  return loss;
2304 }
2305 
2307 //
2308 // Returns random MM-Cerenkov energy loss over step
2309 
2310 G4double G4PAIxSection::GetStepMMLoss( G4double step )
2311 {
2312  G4long numOfCollisions;
2313  G4double meanNumber, loss = 0.0;
2314 
2315  // G4cout<<" G4PAIxSection::GetStepMMLoss "<<G4endl;
2316 
2317  meanNumber = fIntegralMM[1]*step;
2318  numOfCollisions = G4Poisson(meanNumber);
2319 
2320  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2321 
2322  while(numOfCollisions)
2323  {
2324  loss += GetMMEnergyTransfer();
2325  numOfCollisions--;
2326  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2327  }
2328  // G4cout<<"PAI MM-Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2329 
2330  return loss;
2331 }
2332 
2334 //
2335 // Returns Cerenkov energy transfer in one collision
2336 
2337 G4double G4PAIxSection::GetCerenkovEnergyTransfer()
2338 {
2339  G4int iTransfer ;
2340 
2341  G4double energyTransfer, position;
2342 
2343  position = fIntegralCerenkov[1]*G4UniformRand();
2344 
2345  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2346  {
2347  if( position >= fIntegralCerenkov[iTransfer] ) break;
2348  }
2349  if(iTransfer > fSplineNumber) iTransfer--;
2350 
2351  energyTransfer = fSplineEnergy[iTransfer];
2352 
2353  if(iTransfer > 1)
2354  {
2355  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2356  }
2357  return energyTransfer;
2358 }
2359 
2361 //
2362 // Returns MM-Cerenkov energy transfer in one collision
2363 
2364 G4double G4PAIxSection::GetMMEnergyTransfer()
2365 {
2366  G4int iTransfer ;
2367 
2368  G4double energyTransfer, position;
2369 
2370  position = fIntegralMM[1]*G4UniformRand();
2371 
2372  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2373  {
2374  if( position >= fIntegralMM[iTransfer] ) break;
2375  }
2376  if(iTransfer > fSplineNumber) iTransfer--;
2377 
2378  energyTransfer = fSplineEnergy[iTransfer];
2379 
2380  if(iTransfer > 1)
2381  {
2382  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2383  }
2384  return energyTransfer;
2385 }
2386 
2388 //
2389 // Returns random plasmon energy loss over step
2390 
2391 G4double G4PAIxSection::GetStepPlasmonLoss( G4double step )
2392 {
2393  G4long numOfCollisions;
2394  G4double meanNumber, loss = 0.0;
2395 
2396  // G4cout<<" G4PAIxSection::GetStepPlasmonLoss "<<G4endl;
2397 
2398  meanNumber = fIntegralPlasmon[1]*step;
2399  numOfCollisions = G4Poisson(meanNumber);
2400 
2401  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2402 
2403  while(numOfCollisions)
2404  {
2405  loss += GetPlasmonEnergyTransfer();
2406  numOfCollisions--;
2407  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2408  }
2409  // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
2410 
2411  return loss;
2412 }
2413 
2415 //
2416 // Returns plasmon energy transfer in one collision
2417 
2418 G4double G4PAIxSection::GetPlasmonEnergyTransfer()
2419 {
2420  G4int iTransfer ;
2421 
2422  G4double energyTransfer, position;
2423 
2424  position = fIntegralPlasmon[1]*G4UniformRand();
2425 
2426  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2427  {
2428  if( position >= fIntegralPlasmon[iTransfer] ) break;
2429  }
2430  if(iTransfer > fSplineNumber) iTransfer--;
2431 
2432  energyTransfer = fSplineEnergy[iTransfer];
2433 
2434  if(iTransfer > 1)
2435  {
2436  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2437  }
2438  return energyTransfer;
2439 }
2440 
2442 //
2443 // Returns random resonance energy loss over step
2444 
2445 G4double G4PAIxSection::GetStepResonanceLoss( G4double step )
2446 {
2447  G4long numOfCollisions;
2448  G4double meanNumber, loss = 0.0;
2449 
2450  // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl;
2451 
2452  meanNumber = fIntegralResonance[1]*step;
2453  numOfCollisions = G4Poisson(meanNumber);
2454 
2455  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2456 
2457  while(numOfCollisions)
2458  {
2459  loss += GetResonanceEnergyTransfer();
2460  numOfCollisions--;
2461  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2462  }
2463  // G4cout<<"PAI resonance loss = "<<loss/keV<<" keV"<<G4endl;
2464 
2465  return loss;
2466 }
2467 
2468 
2470 //
2471 // Returns resonance energy transfer in one collision
2472 
2473 G4double G4PAIxSection::GetResonanceEnergyTransfer()
2474 {
2475  G4int iTransfer ;
2476 
2477  G4double energyTransfer, position;
2478 
2479  position = fIntegralResonance[1]*G4UniformRand();
2480 
2481  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2482  {
2483  if( position >= fIntegralResonance[iTransfer] ) break;
2484  }
2485  if(iTransfer > fSplineNumber) iTransfer--;
2486 
2487  energyTransfer = fSplineEnergy[iTransfer];
2488 
2489  if(iTransfer > 1)
2490  {
2491  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2492  }
2493  return energyTransfer;
2494 }
2495 
2496 
2498 //
2499 // Returns Rutherford energy transfer in one collision
2500 
2501 G4double G4PAIxSection::GetRutherfordEnergyTransfer()
2502 {
2503  G4int iTransfer ;
2504 
2505  G4double energyTransfer, position;
2506 
2507  position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand();
2508 
2509  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2510  {
2511  if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break;
2512  }
2513  if(iTransfer > fSplineNumber) iTransfer--;
2514 
2515  energyTransfer = fSplineEnergy[iTransfer];
2516 
2517  if(iTransfer > 1)
2518  {
2519  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2520  }
2521  return energyTransfer;
2522 }
2523 
2525 //
2526 
2527 void G4PAIxSection::CallError(G4int i, const G4String& methodName) const
2528 {
2529  G4String head = "G4PAIxSection::" + methodName + "()";
2530  G4ExceptionDescription ed;
2531  ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
2532  G4Exception(head,"pai001",FatalException,ed);
2533 }
2534 
2536 //
2537 // Init array of Lorentz factors
2538 //
2539 
2540 G4int G4PAIxSection::fNumberOfGammas = 111;
2541 
2542 const G4double G4PAIxSection::fLorentzFactor[112] = // fNumberOfGammas+1
2543 {
2544 0.0,
2545 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
2546 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
2547 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
2548 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
2549 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
2550 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
2551 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
2552 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
2553 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
2554 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
2555 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
2556 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
2557 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
2558 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
2559 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
2560 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
2561 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
2562 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
2563 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
2564 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
2565 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
2566 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
2567 1.065799e+05
2568 };
2569 
2571 //
2572 // The number of gamma for creation of spline (near ion-min , G ~ 4 )
2573 //
2574 
2575 const
2576 G4int G4PAIxSection::fRefGammaNumber = 29;
2577 
2578 
2579 //
2580 // end of G4PAIxSection implementation file
2581 //
2583 
G4SandiaTable::SandiaMixing
G4int SandiaMixing(G4int Z[], const G4double *fractionW, G4int el, G4int mi)
Definition: G4SandiaTable.cc:865
x
Float_t x
Definition: compare.C:6
G4PAIxSection::fError
static const G4double fError
Definition: G4PAIxSection.hh:223
G4PAIxSection::fRefGammaNumber
static const G4int fRefGammaNumber
Definition: G4PAIxSection.hh:229
G4PAIxSection::PAIdNdxPlasmon
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1412
G4PAIxSection::IntegralResonance
void IntegralResonance()
Definition: G4PAIxSection.cc:1613
G4PAIxSection::IntegralPlasmon
void IntegralPlasmon()
Definition: G4PAIxSection.cc:1586
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4PAIxSection::SumOverBordPlasmon
G4double SumOverBordPlasmon(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:2122
G4PAIxSection::PAIdNdxMM
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1362
G4PAIxSection::SumOverInterPlasmon
G4double SumOverInterPlasmon(G4int intervalNumber)
Definition: G4PAIxSection.cc:1795
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:593
keV
static constexpr double keV
Definition: G4SIunits.hh:216
G4PAIxSection::GetRutherfordEnergyTransfer
G4double GetRutherfordEnergyTransfer()
Definition: G4PAIxSection.cc:2501
G4PAIxSection::GetStepEnergyLoss
G4double GetStepEnergyLoss(G4double step)
Definition: G4PAIxSection.cc:2229
G4PAIxSection::GetResonanceEnergyTransfer
G4double GetResonanceEnergyTransfer()
Definition: G4PAIxSection.cc:2473
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:262
G4endl
#define G4endl
Definition: G4ios.hh:61
y
Float_t y
Definition: compare.C:6
CLHEP::hbarc
static constexpr double hbarc
Definition: PhysicalConstants.h:65
G4PAIxSection::IntegralPAIxSection
void IntegralPAIxSection()
Definition: G4PAIxSection.cc:1492
G4PAIxSection::GetStepMMLoss
G4double GetStepMMLoss(G4double step)
Definition: G4PAIxSection.cc:2310
G4PAIxSection::SumOverInterval
G4double SumOverInterval(G4int intervalNumber)
Definition: G4PAIxSection.cc:1640
G4PAIxSection::NormShift
void NormShift(G4double betaGammaSq)
Definition: G4PAIxSection.cc:875
G4PAIxSection::GetStepCerenkovLoss
G4double GetStepCerenkovLoss(G4double step)
Definition: G4PAIxSection.cc:2283
beta
Double_t beta
Definition: TG43_relative_dose.C:15
G4PAIxSection::fLorentzFactor
static const G4double fLorentzFactor[112]
Definition: G4PAIxSection.hh:226
G4PAIxSection::CallError
void CallError(G4int i, const G4String &methodName) const
Definition: G4PAIxSection.cc:2527
G4PAIxSection::Initialize
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
Definition: G4PAIxSection.cc:608
G4SandiaTable::GetMaxInterval
G4int GetMaxInterval() const
Definition: G4SandiaTable.cc:716
G4PAIxSection::GetEnergyTransfer
G4double GetEnergyTransfer()
Definition: G4PAIxSection.cc:2256
y2
Float_t y2[n_points_geant4]
Definition: compare.C:26
G4PAIxSection::ImPartDielectricConst
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1078
G4SandiaTable::GetLowerI1
G4bool GetLowerI1()
Definition: G4SandiaTable.hh:153
g
static constexpr double g
Definition: G4SIunits.hh:183
G4double
double G4double
Definition: G4Types.hh:76
G4PAIxSection::GetMMEnergyTransfer
G4double GetMMEnergyTransfer()
Definition: G4PAIxSection.cc:2364
position
#define position
Definition: xmlparse.cc:622
cm2
static constexpr double cm2
Definition: G4SIunits.hh:120
c2
TCanvas * c2
Definition: plot_hist.C:75
G4long
long G4long
Definition: G4Types.hh:80
pyG4EmCalculator::range
const G4ParticleDefinition const G4Material *G4double range
Definition: pyG4EmCalculator.cc:185
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4PAIxSection::RePartDielectricConst
G4double RePartDielectricConst(G4double energy)
Definition: G4PAIxSection.cc:1173
G4PAIxSection::SumOverBordResonance
G4double SumOverBordResonance(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:2176
G4PAIxSection::SumOverBordMM
G4double SumOverBordMM(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:2057
G4PAIxSection::fNumberOfGammas
static G4int fNumberOfGammas
Definition: G4PAIxSection.hh:225
G4PAIxSection::GetElectronRange
G4double GetElectronRange(G4double energy)
Definition: G4PAIxSection.cc:1131
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
energy
double energy
Definition: plottest35.C:25
d
Float_t d
Definition: advanced/microbeam/plot.C:236
G4SandiaTable::SandiaIntervals
G4int SandiaIntervals(G4int Z[], G4int el)
Definition: G4SandiaTable.cc:769
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4PAIxSection::SumOverIntervaldEdx
G4double SumOverIntervaldEdx(G4int intervalNumber)
Definition: G4PAIxSection.cc:1687
G4PAIxSection::SplainPAI
void SplainPAI(G4double betaGammaSq)
Definition: G4PAIxSection.cc:952
G4PAIxSection::SumOverBordCerenkov
G4double SumOverBordCerenkov(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1992
G4PAIxSection::DifPAIxSection
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1227
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int) const
Definition: G4SandiaTable.cc:1034
G4PAIxSection::fDelta
static const G4double fDelta
Definition: G4PAIxSection.hh:222
result
G4double G4ParticleHPJENDLHEData::G4double result
Definition: G4ParticleHPJENDLHEData.cc:257
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:43
G4PAIxSection::GetStepPlasmonLoss
G4double GetStepPlasmonLoss(G4double step)
Definition: G4PAIxSection.cc:2391
G4SandiaTable::GetPhotoAbsorpCof
G4double GetPhotoAbsorpCof(G4int i, G4int j) const
Definition: G4SandiaTable.cc:742
G4PAIxSection::SumOverInterResonance
G4double SumOverInterResonance(G4int intervalNumber)
Definition: G4PAIxSection.cc:1830
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4PAIxSection::RutherfordIntegral
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
Definition: G4PAIxSection.cc:1057
G4int
int G4int
Definition: G4Types.hh:78
G4PAIxSection::PAIdNdxResonance
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1452
G4PAIxSection::GetPlasmonEnergyTransfer
G4double GetPlasmonEnergyTransfer()
Definition: G4PAIxSection.cc:2418
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
G4cout
G4GLOB_DLL std::ostream G4cout
G4PAIxSection::SumOverInterCerenkov
G4double SumOverInterCerenkov(G4int intervalNumber)
Definition: G4PAIxSection.cc:1721
G4PAIxSection::ComputeLowEnergyCof
void ComputeLowEnergyCof()
Definition: G4PAIxSection.cc:781
G4PAIxSection::fMaxSplineSize
static const G4int fMaxSplineSize
Definition: G4PAIxSection.hh:256
G4Poisson
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4PAIxSection::PAIdNdxCerenkov
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1304
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4PAIxSection::GetCerenkovEnergyTransfer
G4double GetCerenkovEnergyTransfer()
Definition: G4PAIxSection.cc:2337
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4PAIxSection::SumOverBorderdEdx
G4double SumOverBorderdEdx(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1937
CLHEP::fine_structure_const
static constexpr double fine_structure_const
Definition: PhysicalConstants.h:96
x2
Float_t x2[n_points_geant4]
Definition: compare.C:26
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4PAIxSection::SumOverInterMM
G4double SumOverInterMM(G4int intervalNumber)
Definition: G4PAIxSection.cc:1757
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:203
G4PAIxSection::GetLorentzFactor
G4double GetLorentzFactor(G4int i) const
Definition: G4PAIxSection.cc:599
G4PAIxSection::IntegralCerenkov
void IntegralCerenkov()
Definition: G4PAIxSection.cc:1524
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:218
G4PAIxSection::GetStepResonanceLoss
G4double GetStepResonanceLoss(G4double step)
Definition: G4PAIxSection.cc:2445
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
G4PAIxSection::SumOverBorder
G4double SumOverBorder(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1864
G4PAIxSection::GetPhotonRange
G4double GetPhotonRange(G4double energy)
Definition: G4PAIxSection.cc:1098
G4Material::GetDensity
G4double GetDensity() const
Definition: G4Material.hh:181
다음에 의해 생성됨 :   doxygen 1.8.5