Geant4  10.01.p03
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 83662 2014-09-08 10:00:20Z gcosmo $
28 // GEANT4 tag $Name: geant4-09-03-ref-06 $
29 //
30 //
31 // G4PAIxSection.cc -- class implementation file
32 //
33 // GEANT 4 class implementation file
34 //
35 // For information related to this code, please, contact
36 // the Geant4 Collaboration.
37 //
38 // R&D: Vladimir.Grichine@cern.ch
39 //
40 // History:
41 //
42 // 13.05.03 V. Grichine, bug fixed for maxEnergyTransfer > max interval energy
43 // 28.05.01 V.Ivanchenko minor changes to provide ANSI -wall compilation
44 // 17.05.01 V. Grichine, low energy extension down to 10*keV of proton
45 // 20.11.98 adapted to a new Material/SandiaTable interface, mma
46 // 11.06.97 V. Grichine, 1st version
47 //
48 
49 #include "G4PAIxSection.hh"
50 
51 #include "globals.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4ios.hh"
55 #include "G4Poisson.hh"
56 #include "G4Material.hh"
57 #include "G4MaterialCutsCouple.hh"
58 #include "G4SandiaTable.hh"
59 
60 using namespace std;
61 
62 /* ******************************************************************
63 
64 // Init array of Lorentz factors
65 
66 const G4double G4PAIxSection::fLorentzFactor[22] =
67 {
68  0.0 , 1.1 , 1.2 , 1.3 , 1.5 , 1.8 , 2.0 ,
69  2.5 , 3.0 , 4.0 , 7.0 , 10.0 , 20.0 , 40.0 ,
70  70.0 , 100.0 , 300.0 , 600.0 , 1000.0 , 3000.0 ,
71  10000.0 , 50000.0
72 };
73 
74 const G4int G4PAIxSection::
75 fRefGammaNumber = 29; // The number of gamma for creation of
76  // spline (9)
77 
78 ***************************************************************** */
79 
80 // Local class constants
81 
82 const G4double G4PAIxSection::fDelta = 0.005; // 0.005 energy shift from interval border
83 const G4double G4PAIxSection::fError = 0.005; // 0.005 error in lin-log approximation
84 
85 const G4int G4PAIxSection::fMaxSplineSize = 1000; // Max size of output spline
86  // arrays
88 //
89 // Constructor
90 //
91 
93 {
94  fSandia = 0;
95  fMatSandiaMatrix = 0;
96  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
97  fIntervalNumber = fSplineNumber = 0;
98  fVerbose = 0;
99 
100  fSplineEnergy = G4DataVector(fMaxSplineSize,0.0);
101  fRePartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
102  fImPartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
103  fIntegralTerm = G4DataVector(fMaxSplineSize,0.0);
104  fDifPAIxSection = G4DataVector(fMaxSplineSize,0.0);
105  fdNdxCerenkov = G4DataVector(fMaxSplineSize,0.0);
106  fdNdxPlasmon = G4DataVector(fMaxSplineSize,0.0);
107  fdNdxMM = G4DataVector(fMaxSplineSize,0.0);
108  fdNdxResonance = G4DataVector(fMaxSplineSize,0.0);
109  fIntegralPAIxSection = G4DataVector(fMaxSplineSize,0.0);
110  fIntegralPAIdEdx = G4DataVector(fMaxSplineSize,0.0);
111  fIntegralCerenkov = G4DataVector(fMaxSplineSize,0.0);
112  fIntegralPlasmon = G4DataVector(fMaxSplineSize,0.0);
113  fIntegralMM = G4DataVector(fMaxSplineSize,0.0);
114  fIntegralResonance = G4DataVector(fMaxSplineSize,0.0);
115 
116  fMaterialIndex = 0;
117 
118  for( G4int i = 0; i < 500; ++i )
119  {
120  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
121  }
122 }
123 
125 //
126 // Constructor
127 //
128 
130 {
131  fDensity = matCC->GetMaterial()->GetDensity();
132  G4int matIndex = matCC->GetMaterial()->GetIndex();
133  fMaterialIndex = matIndex;
134 
135  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
136  fSandia = (*theMaterialTable)[matIndex]->GetSandiaTable();
137 
138  fVerbose = 0;
139 
140  G4int i, j;
141  fMatSandiaMatrix = new G4OrderedTable();
142 
143  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
144  {
145  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
146  }
147  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
148  {
149  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
150 
151  for(j = 1; j < 5; j++)
152  {
153  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
154  }
155  }
156  ComputeLowEnergyCof();
157  // fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0;
158 }
159 
161 
163  G4double maxEnergyTransfer)
164 {
165  fSandia = 0;
166  fMatSandiaMatrix = 0;
167  fVerbose = 0;
168  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
169  G4int i, j;
170 
171  fMaterialIndex = materialIndex;
172  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
173  fElectronDensity = (*theMaterialTable)[materialIndex]->
174  GetElectronDensity();
175  fIntervalNumber = (*theMaterialTable)[materialIndex]->
176  GetSandiaTable()->GetMatNbOfIntervals();
177  fIntervalNumber--;
178  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
179 
180  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
181  fA1 = G4DataVector(fIntervalNumber+2,0.0);
182  fA2 = G4DataVector(fIntervalNumber+2,0.0);
183  fA3 = G4DataVector(fIntervalNumber+2,0.0);
184  fA4 = G4DataVector(fIntervalNumber+2,0.0);
185 
186  for(i = 1; i <= fIntervalNumber; i++ )
187  {
188  if(((*theMaterialTable)[materialIndex]->
189  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
190  i > fIntervalNumber )
191  {
192  fEnergyInterval[i] = maxEnergyTransfer;
193  fIntervalNumber = i;
194  break;
195  }
196  fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
197  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
198  fA1[i] = (*theMaterialTable)[materialIndex]->
199  GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
200  fA2[i] = (*theMaterialTable)[materialIndex]->
201  GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
202  fA3[i] = (*theMaterialTable)[materialIndex]->
203  GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
204  fA4[i] = (*theMaterialTable)[materialIndex]->
205  GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
206  // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
207  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
208  }
209  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
210  {
211  fIntervalNumber++;
212  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
213  }
214 
215  // Now checking, if two borders are too close together
216 
217  for(i=1;i<fIntervalNumber;i++)
218  {
219  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
220  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
221  {
222  continue;
223  }
224  else
225  {
226  for(j=i;j<fIntervalNumber;j++)
227  {
228  fEnergyInterval[j] = fEnergyInterval[j+1];
229  fA1[j] = fA1[j+1];
230  fA2[j] = fA2[j+1];
231  fA3[j] = fA3[j+1];
232  fA4[j] = fA4[j+1];
233  }
234  fIntervalNumber--;
235  i--;
236  }
237  }
238 
239 
240  /* *********************************
241 
242  fSplineEnergy = new G4double[fMaxSplineSize];
243  fRePartDielectricConst = new G4double[fMaxSplineSize];
244  fImPartDielectricConst = new G4double[fMaxSplineSize];
245  fIntegralTerm = new G4double[fMaxSplineSize];
246  fDifPAIxSection = new G4double[fMaxSplineSize];
247  fIntegralPAIxSection = new G4double[fMaxSplineSize];
248 
249  for(i=0;i<fMaxSplineSize;i++)
250  {
251  fSplineEnergy[i] = 0.0;
252  fRePartDielectricConst[i] = 0.0;
253  fImPartDielectricConst[i] = 0.0;
254  fIntegralTerm[i] = 0.0;
255  fDifPAIxSection[i] = 0.0;
256  fIntegralPAIxSection[i] = 0.0;
257  }
258  ************************************************** */
259  ComputeLowEnergyCof();
260  InitPAI(); // create arrays allocated above
261  /*
262  delete[] fEnergyInterval;
263  delete[] fA1;
264  delete[] fA2;
265  delete[] fA3;
266  delete[] fA4;
267  */
268 }
269 
271 //
272 // Constructor called from G4PAIPhotonModel !!!
273 
275  G4double maxEnergyTransfer,
276  G4double betaGammaSq,
277  G4double** photoAbsCof,
278  G4int intNumber )
279 {
280  fSandia = 0;
281  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
282  fIntervalNumber = fSplineNumber = 0;
283  fVerbose = 0;
284 
285  fSplineEnergy = G4DataVector(500,0.0);
286  fRePartDielectricConst = G4DataVector(500,0.0);
287  fImPartDielectricConst = G4DataVector(500,0.0);
288  fIntegralTerm = G4DataVector(500,0.0);
289  fDifPAIxSection = G4DataVector(500,0.0);
290  fdNdxCerenkov = G4DataVector(500,0.0);
291  fdNdxPlasmon = G4DataVector(500,0.0);
292  fdNdxMM = G4DataVector(500,0.0);
293  fdNdxResonance = G4DataVector(500,0.0);
294  fIntegralPAIxSection = G4DataVector(500,0.0);
295  fIntegralPAIdEdx = G4DataVector(500,0.0);
296  fIntegralCerenkov = G4DataVector(500,0.0);
297  fIntegralPlasmon = G4DataVector(500,0.0);
298  fIntegralMM = G4DataVector(500,0.0);
299  fIntegralResonance = G4DataVector(500,0.0);
300 
301  for( G4int i = 0; i < 500; ++i )
302  {
303  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
304  }
305 
306  fSandia = 0;
307  fMatSandiaMatrix = 0;
308  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
309  G4int i, j;
310 
311  fMaterialIndex = materialIndex;
312  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
313  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
314 
315  fIntervalNumber = intNumber;
316  fIntervalNumber--;
317  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
318 
319  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
320  fA1 = G4DataVector(fIntervalNumber+2,0.0);
321  fA2 = G4DataVector(fIntervalNumber+2,0.0);
322  fA3 = G4DataVector(fIntervalNumber+2,0.0);
323  fA4 = G4DataVector(fIntervalNumber+2,0.0);
324 
325 
326  /*
327  fEnergyInterval = new G4double[fIntervalNumber+2];
328  fA1 = new G4double[fIntervalNumber+2];
329  fA2 = new G4double[fIntervalNumber+2];
330  fA3 = new G4double[fIntervalNumber+2];
331  fA4 = new G4double[fIntervalNumber+2];
332  */
333  for( i = 1; i <= fIntervalNumber; i++ )
334  {
335  if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
336  i > fIntervalNumber )
337  {
338  fEnergyInterval[i] = maxEnergyTransfer;
339  fIntervalNumber = i;
340  break;
341  }
342  fEnergyInterval[i] = photoAbsCof[i-1][0];
343  fA1[i] = photoAbsCof[i-1][1];
344  fA2[i] = photoAbsCof[i-1][2];
345  fA3[i] = photoAbsCof[i-1][3];
346  fA4[i] = photoAbsCof[i-1][4];
347  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
348  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
349  }
350  // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
351 
352  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
353  {
354  fIntervalNumber++;
355  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
356  }
357  // G4cout<<"after check of max energy transfer"<<G4endl;
358 
359  for( i = 1; i <= fIntervalNumber; i++ )
360  {
361  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
362  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
363  }
364  // Now checking, if two borders are too close together
365 
366  for( i = 1; i < fIntervalNumber; i++ )
367  {
368  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
369  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
370  {
371  continue;
372  }
373  else
374  {
375  for(j=i;j<fIntervalNumber;j++)
376  {
377  fEnergyInterval[j] = fEnergyInterval[j+1];
378  fA1[j] = fA1[j+1];
379  fA2[j] = fA2[j+1];
380  fA3[j] = fA3[j+1];
381  fA4[j] = fA4[j+1];
382  }
383  fIntervalNumber--;
384  i--;
385  }
386  }
387  // G4cout<<"after check of close borders"<<G4endl;
388 
389  for( i = 1; i <= fIntervalNumber; i++ )
390  {
391  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
392  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
393  }
394 
395  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
396 
397  ComputeLowEnergyCof();
398  G4double betaGammaSqRef =
399  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
400 
401  NormShift(betaGammaSqRef);
402  SplainPAI(betaGammaSqRef);
403 
404  // Preparation of integral PAI cross section for input betaGammaSq
405 
406  for(i = 1; i <= fSplineNumber; i++)
407  {
408  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
409  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
410  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
411  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
412  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
413 
414  // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
415  // <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
416  }
417  IntegralCerenkov();
418  IntegralMM();
419  IntegralPlasmon();
420  IntegralResonance();
421  IntegralPAIxSection();
422  /*
423  delete[] fEnergyInterval;
424  delete[] fA1;
425  delete[] fA2;
426  delete[] fA3;
427  delete[] fA4;
428  */
429 }
430 
432 //
433 // Test Constructor with beta*gamma square value
434 
436  G4double maxEnergyTransfer,
437  G4double betaGammaSq )
438 {
439  fSandia = 0;
440  fMatSandiaMatrix = 0;
441  fVerbose = 0;
442  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
443 
444  G4int i, j, numberOfElements;
445 
446  fMaterialIndex = materialIndex;
447  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
448  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
449  numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements();
450 
451  G4int* thisMaterialZ = new G4int[numberOfElements];
452 
453  for( i = 0; i < numberOfElements; i++ )
454  {
455  thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
456  GetElement(i)->GetZ();
457  }
458  // fSandia = new G4SandiaTable(materialIndex);
459  fSandia = (*theMaterialTable)[materialIndex]->GetSandiaTable();
460  G4SandiaTable thisMaterialSandiaTable(materialIndex);
461  fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals(thisMaterialZ,
462  numberOfElements);
463  fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
464  ( thisMaterialZ ,
465  (*theMaterialTable)[materialIndex]->GetFractionVector() ,
466  numberOfElements,fIntervalNumber);
467 
468  fIntervalNumber--;
469 
470  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
471  fA1 = G4DataVector(fIntervalNumber+2,0.0);
472  fA2 = G4DataVector(fIntervalNumber+2,0.0);
473  fA3 = G4DataVector(fIntervalNumber+2,0.0);
474  fA4 = G4DataVector(fIntervalNumber+2,0.0);
475 
476  /*
477  fEnergyInterval = new G4double[fIntervalNumber+2];
478  fA1 = new G4double[fIntervalNumber+2];
479  fA2 = new G4double[fIntervalNumber+2];
480  fA3 = new G4double[fIntervalNumber+2];
481  fA4 = new G4double[fIntervalNumber+2];
482  */
483  for( i = 1; i <= fIntervalNumber; i++ )
484  {
485  if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
486  i > fIntervalNumber)
487  {
488  fEnergyInterval[i] = maxEnergyTransfer;
489  fIntervalNumber = i;
490  break;
491  }
492  fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0);
493  fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity;
494  fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity;
495  fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity;
496  fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity;
497 
498  }
499  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
500  {
501  fIntervalNumber++;
502  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
503  fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
504  fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
505  fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
506  fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
507  }
508  for(i=1;i<=fIntervalNumber;i++)
509  {
510  // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
511  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
512  }
513  // Now checking, if two borders are too close together
514 
515  for( i = 1; i < fIntervalNumber; i++ )
516  {
517  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
518  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
519  {
520  continue;
521  }
522  else
523  {
524  for( j = i; j < fIntervalNumber; j++ )
525  {
526  fEnergyInterval[j] = fEnergyInterval[j+1];
527  fA1[j] = fA1[j+1];
528  fA2[j] = fA2[j+1];
529  fA3[j] = fA3[j+1];
530  fA4[j] = fA4[j+1];
531  }
532  fIntervalNumber--;
533  i--;
534  }
535  }
536 
537  /* *********************************
538  fSplineEnergy = new G4double[fMaxSplineSize];
539  fRePartDielectricConst = new G4double[fMaxSplineSize];
540  fImPartDielectricConst = new G4double[fMaxSplineSize];
541  fIntegralTerm = new G4double[fMaxSplineSize];
542  fDifPAIxSection = new G4double[fMaxSplineSize];
543  fIntegralPAIxSection = new G4double[fMaxSplineSize];
544 
545  for(i=0;i<fMaxSplineSize;i++)
546  {
547  fSplineEnergy[i] = 0.0;
548  fRePartDielectricConst[i] = 0.0;
549  fImPartDielectricConst[i] = 0.0;
550  fIntegralTerm[i] = 0.0;
551  fDifPAIxSection[i] = 0.0;
552  fIntegralPAIxSection[i] = 0.0;
553  }
554  */
555 
556  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
557 
558  ComputeLowEnergyCof();
559  G4double betaGammaSqRef =
560  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
561 
562  NormShift(betaGammaSqRef);
563  SplainPAI(betaGammaSqRef);
564 
565  // Preparation of integral PAI cross section for input betaGammaSq
566 
567  for(i = 1; i <= fSplineNumber; i++)
568  {
569  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
570  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
571  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
572  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
573  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
574  }
575  IntegralPAIxSection();
576  IntegralCerenkov();
577  IntegralMM();
578  IntegralPlasmon();
579  IntegralResonance();
580 }
581 
583 //
584 // Destructor
585 
587 {
588  /* ************************
589  delete[] fSplineEnergy ;
590  delete[] fRePartDielectricConst;
591  delete[] fImPartDielectricConst;
592  delete[] fIntegralTerm ;
593  delete[] fDifPAIxSection ;
594  delete[] fIntegralPAIxSection ;
595  */
596  delete fMatSandiaMatrix;
597 }
598 
599 
600 
601 
603 //
604 // Constructor with beta*gamma square value called from G4PAIPhotModel/Data
605 
606 void G4PAIxSection::Initialize( const G4Material* material,
607  G4double maxEnergyTransfer,
608  G4double betaGammaSq,
609  G4SandiaTable* sandia)
610 {
611  if(fVerbose > 0)
612  {
613  G4cout<<G4endl;
614  G4cout<<"G4PAIxSection::Initialize(...,G4SandiaTable* sandia)"<<G4endl;
615  G4cout<<G4endl;
616  }
617  G4int i, j;
618 
619  fSandia = sandia;
620  fIntervalNumber = sandia->GetMaxInterval();
621  fDensity = material->GetDensity();
622  fElectronDensity = material->GetElectronDensity();
623 
624  // fIntervalNumber--;
625 
626  if( fVerbose > 0 )
627  {
628  G4cout<<"fDensity = "<<fDensity<<"\t"<<fElectronDensity<<"\t fIntervalNumber = "<<fIntervalNumber<<G4endl;
629  }
630  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
631  fA1 = G4DataVector(fIntervalNumber+2,0.0);
632  fA2 = G4DataVector(fIntervalNumber+2,0.0);
633  fA3 = G4DataVector(fIntervalNumber+2,0.0);
634  fA4 = G4DataVector(fIntervalNumber+2,0.0);
635 
636  for( i = 1; i <= fIntervalNumber; i++ )
637  {
638  if ( sandia->GetSandiaMatTablePAI(i-1,0) < 1.*eV && sandia->GetLowerI1() == false)
639  {
640  fIntervalNumber--;
641  continue;
642  }
643  if( ( sandia->GetSandiaMatTablePAI(i-1,0) >= maxEnergyTransfer ) || i >= fIntervalNumber )
644  {
645  fEnergyInterval[i] = maxEnergyTransfer;
646  fIntervalNumber = i;
647  break;
648  }
649  fEnergyInterval[i] = sandia->GetSandiaMatTablePAI(i-1,0);
650  fA1[i] = sandia->GetSandiaMatTablePAI(i-1,1);
651  fA2[i] = sandia->GetSandiaMatTablePAI(i-1,2);
652  fA3[i] = sandia->GetSandiaMatTablePAI(i-1,3);
653  fA4[i] = sandia->GetSandiaMatTablePAI(i-1,4);
654 
655  if( fVerbose > 0 )
656  {
657  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
658  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
659  }
660  }
661  if( fVerbose > 0 ) G4cout<<"last i = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
662 
663  if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer )
664  {
665  fIntervalNumber++;
666  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
667  }
668  if( fVerbose > 0 )
669  {
670  for( i = 1; i <= fIntervalNumber; i++ )
671  {
672  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
673  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
674  }
675  }
676  if( fVerbose > 0 ) G4cout<<"Now checking, if two borders are too close together"<<G4endl;
677 
678  for( i = 1; i < fIntervalNumber; i++ )
679  {
680  if( fEnergyInterval[i+1]-fEnergyInterval[i] >
681  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]) && fEnergyInterval[i] > 0.) continue;
682  else
683  {
684  if( fVerbose > 0 ) G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fEnergyInterval[i+1]/keV;
685 
686  for( j = i; j < fIntervalNumber; j++ )
687  {
688  fEnergyInterval[j] = fEnergyInterval[j+1];
689  fA1[j] = fA1[j+1];
690  fA2[j] = fA2[j+1];
691  fA3[j] = fA3[j+1];
692  fA4[j] = fA4[j+1];
693  }
694  fIntervalNumber--;
695  i--;
696  }
697  }
698  if( fVerbose > 0 )
699  {
700  for( i = 1; i <= fIntervalNumber; i++ )
701  {
702  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
703  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
704  }
705  }
706  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
707 
708  ComputeLowEnergyCof(material);
709 
710  G4double betaGammaSqRef =
711  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
712 
713  NormShift(betaGammaSqRef);
714  SplainPAI(betaGammaSqRef);
715 
716  // Preparation of integral PAI cross section for input betaGammaSq
717 
718  for( i = 1; i <= fSplineNumber; i++ )
719  {
720  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
721 
722 
723  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
724  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
725  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
726  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
727  }
728  IntegralPAIxSection();
729  IntegralCerenkov();
730  IntegralMM();
731  IntegralPlasmon();
732  IntegralResonance();
733 
734  for( i = 1; i <= fSplineNumber; i++ )
735  {
736  if(fVerbose>0) G4cout<<i<<"; w = "<<fSplineEnergy[i]/keV<<" keV; dN/dx_>w = "<<fIntegralPAIxSection[i]<<" 1/mm"<<G4endl;
737  }
738 }
739 
740 
742 //
743 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
744 //
745 
747 {
748  G4int i, numberOfElements = material->GetNumberOfElements();
749  G4double sumZ = 0., sumCof = 0.;
750 
751  static const G4double p0 = 1.20923e+00;
752  static const G4double p1 = 3.53256e-01;
753  static const G4double p2 = -1.45052e-03;
754 
755  G4double* thisMaterialZ = new G4double[numberOfElements];
756  G4double* thisMaterialCof = new G4double[numberOfElements];
757 
758  for( i = 0; i < numberOfElements; i++ )
759  {
760  thisMaterialZ[i] = material->GetElement(i)->GetZ();
761  sumZ += thisMaterialZ[i];
762  thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
763  }
764  for( i = 0; i < numberOfElements; i++ )
765  {
766  sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
767  }
768  fLowEnergyCof = sumCof;
769  delete [] thisMaterialZ;
770  delete [] thisMaterialCof;
771  // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
772 }
773 
775 //
776 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
777 //
778 
780 {
781  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
782  G4int i, numberOfElements = (*theMaterialTable)[fMaterialIndex]->GetNumberOfElements();
783  G4double sumZ = 0., sumCof = 0.;
784 
785  const G4double p0 = 1.20923e+00;
786  const G4double p1 = 3.53256e-01;
787  const G4double p2 = -1.45052e-03;
788 
789  G4double* thisMaterialZ = new G4double[numberOfElements];
790  G4double* thisMaterialCof = new G4double[numberOfElements];
791 
792  for( i = 0; i < numberOfElements; i++ )
793  {
794  thisMaterialZ[i] = (*theMaterialTable)[fMaterialIndex]->GetElement(i)->GetZ();
795  sumZ += thisMaterialZ[i];
796  thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
797  }
798  for( i = 0; i < numberOfElements; i++ )
799  {
800  sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
801  }
802  fLowEnergyCof = sumCof;
803  // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
804  delete [] thisMaterialZ;
805  delete [] thisMaterialCof;
806 }
807 
809 //
810 // General control function for class G4PAIxSection
811 //
812 
814 {
815  G4int i;
816  G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
817  fLorentzFactor[fRefGammaNumber] - 1;
818 
819  // Preparation of integral PAI cross section for reference gamma
820 
821  NormShift(betaGammaSq);
822  SplainPAI(betaGammaSq);
823 
824  IntegralPAIxSection();
825  IntegralCerenkov();
826  IntegralMM();
827  IntegralPlasmon();
828  IntegralResonance();
829 
830  for(i = 0; i<= fSplineNumber; i++)
831  {
832  fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i];
833  if(i != 0)
834  {
835  fPAItable[i][0] = fSplineEnergy[i];
836  }
837  }
838  fPAItable[0][0] = fSplineNumber;
839 
840  for(G4int j = 1; j < 112; j++) // for other gammas
841  {
842  if( j == fRefGammaNumber ) continue;
843 
844  betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
845 
846  for(i = 1; i <= fSplineNumber; i++)
847  {
848  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
849  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
850  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
851  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
852  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
853  }
854  IntegralPAIxSection();
855  IntegralCerenkov();
856  IntegralMM();
857  IntegralPlasmon();
858  IntegralResonance();
859 
860  for(i = 0; i <= fSplineNumber; i++)
861  {
862  fPAItable[i][j] = fIntegralPAIxSection[i];
863  }
864  }
865 
866 }
867 
869 //
870 // Shifting from borders to intervals Creation of first energy points
871 //
872 
874 {
875  G4int i, j;
876 
877  if(fVerbose>0) G4cout<<" G4PAIxSection::NormShift call "<<G4endl;
878 
879 
880  for( i = 1; i <= fIntervalNumber-1; i++ )
881  {
882  for( j = 1; j <= 2; j++ )
883  {
884  fSplineNumber = (i-1)*2 + j;
885 
886  if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
887  else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
888  if(fVerbose>0) G4cout<<"cn = "<<fSplineNumber<<"; "<<"w = "<<fSplineEnergy[fSplineNumber]/keV<<" keV"<<G4endl;
889  }
890  }
891  fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
892 
893  j = 1;
894 
895  for( i = 2; i <= fSplineNumber; i++ )
896  {
897  if( fSplineEnergy[i]<fEnergyInterval[j+1] )
898  {
899  fIntegralTerm[i] = fIntegralTerm[i-1] +
900  RutherfordIntegral(j,fSplineEnergy[i-1],
901  fSplineEnergy[i] );
902  }
903  else
904  {
905  G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
906  fEnergyInterval[j+1] );
907  j++;
908  fIntegralTerm[i] = fIntegralTerm[i-1] + x +
909  RutherfordIntegral(j,fEnergyInterval[j],
910  fSplineEnergy[i] );
911  }
912  if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV \t"<<fIntegralTerm[i]<<"\n"<<G4endl;
913  }
914  fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
915  fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
916 
917  // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl;
918 
919  // Calculation of PAI differrential cross-section (1/(keV*cm))
920  // in the energy points near borders of energy intervals
921 
922  for(G4int k = 1; k <= fIntervalNumber-1; k++ )
923  {
924  for( j = 1; j <= 2; j++ )
925  {
926  i = (k-1)*2 + j;
927  fImPartDielectricConst[i] = fNormalizationCof*
928  ImPartDielectricConst(k,fSplineEnergy[i]);
929  fRePartDielectricConst[i] = fNormalizationCof*
930  RePartDielectricConst(fSplineEnergy[i]);
931  fIntegralTerm[i] *= fNormalizationCof;
932 
933  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
934  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
935  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
936  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
937  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
938  if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV, xsc = "<<fDifPAIxSection[i]<<"\n"<<G4endl;
939  }
940  }
941 
942 } // end of NormShift
943 
945 //
946 // Creation of new energy points as geometrical mean of existing
947 // one, calculation PAI_cs for them, while the error of logarithmic
948 // linear approximation would be smaller than 'fError'
949 
951 {
952  G4int j, k = 1, i = 1;
953 
954  if(fVerbose>0) G4cout<<" G4PAIxSection::SplainPAI call "<<G4endl;
955 
956  while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
957  {
958  // if( std::abs(fSplineEnergy[i+1]-fEnergyInterval[k+1]) > (fSplineEnergy[i+1]+fEnergyInterval[k+1])*5.e-7 )
959  if( fSplineEnergy[i+1] > fEnergyInterval[k+1] )
960  {
961  k++; // Here next energy point is in next energy interval
962  i++;
963  if(fVerbose>0) G4cout<<" in if: i = "<<i<<"; k = "<<k<<G4endl;
964  continue;
965  }
966  if(fVerbose>0) G4cout<<" out if: i = "<<i<<"; k = "<<k<<G4endl;
967 
968  // Shifting of arrayes for inserting the geometrical
969  // average of 'i' and 'i+1' energy points to 'i+1' place
970  fSplineNumber++;
971 
972  for( j = fSplineNumber; j >= i+2; j-- )
973  {
974  fSplineEnergy[j] = fSplineEnergy[j-1];
975  fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
976  fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
977  fIntegralTerm[j] = fIntegralTerm[j-1];
978 
979  fDifPAIxSection[j] = fDifPAIxSection[j-1];
980  fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
981  fdNdxMM[j] = fdNdxMM[j-1];
982  fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
983  fdNdxResonance[j] = fdNdxResonance[j-1];
984  }
985  G4double x1 = fSplineEnergy[i];
986  G4double x2 = fSplineEnergy[i+1];
987  G4double yy1 = fDifPAIxSection[i];
988  G4double y2 = fDifPAIxSection[i+1];
989 
990  if(fVerbose>0) G4cout<<"Spline: x1 = "<<x1<<"; x2 = "<<x2<<", yy1 = "<<yy1<<"; y2 = "<<y2<<G4endl;
991 
992 
993  G4double en1 = sqrt(x1*x2);
994  // G4double en1 = 0.5*(x1 + x2);
995 
996 
997  fSplineEnergy[i+1] = en1;
998 
999  // Calculation of logarithmic linear approximation
1000  // in this (enr) energy point, which number is 'i+1' now
1001 
1002  G4double a = log10(y2/yy1)/log10(x2/x1);
1003  G4double b = log10(yy1) - a*log10(x1);
1004  G4double y = a*log10(en1) + b;
1005 
1006  y = pow(10.,y);
1007 
1008  // Calculation of the PAI dif. cross-section at this point
1009 
1010  fImPartDielectricConst[i+1] = fNormalizationCof*
1011  ImPartDielectricConst(k,fSplineEnergy[i+1]);
1012  fRePartDielectricConst[i+1] = fNormalizationCof*
1013  RePartDielectricConst(fSplineEnergy[i+1]);
1014  fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
1015  RutherfordIntegral(k,fSplineEnergy[i],
1016  fSplineEnergy[i+1]);
1017 
1018  fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
1019  fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
1020  fdNdxMM[i+1] = PAIdNdxMM(i+1,betaGammaSq);
1021  fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
1022  fdNdxResonance[i+1] = PAIdNdxResonance(i+1,betaGammaSq);
1023 
1024  // Condition for next division of this segment or to pass
1025 
1026  if(fVerbose>0) G4cout<<"Spline, a = "<<a<<"; b = "<<b<<"; new xsc = "<<y<<"; compxsc = "<<fDifPAIxSection[i+1]<<G4endl;
1027 
1028  // to higher energies
1029 
1030  G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
1031 
1032  G4double delta = 2.*(fSplineEnergy[i+1]-fSplineEnergy[i])/(fSplineEnergy[i+1]+fSplineEnergy[i]);
1033 
1034  if( x < 0 )
1035  {
1036  x = -x;
1037  }
1038  if( x > fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*fDelta )
1039  {
1040  continue; // next division
1041  }
1042  i += 2; // pass to next segment
1043 
1044  } // 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 
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 
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 
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 
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 
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 
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 
1302  G4double betaGammaSq )
1303 {
1304  G4double logarithm, x3, x5, argument, modul2, dNdxC;
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 
1360  G4double betaGammaSq )
1361 {
1362  G4double logarithm, x3, x5, argument, dNdxC;
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 
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 
1450  G4double betaGammaSq )
1451 {
1452  G4double resonance, modul2, dNdxP;
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 
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 
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 
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 
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 
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 
1638 {
1639  G4double x0,x1,y0,yy1,a,b,c,result;
1640 
1641  x0 = fSplineEnergy[i];
1642  x1 = fSplineEnergy[i+1];
1643  if(fVerbose>0) G4cout<<"SumOverInterval i= " << i << " x0 = "<<x0<<"; x1 = "<<x1<<G4endl;
1644 
1645  if( x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1646 
1647  y0 = fDifPAIxSection[i];
1648  yy1 = fDifPAIxSection[i+1];
1649 
1650  if(fVerbose>0) G4cout<<"x0 = "<<x0<<"; x1 = "<<x1<<", y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1651 
1652  c = x1/x0;
1653  a = log10(yy1/y0)/log10(c);
1654 
1655  if(fVerbose>0) G4cout<<"SumOverInterval, a = "<<a<<"; c = "<<c<<G4endl;
1656 
1657  // b = log10(y0) - a*log10(x0);
1658  b = y0/pow(x0,a);
1659  a += 1.;
1660  if( std::fabs(a) < 1.e-6 )
1661  {
1662  result = b*log(x1/x0);
1663  }
1664  else
1665  {
1666  result = y0*(x1*pow(c,a-1) - x0)/a;
1667  }
1668  a += 1.;
1669  if( std::fabs(a) < 1.e-6 )
1670  {
1671  fIntegralPAIxSection[0] += b*log(x1/x0);
1672  }
1673  else
1674  {
1675  fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1676  }
1677  if(fVerbose>0) G4cout<<"SumOverInterval, result = "<<result<<G4endl;
1678  return result;
1679 
1680 } // end of SumOverInterval
1681 
1683 
1685 {
1686  G4double x0,x1,y0,yy1,a,b,c,result;
1687 
1688  x0 = fSplineEnergy[i];
1689  x1 = fSplineEnergy[i+1];
1690 
1691  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1692 
1693  y0 = fDifPAIxSection[i];
1694  yy1 = fDifPAIxSection[i+1];
1695  c = x1/x0;
1696  a = log10(yy1/y0)/log10(c);
1697  // b = log10(y0) - a*log10(x0);
1698  b = y0/pow(x0,a);
1699  a += 2;
1700  if(a == 0)
1701  {
1702  result = b*log(x1/x0);
1703  }
1704  else
1705  {
1706  result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1707  }
1708  return result;
1709 
1710 } // end of SumOverInterval
1711 
1713 //
1714 // Calculation the PAI Cerenkov integral cross-section inside
1715 // of interval of continuous values of photo-ionisation Cerenkov
1716 // cross-section. Parameter 'i' is the number of interval.
1717 
1719 {
1720  G4double x0,x1,y0,yy1,a,b,c,result;
1721 
1722  x0 = fSplineEnergy[i];
1723  x1 = fSplineEnergy[i+1];
1724 
1725  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1726 
1727  y0 = fdNdxCerenkov[i];
1728  yy1 = fdNdxCerenkov[i+1];
1729  // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1730  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1731 
1732  c = x1/x0;
1733  a = log10(yy1/y0)/log10(c);
1734  b = y0/pow(x0,a);
1735 
1736  a += 1.0;
1737  if(a == 0) result = b*log(c);
1738  else result = y0*(x1*pow(c,a-1) - x0)/a;
1739  a += 1.0;
1740 
1741  if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
1742  else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1743  // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1744  return result;
1745 
1746 } // end of SumOverInterCerenkov
1747 
1749 //
1750 // Calculation the PAI MM-Cerenkov integral cross-section inside
1751 // of interval of continuous values of photo-ionisation Cerenkov
1752 // cross-section. Parameter 'i' is the number of interval.
1753 
1755 {
1756  G4double x0,x1,y0,yy1,a,b,c,result;
1757 
1758  x0 = fSplineEnergy[i];
1759  x1 = fSplineEnergy[i+1];
1760 
1761  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1762 
1763  y0 = fdNdxMM[i];
1764  yy1 = fdNdxMM[i+1];
1765  //G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1766  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1767 
1768  c = x1/x0;
1769  //G4cout<<" c = "<<c<< " yy1/y0= " << yy1/y0 <<G4endl;
1770  a = log10(yy1/y0)/log10(c);
1771  if(a > 10.0) return 0.;
1772  b = y0/pow(x0,a);
1773 
1774  a += 1.0;
1775  if(a == 0) result = b*log(c);
1776  else result = y0*(x1*pow(c,a-1) - x0)/a;
1777  a += 1.0;
1778 
1779  if( a == 0 ) fIntegralMM[0] += b*log(c);
1780  else fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1781  //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1782  return result;
1783 
1784 } // end of SumOverInterMM
1785 
1787 //
1788 // Calculation the PAI Plasmon integral cross-section inside
1789 // of interval of continuous values of photo-ionisation Plasmon
1790 // cross-section. Parameter 'i' is the number of interval.
1791 
1793 {
1794  G4double x0,x1,y0,yy1,a,b,c,result;
1795 
1796  x0 = fSplineEnergy[i];
1797  x1 = fSplineEnergy[i+1];
1798 
1799  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1800 
1801  y0 = fdNdxPlasmon[i];
1802  yy1 = fdNdxPlasmon[i+1];
1803  c =x1/x0;
1804  a = log10(yy1/y0)/log10(c);
1805  if(a > 10.0) return 0.;
1806  // b = log10(y0) - a*log10(x0);
1807  b = y0/pow(x0,a);
1808 
1809  a += 1.0;
1810  if(a == 0) result = b*log(x1/x0);
1811  else result = y0*(x1*pow(c,a-1) - x0)/a;
1812  a += 1.0;
1813 
1814  if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
1815  else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1816 
1817  return result;
1818 
1819 } // end of SumOverInterPlasmon
1820 
1822 //
1823 // Calculation the PAI resonance integral cross-section inside
1824 // of interval of continuous values of photo-ionisation resonance
1825 // cross-section. Parameter 'i' is the number of interval.
1826 
1828 {
1829  G4double x0,x1,y0,yy1,a,b,c,result;
1830 
1831  x0 = fSplineEnergy[i];
1832  x1 = fSplineEnergy[i+1];
1833 
1834  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1835 
1836  y0 = fdNdxResonance[i];
1837  yy1 = fdNdxResonance[i+1];
1838  c =x1/x0;
1839  a = log10(yy1/y0)/log10(c);
1840  if(a > 10.0) return 0.;
1841  // b = log10(y0) - a*log10(x0);
1842  b = y0/pow(x0,a);
1843 
1844  a += 1.0;
1845  if(a == 0) result = b*log(x1/x0);
1846  else result = y0*(x1*pow(c,a-1) - x0)/a;
1847  a += 1.0;
1848 
1849  if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0);
1850  else fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1851 
1852  return result;
1853 
1854 } // end of SumOverInterResonance
1855 
1857 //
1858 // Integration of PAI cross-section for the case of
1859 // passing across border between intervals
1860 
1862  G4double en0 )
1863 {
1864  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1865 
1866  e0 = en0;
1867  x0 = fSplineEnergy[i];
1868  x1 = fSplineEnergy[i+1];
1869  y0 = fDifPAIxSection[i];
1870  yy1 = fDifPAIxSection[i+1];
1871 
1872  //c = x1/x0;
1873  d = e0/x0;
1874  a = log10(yy1/y0)/log10(x1/x0);
1875  if(a > 10.0) return 0.;
1876 
1877  if(fVerbose>0) G4cout<<"SumOverBorder, a = "<<a<<G4endl;
1878 
1879  // b0 = log10(y0) - a*log10(x0);
1880  b = y0/pow(x0,a); // pow(10.,b);
1881 
1882  a += 1.;
1883  if( std::fabs(a) < 1.e-6 )
1884  {
1885  result = b*log(x0/e0);
1886  }
1887  else
1888  {
1889  result = y0*(x0 - e0*pow(d,a-1))/a;
1890  }
1891  a += 1.;
1892  if( std::fabs(a) < 1.e-6 )
1893  {
1894  fIntegralPAIxSection[0] += b*log(x0/e0);
1895  }
1896  else
1897  {
1898  fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1899  }
1900  x0 = fSplineEnergy[i - 1];
1901  x1 = fSplineEnergy[i - 2];
1902  y0 = fDifPAIxSection[i - 1];
1903  yy1 = fDifPAIxSection[i - 2];
1904 
1905  //c = x1/x0;
1906  d = e0/x0;
1907  a = log10(yy1/y0)/log10(x1/x0);
1908  // b0 = log10(y0) - a*log10(x0);
1909  b = y0/pow(x0,a);
1910  a += 1.;
1911  if( std::fabs(a) < 1.e-6 )
1912  {
1913  result += b*log(e0/x0);
1914  }
1915  else
1916  {
1917  result += y0*(e0*pow(d,a-1) - x0)/a;
1918  }
1919  a += 1.;
1920  if( std::fabs(a) < 1.e-6 )
1921  {
1922  fIntegralPAIxSection[0] += b*log(e0/x0);
1923  }
1924  else
1925  {
1926  fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1927  }
1928  return result;
1929 
1930 }
1931 
1933 
1935  G4double en0 )
1936 {
1937  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1938 
1939  e0 = en0;
1940  x0 = fSplineEnergy[i];
1941  x1 = fSplineEnergy[i+1];
1942  y0 = fDifPAIxSection[i];
1943  yy1 = fDifPAIxSection[i+1];
1944 
1945  //c = x1/x0;
1946  d = e0/x0;
1947  a = log10(yy1/y0)/log10(x1/x0);
1948  if(a > 10.0) return 0.;
1949  // b0 = log10(y0) - a*log10(x0);
1950  b = y0/pow(x0,a); // pow(10.,b);
1951 
1952  a += 2;
1953  if(a == 0)
1954  {
1955  result = b*log(x0/e0);
1956  }
1957  else
1958  {
1959  result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1960  }
1961  x0 = fSplineEnergy[i - 1];
1962  x1 = fSplineEnergy[i - 2];
1963  y0 = fDifPAIxSection[i - 1];
1964  yy1 = fDifPAIxSection[i - 2];
1965 
1966  // c = x1/x0;
1967  d = e0/x0;
1968  a = log10(yy1/y0)/log10(x1/x0);
1969  // b0 = log10(y0) - a*log10(x0);
1970  b = y0/pow(x0,a);
1971  a += 2;
1972  if(a == 0)
1973  {
1974  result += b*log(e0/x0);
1975  }
1976  else
1977  {
1978  result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1979  }
1980  return result;
1981 
1982 }
1983 
1985 //
1986 // Integration of Cerenkov cross-section for the case of
1987 // passing across border between intervals
1988 
1990  G4double en0 )
1991 {
1992  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1993 
1994  e0 = en0;
1995  x0 = fSplineEnergy[i];
1996  x1 = fSplineEnergy[i+1];
1997  y0 = fdNdxCerenkov[i];
1998  yy1 = fdNdxCerenkov[i+1];
1999 
2000  // G4cout<<G4endl;
2001  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2002  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2003  c = x1/x0;
2004  d = e0/x0;
2005  a = log10(yy1/y0)/log10(c);
2006  if(a > 10.0) return 0.;
2007  // b0 = log10(y0) - a*log10(x0);
2008  b = y0/pow(x0,a); // pow(10.,b0);
2009 
2010  a += 1.0;
2011  if( a == 0 ) result = b*log(x0/e0);
2012  else result = y0*(x0 - e0*pow(d,a-1))/a;
2013  a += 1.0;
2014 
2015  if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
2016  else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2017 
2018 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2019 
2020  x0 = fSplineEnergy[i - 1];
2021  x1 = fSplineEnergy[i - 2];
2022  y0 = fdNdxCerenkov[i - 1];
2023  yy1 = fdNdxCerenkov[i - 2];
2024 
2025  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2026  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2027 
2028  c = x1/x0;
2029  d = e0/x0;
2030  a = log10(yy1/y0)/log10(x1/x0);
2031  // b0 = log10(y0) - a*log10(x0);
2032  b = y0/pow(x0,a); // pow(10.,b0);
2033 
2034  a += 1.0;
2035  if( a == 0 ) result += b*log(e0/x0);
2036  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2037  a += 1.0;
2038 
2039  if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
2040  else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2041 
2042  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2043  // <<b<<"; result = "<<result<<G4endl;
2044 
2045  return result;
2046 
2047 }
2048 
2050 //
2051 // Integration of MM-Cerenkov cross-section for the case of
2052 // passing across border between intervals
2053 
2055  G4double en0 )
2056 {
2057  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
2058 
2059  e0 = en0;
2060  x0 = fSplineEnergy[i];
2061  x1 = fSplineEnergy[i+1];
2062  y0 = fdNdxMM[i];
2063  yy1 = fdNdxMM[i+1];
2064 
2065  // G4cout<<G4endl;
2066  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2067  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2068  c = x1/x0;
2069  d = e0/x0;
2070  a = log10(yy1/y0)/log10(c);
2071  if(a > 10.0) return 0.;
2072  // b0 = log10(y0) - a*log10(x0);
2073  b = y0/pow(x0,a); // pow(10.,b0);
2074 
2075  a += 1.0;
2076  if( a == 0 ) result = b*log(x0/e0);
2077  else result = y0*(x0 - e0*pow(d,a-1))/a;
2078  a += 1.0;
2079 
2080  if( a == 0 ) fIntegralMM[0] += b*log(x0/e0);
2081  else fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2082 
2083 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2084 
2085  x0 = fSplineEnergy[i - 1];
2086  x1 = fSplineEnergy[i - 2];
2087  y0 = fdNdxMM[i - 1];
2088  yy1 = fdNdxMM[i - 2];
2089 
2090  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2091  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2092 
2093  c = x1/x0;
2094  d = e0/x0;
2095  a = log10(yy1/y0)/log10(x1/x0);
2096  // b0 = log10(y0) - a*log10(x0);
2097  b = y0/pow(x0,a); // pow(10.,b0);
2098 
2099  a += 1.0;
2100  if( a == 0 ) result += b*log(e0/x0);
2101  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2102  a += 1.0;
2103 
2104  if( a == 0 ) fIntegralMM[0] += b*log(e0/x0);
2105  else fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2106 
2107  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2108  // <<b<<"; result = "<<result<<G4endl;
2109 
2110  return result;
2111 
2112 }
2113 
2115 //
2116 // Integration of Plasmon cross-section for the case of
2117 // passing across border between intervals
2118 
2120  G4double en0 )
2121 {
2122  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2123 
2124  e0 = en0;
2125  x0 = fSplineEnergy[i];
2126  x1 = fSplineEnergy[i+1];
2127  y0 = fdNdxPlasmon[i];
2128  yy1 = fdNdxPlasmon[i+1];
2129 
2130  c = x1/x0;
2131  d = e0/x0;
2132  a = log10(yy1/y0)/log10(c);
2133  if(a > 10.0) return 0.;
2134  // b0 = log10(y0) - a*log10(x0);
2135  b = y0/pow(x0,a); //pow(10.,b);
2136 
2137  a += 1.0;
2138  if( a == 0 ) result = b*log(x0/e0);
2139  else result = y0*(x0 - e0*pow(d,a-1))/a;
2140  a += 1.0;
2141 
2142  if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
2143  else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2144 
2145  x0 = fSplineEnergy[i - 1];
2146  x1 = fSplineEnergy[i - 2];
2147  y0 = fdNdxPlasmon[i - 1];
2148  yy1 = fdNdxPlasmon[i - 2];
2149 
2150  c = x1/x0;
2151  d = e0/x0;
2152  a = log10(yy1/y0)/log10(c);
2153  // b0 = log10(y0) - a*log10(x0);
2154  b = y0/pow(x0,a);// pow(10.,b0);
2155 
2156  a += 1.0;
2157  if( a == 0 ) result += b*log(e0/x0);
2158  else result += y0*(e0*pow(d,a-1) - x0)/a;
2159  a += 1.0;
2160 
2161  if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
2162  else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2163 
2164  return result;
2165 
2166 }
2167 
2169 //
2170 // Integration of resonance cross-section for the case of
2171 // passing across border between intervals
2172 
2174  G4double en0 )
2175 {
2176  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2177 
2178  e0 = en0;
2179  x0 = fSplineEnergy[i];
2180  x1 = fSplineEnergy[i+1];
2181  y0 = fdNdxResonance[i];
2182  yy1 = fdNdxResonance[i+1];
2183 
2184  c = x1/x0;
2185  d = e0/x0;
2186  a = log10(yy1/y0)/log10(c);
2187  if(a > 10.0) return 0.;
2188  // b0 = log10(y0) - a*log10(x0);
2189  b = y0/pow(x0,a); //pow(10.,b);
2190 
2191  a += 1.0;
2192  if( a == 0 ) result = b*log(x0/e0);
2193  else result = y0*(x0 - e0*pow(d,a-1))/a;
2194  a += 1.0;
2195 
2196  if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0);
2197  else fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2198 
2199  x0 = fSplineEnergy[i - 1];
2200  x1 = fSplineEnergy[i - 2];
2201  y0 = fdNdxResonance[i - 1];
2202  yy1 = fdNdxResonance[i - 2];
2203 
2204  c = x1/x0;
2205  d = e0/x0;
2206  a = log10(yy1/y0)/log10(c);
2207  // b0 = log10(y0) - a*log10(x0);
2208  b = y0/pow(x0,a);// pow(10.,b0);
2209 
2210  a += 1.0;
2211  if( a == 0 ) result += b*log(e0/x0);
2212  else result += y0*(e0*pow(d,a-1) - x0)/a;
2213  a += 1.0;
2214 
2215  if( a == 0 ) fIntegralResonance[0] += b*log(e0/x0);
2216  else fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2217 
2218  return result;
2219 
2220 }
2221 
2223 //
2224 // Returns random PAI-total energy loss over step
2225 
2227 {
2228  G4long numOfCollisions;
2229  G4double meanNumber, loss = 0.0;
2230 
2231  // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl;
2232 
2233  meanNumber = fIntegralPAIxSection[1]*step;
2234  numOfCollisions = G4Poisson(meanNumber);
2235 
2236  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2237 
2238  while(numOfCollisions)
2239  {
2240  loss += GetEnergyTransfer();
2241  numOfCollisions--;
2242  }
2243  // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
2244 
2245  return loss;
2246 }
2247 
2249 //
2250 // Returns random PAI-total energy transfer in one collision
2251 
2253 {
2254  G4int iTransfer ;
2255 
2256  G4double energyTransfer, position;
2257 
2258  position = fIntegralPAIxSection[1]*G4UniformRand();
2259 
2260  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2261  {
2262  if( position >= fIntegralPAIxSection[iTransfer] ) break;
2263  }
2264  if(iTransfer > fSplineNumber) iTransfer--;
2265 
2266  energyTransfer = fSplineEnergy[iTransfer];
2267 
2268  if(iTransfer > 1)
2269  {
2270  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2271  }
2272  return energyTransfer;
2273 }
2274 
2276 //
2277 // Returns random Cerenkov energy loss over step
2278 
2280 {
2281  G4long numOfCollisions;
2282  G4double meanNumber, loss = 0.0;
2283 
2284  // G4cout<<" G4PAIxSection::GetStepCerenkovLoss "<<G4endl;
2285 
2286  meanNumber = fIntegralCerenkov[1]*step;
2287  numOfCollisions = G4Poisson(meanNumber);
2288 
2289  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2290 
2291  while(numOfCollisions)
2292  {
2293  loss += GetCerenkovEnergyTransfer();
2294  numOfCollisions--;
2295  }
2296  // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2297 
2298  return loss;
2299 }
2300 
2302 //
2303 // Returns random MM-Cerenkov energy loss over step
2304 
2306 {
2307  G4long numOfCollisions;
2308  G4double meanNumber, loss = 0.0;
2309 
2310  // G4cout<<" G4PAIxSection::GetStepMMLoss "<<G4endl;
2311 
2312  meanNumber = fIntegralMM[1]*step;
2313  numOfCollisions = G4Poisson(meanNumber);
2314 
2315  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2316 
2317  while(numOfCollisions)
2318  {
2319  loss += GetMMEnergyTransfer();
2320  numOfCollisions--;
2321  }
2322  // G4cout<<"PAI MM-Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2323 
2324  return loss;
2325 }
2326 
2328 //
2329 // Returns Cerenkov energy transfer in one collision
2330 
2332 {
2333  G4int iTransfer ;
2334 
2335  G4double energyTransfer, position;
2336 
2337  position = fIntegralCerenkov[1]*G4UniformRand();
2338 
2339  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2340  {
2341  if( position >= fIntegralCerenkov[iTransfer] ) break;
2342  }
2343  if(iTransfer > fSplineNumber) iTransfer--;
2344 
2345  energyTransfer = fSplineEnergy[iTransfer];
2346 
2347  if(iTransfer > 1)
2348  {
2349  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2350  }
2351  return energyTransfer;
2352 }
2353 
2355 //
2356 // Returns MM-Cerenkov energy transfer in one collision
2357 
2359 {
2360  G4int iTransfer ;
2361 
2362  G4double energyTransfer, position;
2363 
2364  position = fIntegralMM[1]*G4UniformRand();
2365 
2366  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2367  {
2368  if( position >= fIntegralMM[iTransfer] ) break;
2369  }
2370  if(iTransfer > fSplineNumber) iTransfer--;
2371 
2372  energyTransfer = fSplineEnergy[iTransfer];
2373 
2374  if(iTransfer > 1)
2375  {
2376  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2377  }
2378  return energyTransfer;
2379 }
2380 
2382 //
2383 // Returns random plasmon energy loss over step
2384 
2386 {
2387  G4long numOfCollisions;
2388  G4double meanNumber, loss = 0.0;
2389 
2390  // G4cout<<" G4PAIxSection::GetStepPlasmonLoss "<<G4endl;
2391 
2392  meanNumber = fIntegralPlasmon[1]*step;
2393  numOfCollisions = G4Poisson(meanNumber);
2394 
2395  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2396 
2397  while(numOfCollisions)
2398  {
2399  loss += GetPlasmonEnergyTransfer();
2400  numOfCollisions--;
2401  }
2402  // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
2403 
2404  return loss;
2405 }
2406 
2408 //
2409 // Returns plasmon energy transfer in one collision
2410 
2412 {
2413  G4int iTransfer ;
2414 
2415  G4double energyTransfer, position;
2416 
2417  position = fIntegralPlasmon[1]*G4UniformRand();
2418 
2419  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2420  {
2421  if( position >= fIntegralPlasmon[iTransfer] ) break;
2422  }
2423  if(iTransfer > fSplineNumber) iTransfer--;
2424 
2425  energyTransfer = fSplineEnergy[iTransfer];
2426 
2427  if(iTransfer > 1)
2428  {
2429  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2430  }
2431  return energyTransfer;
2432 }
2433 
2435 //
2436 // Returns random resonance energy loss over step
2437 
2439 {
2440  G4long numOfCollisions;
2441  G4double meanNumber, loss = 0.0;
2442 
2443  // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl;
2444 
2445  meanNumber = fIntegralResonance[1]*step;
2446  numOfCollisions = G4Poisson(meanNumber);
2447 
2448  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2449 
2450  while(numOfCollisions)
2451  {
2452  loss += GetResonanceEnergyTransfer();
2453  numOfCollisions--;
2454  }
2455  // G4cout<<"PAI resonance loss = "<<loss/keV<<" keV"<<G4endl;
2456 
2457  return loss;
2458 }
2459 
2460 
2462 //
2463 // Returns resonance energy transfer in one collision
2464 
2466 {
2467  G4int iTransfer ;
2468 
2469  G4double energyTransfer, position;
2470 
2471  position = fIntegralResonance[1]*G4UniformRand();
2472 
2473  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2474  {
2475  if( position >= fIntegralResonance[iTransfer] ) break;
2476  }
2477  if(iTransfer > fSplineNumber) iTransfer--;
2478 
2479  energyTransfer = fSplineEnergy[iTransfer];
2480 
2481  if(iTransfer > 1)
2482  {
2483  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2484  }
2485  return energyTransfer;
2486 }
2487 
2488 
2490 //
2491 // Returns Rutherford energy transfer in one collision
2492 
2494 {
2495  G4int iTransfer ;
2496 
2497  G4double energyTransfer, position;
2498 
2499  position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand();
2500 
2501  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2502  {
2503  if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break;
2504  }
2505  if(iTransfer > fSplineNumber) iTransfer--;
2506 
2507  energyTransfer = fSplineEnergy[iTransfer];
2508 
2509  if(iTransfer > 1)
2510  {
2511  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2512  }
2513  return energyTransfer;
2514 }
2515 
2517 //
2518 
2519 void G4PAIxSection::CallError(G4int i, const G4String& methodName) const
2520 {
2521  G4String head = "G4PAIxSection::" + methodName + "()";
2523  ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
2524  G4Exception(head,"pai001",FatalException,ed);
2525 }
2526 
2528 //
2529 // Init array of Lorentz factors
2530 //
2531 
2533 
2534 const G4double G4PAIxSection::fLorentzFactor[112] = // fNumberOfGammas+1
2535 {
2536 0.0,
2537 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
2538 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
2539 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
2540 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
2541 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
2542 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
2543 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
2544 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
2545 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
2546 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
2547 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
2548 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
2549 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
2550 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
2551 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
2552 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
2553 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
2554 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
2555 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
2556 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
2557 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
2558 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
2559 1.065799e+05
2560 };
2561 
2563 //
2564 // The number of gamma for creation of spline (near ion-min , G ~ 4 )
2565 //
2566 
2567 const
2569 
2570 
2571 //
2572 // end of G4PAIxSection implementation file
2573 //
2575 
void IntegralPAIxSection()
void SplainPAI(G4double betaGammaSq)
G4double GetRutherfordEnergyTransfer()
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
void IntegralCerenkov()
static c2_factory< G4double > c2
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4double GetPlasmonEnergyTransfer()
static const double cm2
Definition: G4SIunits.hh:107
void IntegralPlasmon()
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void IntegralResonance()
G4double GetStepResonanceLoss(G4double step)
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
G4double SumOverBordPlasmon(G4int intervalNumber, G4double energy)
size_t GetIndex() const
Definition: G4Material.hh:262
G4double SumOverInterMM(G4int intervalNumber)
static const G4double betaBohr
const G4double pi
G4double GetZ() const
Definition: G4Element.hh:131
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:588
std::vector< G4Material * > G4MaterialTable
static const G4double cofBetaBohr
G4double GetDensity() const
Definition: G4Material.hh:180
long G4long
Definition: G4Types.hh:80
void NormShift(G4double betaGammaSq)
G4int GetMaxInterval() const
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
G4double SumOverBordResonance(G4int intervalNumber, G4double energy)
static G4int fNumberOfGammas
G4double GetSandiaMatTablePAI(G4int, G4int) const
G4double SumOverInterPlasmon(G4int intervalNumber)
G4double a
Definition: TRTMaterials.hh:39
G4double SumOverInterval(G4int intervalNumber)
G4double RePartDielectricConst(G4double energy)
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:202
void ComputeLowEnergyCof()
int G4int
Definition: G4Types.hh:78
void CallError(G4int i, const G4String &methodName) const
static const G4int fMaxSplineSize
G4int SandiaMixing(G4int Z[], const G4double *fractionW, G4int el, G4int mi)
static const G4double betaBohr4
G4double GetResonanceEnergyTransfer()
#define position
Definition: xmlparse.cc:605
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
G4bool GetLowerI1()
G4double GetPhotonRange(G4double energy)
G4double SumOverBordMM(G4int intervalNumber, G4double energy)
static const G4double betaBohr2
static const G4double c3
const G4double p2
const G4double p1
G4double GetCerenkovEnergyTransfer()
static const G4int fRefGammaNumber
static const G4double fError
static const G4double c1
G4double SumOverInterResonance(G4int intervalNumber)
G4int SandiaIntervals(G4int Z[], G4int el)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetStepPlasmonLoss(G4double step)
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
static const double eV
Definition: G4SIunits.hh:194
const G4double p0
G4double GetStepEnergyLoss(G4double step)
G4double SumOverBordCerenkov(G4int intervalNumber, G4double energy)
G4double GetMMEnergyTransfer()
G4double SumOverBorder(G4int intervalNumber, G4double energy)
G4double energy(const ThreeVector &p, const G4double m)
static const double g
Definition: G4SIunits.hh:162
static const G4double fDelta
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
G4double SumOverInterCerenkov(G4int intervalNumber)
#define DBL_MIN
Definition: templates.hh:75
G4double GetPhotoAbsorpCof(G4int i, G4int j) const
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
static const double keV
Definition: G4SIunits.hh:195
double G4double
Definition: G4Types.hh:76
G4double GetEnergyTransfer()
G4double GetStepMMLoss(G4double step)
G4double SumOverIntervaldEdx(G4int intervalNumber)
#define DBL_MAX
Definition: templates.hh:83
G4double GetStepCerenkovLoss(G4double step)
G4double SumOverBorderdEdx(G4int intervalNumber, G4double energy)
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)
static const G4double fLorentzFactor[112]
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
const G4Material * GetMaterial() const
G4double GetElectronRange(G4double energy)