ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIxSection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
28 //
29 // G4PAIxSection.cc -- class implementation file
30 //
31 // GEANT 4 class implementation file
32 //
33 // For information related to this code, please, contact
34 // the Geant4 Collaboration.
35 //
36 // R&D: Vladimir.Grichine@cern.ch
37 //
38 // History:
39 //
40 // 13.05.03 V. Grichine, bug fixed for maxEnergyTransfer > max interval energy
41 // 28.05.01 V.Ivanchenko minor changes to provide ANSI -wall compilation
42 // 17.05.01 V. Grichine, low energy extension down to 10*keV of proton
43 // 20.11.98 adapted to a new Material/SandiaTable interface, mma
44 // 11.06.97 V. Grichine, 1st version
45 //
46 
47 #include "G4PAIxSection.hh"
48 
49 #include "globals.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4ios.hh"
53 #include "G4Poisson.hh"
54 #include "G4Material.hh"
55 #include "G4MaterialCutsCouple.hh"
56 #include "G4SandiaTable.hh"
57 
58 using namespace std;
59 
60 /* ******************************************************************
61 
62 // Init array of Lorentz factors
63 
64 const G4double G4PAIxSection::fLorentzFactor[22] =
65 {
66  0.0 , 1.1 , 1.2 , 1.3 , 1.5 , 1.8 , 2.0 ,
67  2.5 , 3.0 , 4.0 , 7.0 , 10.0 , 20.0 , 40.0 ,
68  70.0 , 100.0 , 300.0 , 600.0 , 1000.0 , 3000.0 ,
69  10000.0 , 50000.0
70 };
71 
72 const G4int G4PAIxSection::
73 fRefGammaNumber = 29; // The number of gamma for creation of
74  // spline (9)
75 
76 ***************************************************************** */
77 
78 // Local class constants
79 
80 const G4double G4PAIxSection::fDelta = 0.005; // 0.005 energy shift from interval border
81 const G4double G4PAIxSection::fError = 0.005; // 0.005 error in lin-log approximation
82 
83 const G4int G4PAIxSection::fMaxSplineSize = 1000; // Max size of output spline
84  // arrays
86 //
87 // Constructor
88 //
89 
90 G4PAIxSection::G4PAIxSection()
91 {
92  fSandia = 0;
93  fMatSandiaMatrix = 0;
94  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
95  fIntervalNumber = fSplineNumber = 0;
96  fVerbose = 0;
97 
98  fSplineEnergy = G4DataVector(fMaxSplineSize,0.0);
99  fRePartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
100  fImPartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
101  fIntegralTerm = G4DataVector(fMaxSplineSize,0.0);
102  fDifPAIxSection = G4DataVector(fMaxSplineSize,0.0);
103  fdNdxCerenkov = G4DataVector(fMaxSplineSize,0.0);
104  fdNdxPlasmon = G4DataVector(fMaxSplineSize,0.0);
105  fdNdxMM = G4DataVector(fMaxSplineSize,0.0);
106  fdNdxResonance = G4DataVector(fMaxSplineSize,0.0);
107  fIntegralPAIxSection = G4DataVector(fMaxSplineSize,0.0);
108  fIntegralPAIdEdx = G4DataVector(fMaxSplineSize,0.0);
109  fIntegralCerenkov = G4DataVector(fMaxSplineSize,0.0);
110  fIntegralPlasmon = G4DataVector(fMaxSplineSize,0.0);
111  fIntegralMM = G4DataVector(fMaxSplineSize,0.0);
112  fIntegralResonance = G4DataVector(fMaxSplineSize,0.0);
113 
114  fMaterialIndex = 0;
115 
116  for( G4int i = 0; i < 500; ++i )
117  {
118  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
119  }
120 }
121 
123 //
124 // Constructor
125 //
126 
127 G4PAIxSection::G4PAIxSection(G4MaterialCutsCouple* matCC)
128 {
129  fDensity = matCC->GetMaterial()->GetDensity();
130  G4int matIndex = matCC->GetMaterial()->GetIndex();
131  fMaterialIndex = matIndex;
132 
133  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
134  fSandia = (*theMaterialTable)[matIndex]->GetSandiaTable();
135 
136  fVerbose = 0;
137 
138  G4int i, j;
139  fMatSandiaMatrix = new G4OrderedTable();
140 
141  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
142  {
143  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
144  }
145  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
146  {
147  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
148 
149  for(j = 1; j < 5; j++)
150  {
151  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
152  }
153  }
154  ComputeLowEnergyCof();
155  // fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0;
156 }
157 
159 
160 G4PAIxSection::G4PAIxSection(G4int materialIndex,
161  G4double maxEnergyTransfer)
162 {
163  fSandia = 0;
164  fMatSandiaMatrix = 0;
165  fVerbose = 0;
166  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
167  G4int i, j;
168 
169  fMaterialIndex = materialIndex;
170  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
171  fElectronDensity = (*theMaterialTable)[materialIndex]->
172  GetElectronDensity();
173  fIntervalNumber = (*theMaterialTable)[materialIndex]->
174  GetSandiaTable()->GetMatNbOfIntervals();
175  fIntervalNumber--;
176  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
177 
178  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
179  fA1 = G4DataVector(fIntervalNumber+2,0.0);
180  fA2 = G4DataVector(fIntervalNumber+2,0.0);
181  fA3 = G4DataVector(fIntervalNumber+2,0.0);
182  fA4 = G4DataVector(fIntervalNumber+2,0.0);
183 
184  for(i = 1; i <= fIntervalNumber; i++ )
185  {
186  if(((*theMaterialTable)[materialIndex]->
187  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
188  i > fIntervalNumber )
189  {
190  fEnergyInterval[i] = maxEnergyTransfer;
191  fIntervalNumber = i;
192  break;
193  }
194  fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
195  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
196  fA1[i] = (*theMaterialTable)[materialIndex]->
197  GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
198  fA2[i] = (*theMaterialTable)[materialIndex]->
199  GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
200  fA3[i] = (*theMaterialTable)[materialIndex]->
201  GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
202  fA4[i] = (*theMaterialTable)[materialIndex]->
203  GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
204  // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
205  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
206  }
207  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
208  {
209  fIntervalNumber++;
210  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
211  }
212 
213  // Now checking, if two borders are too close together
214 
215  for(i=1;i<fIntervalNumber;i++)
216  {
217  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
218  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
219  {
220  continue;
221  }
222  else
223  {
224  for(j=i;j<fIntervalNumber;j++)
225  {
226  fEnergyInterval[j] = fEnergyInterval[j+1];
227  fA1[j] = fA1[j+1];
228  fA2[j] = fA2[j+1];
229  fA3[j] = fA3[j+1];
230  fA4[j] = fA4[j+1];
231  }
232  fIntervalNumber--;
233  i--;
234  }
235  }
236 
237 
238  /* *********************************
239 
240  fSplineEnergy = new G4double[fMaxSplineSize];
241  fRePartDielectricConst = new G4double[fMaxSplineSize];
242  fImPartDielectricConst = new G4double[fMaxSplineSize];
243  fIntegralTerm = new G4double[fMaxSplineSize];
244  fDifPAIxSection = new G4double[fMaxSplineSize];
245  fIntegralPAIxSection = new G4double[fMaxSplineSize];
246 
247  for(i=0;i<fMaxSplineSize;i++)
248  {
249  fSplineEnergy[i] = 0.0;
250  fRePartDielectricConst[i] = 0.0;
251  fImPartDielectricConst[i] = 0.0;
252  fIntegralTerm[i] = 0.0;
253  fDifPAIxSection[i] = 0.0;
254  fIntegralPAIxSection[i] = 0.0;
255  }
256  ************************************************** */
257  ComputeLowEnergyCof();
258  InitPAI(); // create arrays allocated above
259  /*
260  delete[] fEnergyInterval;
261  delete[] fA1;
262  delete[] fA2;
263  delete[] fA3;
264  delete[] fA4;
265  */
266 }
267 
269 //
270 // Constructor called from G4PAIPhotonModel !!!
271 
272 G4PAIxSection::G4PAIxSection( G4int materialIndex,
273  G4double maxEnergyTransfer,
274  G4double betaGammaSq,
275  G4double** photoAbsCof,
276  G4int intNumber )
277 {
278  fSandia = 0;
279  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
280  fIntervalNumber = fSplineNumber = 0;
281  fVerbose = 0;
282 
283  fSplineEnergy = G4DataVector(500,0.0);
284  fRePartDielectricConst = G4DataVector(500,0.0);
285  fImPartDielectricConst = G4DataVector(500,0.0);
286  fIntegralTerm = G4DataVector(500,0.0);
287  fDifPAIxSection = G4DataVector(500,0.0);
288  fdNdxCerenkov = G4DataVector(500,0.0);
289  fdNdxPlasmon = G4DataVector(500,0.0);
290  fdNdxMM = G4DataVector(500,0.0);
291  fdNdxResonance = G4DataVector(500,0.0);
292  fIntegralPAIxSection = G4DataVector(500,0.0);
293  fIntegralPAIdEdx = G4DataVector(500,0.0);
294  fIntegralCerenkov = G4DataVector(500,0.0);
295  fIntegralPlasmon = G4DataVector(500,0.0);
296  fIntegralMM = G4DataVector(500,0.0);
297  fIntegralResonance = G4DataVector(500,0.0);
298 
299  for( G4int i = 0; i < 500; ++i )
300  {
301  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
302  }
303 
304  fSandia = 0;
305  fMatSandiaMatrix = 0;
306  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
307  G4int i, j;
308 
309  fMaterialIndex = materialIndex;
310  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
311  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
312 
313  fIntervalNumber = intNumber;
314  fIntervalNumber--;
315  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
316 
317  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
318  fA1 = G4DataVector(fIntervalNumber+2,0.0);
319  fA2 = G4DataVector(fIntervalNumber+2,0.0);
320  fA3 = G4DataVector(fIntervalNumber+2,0.0);
321  fA4 = G4DataVector(fIntervalNumber+2,0.0);
322 
323 
324  /*
325  fEnergyInterval = new G4double[fIntervalNumber+2];
326  fA1 = new G4double[fIntervalNumber+2];
327  fA2 = new G4double[fIntervalNumber+2];
328  fA3 = new G4double[fIntervalNumber+2];
329  fA4 = new G4double[fIntervalNumber+2];
330  */
331  for( i = 1; i <= fIntervalNumber; i++ )
332  {
333  if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
334  i > fIntervalNumber )
335  {
336  fEnergyInterval[i] = maxEnergyTransfer;
337  fIntervalNumber = i;
338  break;
339  }
340  fEnergyInterval[i] = photoAbsCof[i-1][0];
341  fA1[i] = photoAbsCof[i-1][1];
342  fA2[i] = photoAbsCof[i-1][2];
343  fA3[i] = photoAbsCof[i-1][3];
344  fA4[i] = photoAbsCof[i-1][4];
345  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
346  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
347  }
348  // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
349 
350  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
351  {
352  fIntervalNumber++;
353  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
354  }
355  // G4cout<<"after check of max energy transfer"<<G4endl;
356 
357  for( i = 1; i <= fIntervalNumber; i++ )
358  {
359  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
360  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
361  }
362  // Now checking, if two borders are too close together
363 
364  for( i = 1; i < fIntervalNumber; i++ )
365  {
366  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
367  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
368  {
369  continue;
370  }
371  else
372  {
373  for(j=i;j<fIntervalNumber;j++)
374  {
375  fEnergyInterval[j] = fEnergyInterval[j+1];
376  fA1[j] = fA1[j+1];
377  fA2[j] = fA2[j+1];
378  fA3[j] = fA3[j+1];
379  fA4[j] = fA4[j+1];
380  }
381  fIntervalNumber--;
382  i--;
383  }
384  }
385  // G4cout<<"after check of close borders"<<G4endl;
386 
387  for( i = 1; i <= fIntervalNumber; i++ )
388  {
389  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
390  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
391  }
392 
393  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
394 
395  ComputeLowEnergyCof();
396  G4double betaGammaSqRef =
397  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
398 
399  NormShift(betaGammaSqRef);
400  SplainPAI(betaGammaSqRef);
401 
402  // Preparation of integral PAI cross section for input betaGammaSq
403 
404  for(i = 1; i <= fSplineNumber; i++)
405  {
406  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
407  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
408  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
409  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
410  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
411 
412  // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
413  // <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
414  }
415  IntegralCerenkov();
416  IntegralMM();
417  IntegralPlasmon();
418  IntegralResonance();
419  IntegralPAIxSection();
420  /*
421  delete[] fEnergyInterval;
422  delete[] fA1;
423  delete[] fA2;
424  delete[] fA3;
425  delete[] fA4;
426  */
427 }
428 
430 //
431 // Test Constructor with beta*gamma square value
432 
433 G4PAIxSection::G4PAIxSection( G4int materialIndex,
434  G4double maxEnergyTransfer,
435  G4double betaGammaSq )
436 {
437  fSandia = 0;
438  fMatSandiaMatrix = 0;
439  fVerbose = 0;
440  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
441 
442  G4int i, j, numberOfElements;
443 
444  fMaterialIndex = materialIndex;
445  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
446  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
447  numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements();
448 
449  G4int* thisMaterialZ = new G4int[numberOfElements];
450 
451  for( i = 0; i < numberOfElements; i++ )
452  {
453  thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
454  GetElement(i)->GetZ();
455  }
456  // fSandia = new G4SandiaTable(materialIndex);
457  fSandia = (*theMaterialTable)[materialIndex]->GetSandiaTable();
458  G4SandiaTable thisMaterialSandiaTable(materialIndex);
459  fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals(thisMaterialZ,
460  numberOfElements);
461  fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
462  ( thisMaterialZ ,
463  (*theMaterialTable)[materialIndex]->GetFractionVector() ,
464  numberOfElements,fIntervalNumber);
465 
466  fIntervalNumber--;
467 
468  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
469  fA1 = G4DataVector(fIntervalNumber+2,0.0);
470  fA2 = G4DataVector(fIntervalNumber+2,0.0);
471  fA3 = G4DataVector(fIntervalNumber+2,0.0);
472  fA4 = G4DataVector(fIntervalNumber+2,0.0);
473 
474  /*
475  fEnergyInterval = new G4double[fIntervalNumber+2];
476  fA1 = new G4double[fIntervalNumber+2];
477  fA2 = new G4double[fIntervalNumber+2];
478  fA3 = new G4double[fIntervalNumber+2];
479  fA4 = new G4double[fIntervalNumber+2];
480  */
481  for( i = 1; i <= fIntervalNumber; i++ )
482  {
483  if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
484  i > fIntervalNumber)
485  {
486  fEnergyInterval[i] = maxEnergyTransfer;
487  fIntervalNumber = i;
488  break;
489  }
490  fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0);
491  fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity;
492  fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity;
493  fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity;
494  fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity;
495 
496  }
497  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
498  {
499  fIntervalNumber++;
500  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
501  fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
502  fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
503  fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
504  fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
505  }
506  for(i=1;i<=fIntervalNumber;i++)
507  {
508  // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
509  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
510  }
511  // Now checking, if two borders are too close together
512 
513  for( i = 1; i < fIntervalNumber; i++ )
514  {
515  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
516  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
517  {
518  continue;
519  }
520  else
521  {
522  for( j = i; j < fIntervalNumber; j++ )
523  {
524  fEnergyInterval[j] = fEnergyInterval[j+1];
525  fA1[j] = fA1[j+1];
526  fA2[j] = fA2[j+1];
527  fA3[j] = fA3[j+1];
528  fA4[j] = fA4[j+1];
529  }
530  fIntervalNumber--;
531  i--;
532  }
533  }
534 
535  /* *********************************
536  fSplineEnergy = new G4double[fMaxSplineSize];
537  fRePartDielectricConst = new G4double[fMaxSplineSize];
538  fImPartDielectricConst = new G4double[fMaxSplineSize];
539  fIntegralTerm = new G4double[fMaxSplineSize];
540  fDifPAIxSection = new G4double[fMaxSplineSize];
541  fIntegralPAIxSection = new G4double[fMaxSplineSize];
542 
543  for(i=0;i<fMaxSplineSize;i++)
544  {
545  fSplineEnergy[i] = 0.0;
546  fRePartDielectricConst[i] = 0.0;
547  fImPartDielectricConst[i] = 0.0;
548  fIntegralTerm[i] = 0.0;
549  fDifPAIxSection[i] = 0.0;
550  fIntegralPAIxSection[i] = 0.0;
551  }
552  */
553 
554  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
555 
556  ComputeLowEnergyCof();
557  G4double betaGammaSqRef =
558  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
559 
560  NormShift(betaGammaSqRef);
561  SplainPAI(betaGammaSqRef);
562 
563  // Preparation of integral PAI cross section for input betaGammaSq
564 
565  for(i = 1; i <= fSplineNumber; i++)
566  {
567  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
568  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
569  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
570  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
571  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
572  }
573  IntegralPAIxSection();
574  IntegralCerenkov();
575  IntegralMM();
576  IntegralPlasmon();
577  IntegralResonance();
578 }
579 
581 //
582 // Destructor
583 
584 G4PAIxSection::~G4PAIxSection()
585 {
586  /* ************************
587  delete[] fSplineEnergy ;
588  delete[] fRePartDielectricConst;
589  delete[] fImPartDielectricConst;
590  delete[] fIntegralTerm ;
591  delete[] fDifPAIxSection ;
592  delete[] fIntegralPAIxSection ;
593  */
594  delete fMatSandiaMatrix;
595 }
596 
597 G4double G4PAIxSection::GetLorentzFactor(G4int j) const
598 {
599  return fLorentzFactor[j];
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 
746 void G4PAIxSection::ComputeLowEnergyCof(const G4Material* material)
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 
779 void G4PAIxSection::ComputeLowEnergyCof()
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 
813 void G4PAIxSection::InitPAI()
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 
873 void G4PAIxSection::NormShift(G4double betaGammaSq)
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 
950 void G4PAIxSection::SplainPAI(G4double betaGammaSq)
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 
1055 G4double G4PAIxSection::RutherfordIntegral( G4int k,
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 
1076 G4double G4PAIxSection::ImPartDielectricConst( G4int k ,
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 
1096 G4double G4PAIxSection::GetPhotonRange( G4double energy1 )
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 
1129 G4double G4PAIxSection::GetElectronRange( G4double energy )
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 
1171 G4double G4PAIxSection::RePartDielectricConst(G4double enb)
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 
1225 G4double G4PAIxSection::DifPAIxSection( G4int i ,
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 
1302 G4double G4PAIxSection::PAIdNdxCerenkov( G4int i ,
1303  G4double betaGammaSq )
1304 {
1305  G4double logarithm, x3, x5, argument, modul2, dNdxC;
1306  G4double be2, betaBohr2, cofBetaBohr;
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 
1360 G4double G4PAIxSection::PAIdNdxMM( G4int i ,
1361  G4double betaGammaSq )
1362 {
1363  G4double logarithm, x3, x5, argument, dNdxC;
1364  G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
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 
1410 G4double G4PAIxSection::PAIdNdxPlasmon( G4int i ,
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 
1450 G4double G4PAIxSection::PAIdNdxResonance( G4int i ,
1451  G4double betaGammaSq )
1452 {
1453  G4double resonance, modul2, dNdxP;
1454  G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
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 
1490 void G4PAIxSection::IntegralPAIxSection()
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 
1522 void G4PAIxSection::IntegralCerenkov()
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 
1553 void G4PAIxSection::IntegralMM()
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 
1584 void G4PAIxSection::IntegralPlasmon()
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 
1611 void G4PAIxSection::IntegralResonance()
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 
1638 G4double G4PAIxSection::SumOverInterval( G4int i )
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::abs(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::abs(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 
1685 G4double G4PAIxSection::SumOverIntervaldEdx( G4int i )
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 
1719 G4double G4PAIxSection::SumOverInterCerenkov( G4int i )
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 
1755 G4double G4PAIxSection::SumOverInterMM( G4int i )
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 
1793 G4double G4PAIxSection::SumOverInterPlasmon( G4int i )
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 
1828 G4double G4PAIxSection::SumOverInterResonance( G4int i )
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 
1862 G4double G4PAIxSection::SumOverBorder( G4int i ,
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::abs(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::abs(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::abs(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::abs(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 
1935 G4double G4PAIxSection::SumOverBorderdEdx( G4int i ,
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 
1990 G4double G4PAIxSection::SumOverBordCerenkov( G4int i ,
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 
2055 G4double G4PAIxSection::SumOverBordMM( G4int i ,
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 
2120 G4double G4PAIxSection::SumOverBordPlasmon( G4int i ,
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 
2174 G4double G4PAIxSection::SumOverBordResonance( G4int i ,
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 
2227 G4double G4PAIxSection::GetStepEnergyLoss( G4double step )
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 
2254 G4double G4PAIxSection::GetEnergyTransfer()
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 
2281 G4double G4PAIxSection::GetStepCerenkovLoss( G4double step )
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 
2308 G4double G4PAIxSection::GetStepMMLoss( G4double step )
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 
2335 G4double G4PAIxSection::GetCerenkovEnergyTransfer()
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 
2362 G4double G4PAIxSection::GetMMEnergyTransfer()
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 
2389 G4double G4PAIxSection::GetStepPlasmonLoss( G4double step )
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 
2416 G4double G4PAIxSection::GetPlasmonEnergyTransfer()
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 
2443 G4double G4PAIxSection::GetStepResonanceLoss( G4double step )
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 
2471 G4double G4PAIxSection::GetResonanceEnergyTransfer()
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 
2499 G4double G4PAIxSection::GetRutherfordEnergyTransfer()
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 + "()";
2528  G4ExceptionDescription ed;
2529  ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
2530  G4Exception(head,"pai001",FatalException,ed);
2531 }
2532 
2534 //
2535 // Init array of Lorentz factors
2536 //
2537 
2538 G4int G4PAIxSection::fNumberOfGammas = 111;
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
2574 G4int G4PAIxSection::fRefGammaNumber = 29;
2575 
2576 
2577 //
2578 // end of G4PAIxSection implementation file
2579 //
2581