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