Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIxSection.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4PAIxSection.cc 91726 2015-08-03 15:41:36Z gcosmo $
28 // GEANT4 tag $Name: geant4-09-03-ref-06 $
29 //
30 //
31 // G4PAIxSection.cc -- class implementation file
32 //
33 // GEANT 4 class implementation file
34 //
35 // For information related to this code, please, contact
36 // the Geant4 Collaboration.
37 //
38 // R&D: Vladimir.Grichine@cern.ch
39 //
40 // History:
41 //
42 // 13.05.03 V. Grichine, bug fixed for maxEnergyTransfer > max interval energy
43 // 28.05.01 V.Ivanchenko minor changes to provide ANSI -wall compilation
44 // 17.05.01 V. Grichine, low energy extension down to 10*keV of proton
45 // 20.11.98 adapted to a new Material/SandiaTable interface, mma
46 // 11.06.97 V. Grichine, 1st version
47 //
48 
49 #include "G4PAIxSection.hh"
50 
51 #include "globals.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4ios.hh"
55 #include "G4Poisson.hh"
56 #include "G4Material.hh"
57 #include "G4MaterialCutsCouple.hh"
58 #include "G4SandiaTable.hh"
59 
60 using namespace std;
61 
62 /* ******************************************************************
63 
64 // Init array of Lorentz factors
65 
66 const G4double G4PAIxSection::fLorentzFactor[22] =
67 {
68  0.0 , 1.1 , 1.2 , 1.3 , 1.5 , 1.8 , 2.0 ,
69  2.5 , 3.0 , 4.0 , 7.0 , 10.0 , 20.0 , 40.0 ,
70  70.0 , 100.0 , 300.0 , 600.0 , 1000.0 , 3000.0 ,
71  10000.0 , 50000.0
72 };
73 
74 const G4int G4PAIxSection::
75 fRefGammaNumber = 29; // The number of gamma for creation of
76  // spline (9)
77 
78 ***************************************************************** */
79 
80 // Local class constants
81 
82 const G4double G4PAIxSection::fDelta = 0.005; // 0.005 energy shift from interval border
83 const G4double G4PAIxSection::fError = 0.005; // 0.005 error in lin-log approximation
84 
85 const G4int G4PAIxSection::fMaxSplineSize = 1000; // Max size of output spline
86  // arrays
88 //
89 // Constructor
90 //
91 
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 
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::fabs(a) < 1.e-6 )
1662  {
1663  result = b*log(x1/x0);
1664  }
1665  else
1666  {
1667  result = y0*(x1*pow(c,a-1) - x0)/a;
1668  }
1669  a += 1.;
1670  if( std::fabs(a) < 1.e-6 )
1671  {
1672  fIntegralPAIxSection[0] += b*log(x1/x0);
1673  }
1674  else
1675  {
1676  fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1677  }
1678  if(fVerbose>0) G4cout<<"SumOverInterval, result = "<<result<<G4endl;
1679  return result;
1680 
1681 } // end of SumOverInterval
1682 
1684 
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::fabs(a) < 1.e-6 )
1885  {
1886  result = b*log(x0/e0);
1887  }
1888  else
1889  {
1890  result = y0*(x0 - e0*pow(d,a-1))/a;
1891  }
1892  a += 1.;
1893  if( std::fabs(a) < 1.e-6 )
1894  {
1895  fIntegralPAIxSection[0] += b*log(x0/e0);
1896  }
1897  else
1898  {
1899  fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1900  }
1901  x0 = fSplineEnergy[i - 1];
1902  x1 = fSplineEnergy[i - 2];
1903  y0 = fDifPAIxSection[i - 1];
1904  yy1 = fDifPAIxSection[i - 2];
1905 
1906  //c = x1/x0;
1907  d = e0/x0;
1908  a = log10(yy1/y0)/log10(x1/x0);
1909  // b0 = log10(y0) - a*log10(x0);
1910  b = y0/pow(x0,a);
1911  a += 1.;
1912  if( std::fabs(a) < 1.e-6 )
1913  {
1914  result += b*log(e0/x0);
1915  }
1916  else
1917  {
1918  result += y0*(e0*pow(d,a-1) - x0)/a;
1919  }
1920  a += 1.;
1921  if( std::fabs(a) < 1.e-6 )
1922  {
1923  fIntegralPAIxSection[0] += b*log(e0/x0);
1924  }
1925  else
1926  {
1927  fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1928  }
1929  return result;
1930 
1931 }
1932 
1934 
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 
result
G4double G4ParticleHPJENDLHEData::G4double result
Definition: G4ParticleHPJENDLHEData.cc:257
G4PAIxSection::IntegralPAIxSection
void IntegralPAIxSection()
Definition: G4PAIxSection.cc:1490
G4PAIxSection::SplainPAI
void SplainPAI(G4double betaGammaSq)
Definition: G4PAIxSection.cc:950
G4PAIxSection::GetRutherfordEnergyTransfer
G4double GetRutherfordEnergyTransfer()
Definition: G4PAIxSection.cc:2499
G4PAIxSection::PAIdNdxResonance
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1450
G4PAIxSection::IntegralCerenkov
void IntegralCerenkov()
Definition: G4PAIxSection.cc:1522
c2
static c2_factory< G4double > c2
Definition: G4ScreenedNuclearRecoil.cc:126
G4Poisson
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4PAIxSection::GetPlasmonEnergyTransfer
G4double GetPlasmonEnergyTransfer()
Definition: G4PAIxSection.cc:2416
G4PAIxSection::IntegralPlasmon
void IntegralPlasmon()
Definition: G4PAIxSection.cc:1584
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4PAIxSection::IntegralResonance
void IntegralResonance()
Definition: G4PAIxSection.cc:1611
cm2
static constexpr double cm2
Definition: G4SIunits.hh:120
G4PAIxSection::GetStepResonanceLoss
G4double GetStepResonanceLoss(G4double step)
Definition: G4PAIxSection.cc:2443
G4PAIxSection::PAIdNdxCerenkov
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1302
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4MaterialCutsCouple.hh
G4PAIxSection::SumOverBordPlasmon
G4double SumOverBordPlasmon(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:2120
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:262
G4PAIxSection::SumOverInterMM
G4double SumOverInterMM(G4int intervalNumber)
Definition: G4PAIxSection.cc:1755
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:41
G4Material::GetDensity
G4double GetDensity() const
Definition: G4Material.hh:180
G4long
long G4long
Definition: G4Types.hh:80
G4PAIxSection::NormShift
void NormShift(G4double betaGammaSq)
Definition: G4PAIxSection.cc:873
G4SandiaTable::GetMaxInterval
G4int GetMaxInterval() const
Definition: G4SandiaTable.cc:716
G4PAIxSection::PAIdNdxPlasmon
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1410
G4PAIxSection::SumOverBordResonance
G4double SumOverBordResonance(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:2174
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int) const
Definition: G4SandiaTable.cc:1034
G4PAIxSection::SumOverInterPlasmon
G4double SumOverInterPlasmon(G4int intervalNumber)
Definition: G4PAIxSection.cc:1793
test.x
tuple x
Definition: test.py:50
G4PAIxSection::SumOverInterval
G4double SumOverInterval(G4int intervalNumber)
Definition: G4PAIxSection.cc:1638
G4PAIxSection::RePartDielectricConst
G4double RePartDielectricConst(G4double energy)
Definition: G4PAIxSection.cc:1171
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:202
G4PAIxSection::ComputeLowEnergyCof
void ComputeLowEnergyCof()
Definition: G4PAIxSection.cc:779
G4int
int G4int
Definition: G4Types.hh:78
G4SandiaTable::SandiaMixing
G4int SandiaMixing(G4int Z[], const G4double *fractionW, G4int el, G4int mi)
Definition: G4SandiaTable.cc:865
eplot.material
string material
Definition: eplot.py:19
G4PAIxSection.hh
g
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
test.b
tuple b
Definition: test.py:12
G4PAIxSection::GetResonanceEnergyTransfer
G4double GetResonanceEnergyTransfer()
Definition: G4PAIxSection.cc:2471
position
#define position
Definition: xmlparse.cc:622
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
pyG4EmCalculator::range
const G4ParticleDefinition const G4Material *G4double range
Definition: pyG4EmCalculator.cc:185
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4PAIxSection::ImPartDielectricConst
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1076
G4SandiaTable::GetLowerI1
G4bool GetLowerI1()
Definition: G4SandiaTable.hh:153
G4PAIxSection::GetPhotonRange
G4double GetPhotonRange(G4double energy)
Definition: G4PAIxSection.cc:1096
G4PAIxSection::SumOverBordMM
G4double SumOverBordMM(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:2055
G4Poisson.hh
G4Material.hh
G4PAIxSection::GetCerenkovEnergyTransfer
G4double GetCerenkovEnergyTransfer()
Definition: G4PAIxSection.cc:2335
eV
static constexpr double eV
Definition: G4SIunits.hh:215
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4PAIxSection::SumOverInterResonance
G4double SumOverInterResonance(G4int intervalNumber)
Definition: G4PAIxSection.cc:1828
G4SandiaTable::SandiaIntervals
G4int SandiaIntervals(G4int Z[], G4int el)
Definition: G4SandiaTable.cc:769
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhysicalConstants.hh
G4PAIxSection::GetStepPlasmonLoss
G4double GetStepPlasmonLoss(G4double step)
Definition: G4PAIxSection.cc:2389
G4PAIxSection::PAIdNdxMM
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1360
G4PAIxSection::GetStepEnergyLoss
G4double GetStepEnergyLoss(G4double step)
Definition: G4PAIxSection.cc:2227
G4PAIxSection::SumOverBordCerenkov
G4double SumOverBordCerenkov(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1990
G4PAIxSection::GetMMEnergyTransfer
G4double GetMMEnergyTransfer()
Definition: G4PAIxSection.cc:2362
G4PAIxSection::SumOverBorder
G4double SumOverBorder(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1862
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
G4PAIxSection::Initialize
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
Definition: G4PAIxSection.cc:606
G4PAIxSection::SumOverInterCerenkov
G4double SumOverInterCerenkov(G4int intervalNumber)
Definition: G4PAIxSection.cc:1719
G4ios.hh
FatalException
Definition: G4ExceptionSeverity.hh:60
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
G4SandiaTable::GetPhotoAbsorpCof
G4double GetPhotoAbsorpCof(G4int i, G4int j) const
Definition: G4SandiaTable.cc:742
G4endl
#define G4endl
Definition: G4ios.hh:61
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4SandiaTable.hh
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
test.c
tuple c
Definition: test.py:13
G4PAIxSection::GetEnergyTransfer
G4double GetEnergyTransfer()
Definition: G4PAIxSection.cc:2254
G4PAIxSection::GetStepMMLoss
G4double GetStepMMLoss(G4double step)
Definition: G4PAIxSection.cc:2308
G4PAIxSection::SumOverIntervaldEdx
G4double SumOverIntervaldEdx(G4int intervalNumber)
Definition: G4PAIxSection.cc:1685
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
keV
static constexpr double keV
Definition: G4SIunits.hh:216
G4PAIxSection::GetStepCerenkovLoss
G4double GetStepCerenkovLoss(G4double step)
Definition: G4PAIxSection.cc:2281
G4PAIxSection::SumOverBorderdEdx
G4double SumOverBorderdEdx(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1935
G4PAIxSection::DifPAIxSection
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1225
plottest35.c1
tuple c1
Definition: plottest35.py:14
G4PAIxSection::RutherfordIntegral
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
Definition: G4PAIxSection.cc:1055
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4PAIxSection::GetElectronRange
G4double GetElectronRange(G4double energy)
Definition: G4PAIxSection.cc:1129