Geant4  10.01.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 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 
92 G4PAIxSection::G4PAIxSection()
93 {
94  fSandia = 0;
95  fMatSandiaMatrix = 0;
96  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
97  fIntervalNumber = fSplineNumber = 0;
98  fVerbose = 0;
99 
100  fSplineEnergy = G4DataVector(fMaxSplineSize,0.0);
101  fRePartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
102  fImPartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
103  fIntegralTerm = G4DataVector(fMaxSplineSize,0.0);
104  fDifPAIxSection = G4DataVector(fMaxSplineSize,0.0);
105  fdNdxCerenkov = G4DataVector(fMaxSplineSize,0.0);
106  fdNdxPlasmon = G4DataVector(fMaxSplineSize,0.0);
107  fdNdxMM = G4DataVector(fMaxSplineSize,0.0);
108  fdNdxResonance = G4DataVector(fMaxSplineSize,0.0);
109  fIntegralPAIxSection = G4DataVector(fMaxSplineSize,0.0);
110  fIntegralPAIdEdx = G4DataVector(fMaxSplineSize,0.0);
111  fIntegralCerenkov = G4DataVector(fMaxSplineSize,0.0);
112  fIntegralPlasmon = G4DataVector(fMaxSplineSize,0.0);
113  fIntegralMM = G4DataVector(fMaxSplineSize,0.0);
114  fIntegralResonance = G4DataVector(fMaxSplineSize,0.0);
115 
116  fMaterialIndex = 0;
117 
118  for( G4int i = 0; i < 500; ++i )
119  {
120  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
121  }
122 }
123 
125 //
126 // Constructor
127 //
128 
129 G4PAIxSection::G4PAIxSection(G4MaterialCutsCouple* matCC)
130 {
131  fDensity = matCC->GetMaterial()->GetDensity();
132  G4int matIndex = matCC->GetMaterial()->GetIndex();
133  fMaterialIndex = matIndex;
134 
135  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
136  fSandia = (*theMaterialTable)[matIndex]->GetSandiaTable();
137 
138  fVerbose = 0;
139 
140  G4int i, j;
141  fMatSandiaMatrix = new G4OrderedTable();
142 
143  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
144  {
145  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
146  }
147  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
148  {
149  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
150 
151  for(j = 1; j < 5; j++)
152  {
153  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
154  }
155  }
156  ComputeLowEnergyCof();
157  // fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0;
158 }
159 
161 
162 G4PAIxSection::G4PAIxSection(G4int materialIndex,
163  G4double maxEnergyTransfer)
164 {
165  fSandia = 0;
166  fMatSandiaMatrix = 0;
167  fVerbose = 0;
168  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
169  G4int i, j;
170 
171  fMaterialIndex = materialIndex;
172  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
173  fElectronDensity = (*theMaterialTable)[materialIndex]->
174  GetElectronDensity();
175  fIntervalNumber = (*theMaterialTable)[materialIndex]->
176  GetSandiaTable()->GetMatNbOfIntervals();
177  fIntervalNumber--;
178  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
179 
180  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
181  fA1 = G4DataVector(fIntervalNumber+2,0.0);
182  fA2 = G4DataVector(fIntervalNumber+2,0.0);
183  fA3 = G4DataVector(fIntervalNumber+2,0.0);
184  fA4 = G4DataVector(fIntervalNumber+2,0.0);
185 
186  for(i = 1; i <= fIntervalNumber; i++ )
187  {
188  if(((*theMaterialTable)[materialIndex]->
189  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
190  i > fIntervalNumber )
191  {
192  fEnergyInterval[i] = maxEnergyTransfer;
193  fIntervalNumber = i;
194  break;
195  }
196  fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
197  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
198  fA1[i] = (*theMaterialTable)[materialIndex]->
199  GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
200  fA2[i] = (*theMaterialTable)[materialIndex]->
201  GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
202  fA3[i] = (*theMaterialTable)[materialIndex]->
203  GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
204  fA4[i] = (*theMaterialTable)[materialIndex]->
205  GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
206  // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
207  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
208  }
209  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
210  {
211  fIntervalNumber++;
212  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
213  }
214 
215  // Now checking, if two borders are too close together
216 
217  for(i=1;i<fIntervalNumber;i++)
218  {
219  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
220  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
221  {
222  continue;
223  }
224  else
225  {
226  for(j=i;j<fIntervalNumber;j++)
227  {
228  fEnergyInterval[j] = fEnergyInterval[j+1];
229  fA1[j] = fA1[j+1];
230  fA2[j] = fA2[j+1];
231  fA3[j] = fA3[j+1];
232  fA4[j] = fA4[j+1];
233  }
234  fIntervalNumber--;
235  i--;
236  }
237  }
238 
239 
240  /* *********************************
241 
242  fSplineEnergy = new G4double[fMaxSplineSize];
243  fRePartDielectricConst = new G4double[fMaxSplineSize];
244  fImPartDielectricConst = new G4double[fMaxSplineSize];
245  fIntegralTerm = new G4double[fMaxSplineSize];
246  fDifPAIxSection = new G4double[fMaxSplineSize];
247  fIntegralPAIxSection = new G4double[fMaxSplineSize];
248 
249  for(i=0;i<fMaxSplineSize;i++)
250  {
251  fSplineEnergy[i] = 0.0;
252  fRePartDielectricConst[i] = 0.0;
253  fImPartDielectricConst[i] = 0.0;
254  fIntegralTerm[i] = 0.0;
255  fDifPAIxSection[i] = 0.0;
256  fIntegralPAIxSection[i] = 0.0;
257  }
258  ************************************************** */
259  ComputeLowEnergyCof();
260  InitPAI(); // create arrays allocated above
261  /*
262  delete[] fEnergyInterval;
263  delete[] fA1;
264  delete[] fA2;
265  delete[] fA3;
266  delete[] fA4;
267  */
268 }
269 
271 //
272 // Constructor called from G4PAIPhotonModel !!!
273 
274 G4PAIxSection::G4PAIxSection( G4int materialIndex,
275  G4double maxEnergyTransfer,
276  G4double betaGammaSq,
277  G4double** photoAbsCof,
278  G4int intNumber )
279 {
280  fSandia = 0;
281  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
282  fIntervalNumber = fSplineNumber = 0;
283  fVerbose = 0;
284 
285  fSplineEnergy = G4DataVector(500,0.0);
286  fRePartDielectricConst = G4DataVector(500,0.0);
287  fImPartDielectricConst = G4DataVector(500,0.0);
288  fIntegralTerm = G4DataVector(500,0.0);
289  fDifPAIxSection = G4DataVector(500,0.0);
290  fdNdxCerenkov = G4DataVector(500,0.0);
291  fdNdxPlasmon = G4DataVector(500,0.0);
292  fdNdxMM = G4DataVector(500,0.0);
293  fdNdxResonance = G4DataVector(500,0.0);
294  fIntegralPAIxSection = G4DataVector(500,0.0);
295  fIntegralPAIdEdx = G4DataVector(500,0.0);
296  fIntegralCerenkov = G4DataVector(500,0.0);
297  fIntegralPlasmon = G4DataVector(500,0.0);
298  fIntegralMM = G4DataVector(500,0.0);
299  fIntegralResonance = G4DataVector(500,0.0);
300 
301  for( G4int i = 0; i < 500; ++i )
302  {
303  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
304  }
305 
306  fSandia = 0;
307  fMatSandiaMatrix = 0;
308  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
309  G4int i, j;
310 
311  fMaterialIndex = materialIndex;
312  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
313  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
314 
315  fIntervalNumber = intNumber;
316  fIntervalNumber--;
317  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
318 
319  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
320  fA1 = G4DataVector(fIntervalNumber+2,0.0);
321  fA2 = G4DataVector(fIntervalNumber+2,0.0);
322  fA3 = G4DataVector(fIntervalNumber+2,0.0);
323  fA4 = G4DataVector(fIntervalNumber+2,0.0);
324 
325 
326  /*
327  fEnergyInterval = new G4double[fIntervalNumber+2];
328  fA1 = new G4double[fIntervalNumber+2];
329  fA2 = new G4double[fIntervalNumber+2];
330  fA3 = new G4double[fIntervalNumber+2];
331  fA4 = new G4double[fIntervalNumber+2];
332  */
333  for( i = 1; i <= fIntervalNumber; i++ )
334  {
335  if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
336  i > fIntervalNumber )
337  {
338  fEnergyInterval[i] = maxEnergyTransfer;
339  fIntervalNumber = i;
340  break;
341  }
342  fEnergyInterval[i] = photoAbsCof[i-1][0];
343  fA1[i] = photoAbsCof[i-1][1];
344  fA2[i] = photoAbsCof[i-1][2];
345  fA3[i] = photoAbsCof[i-1][3];
346  fA4[i] = photoAbsCof[i-1][4];
347  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
348  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
349  }
350  // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
351 
352  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
353  {
354  fIntervalNumber++;
355  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
356  }
357  // G4cout<<"after check of max energy transfer"<<G4endl;
358 
359  for( i = 1; i <= fIntervalNumber; i++ )
360  {
361  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
362  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
363  }
364  // Now checking, if two borders are too close together
365 
366  for( i = 1; i < fIntervalNumber; i++ )
367  {
368  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
369  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
370  {
371  continue;
372  }
373  else
374  {
375  for(j=i;j<fIntervalNumber;j++)
376  {
377  fEnergyInterval[j] = fEnergyInterval[j+1];
378  fA1[j] = fA1[j+1];
379  fA2[j] = fA2[j+1];
380  fA3[j] = fA3[j+1];
381  fA4[j] = fA4[j+1];
382  }
383  fIntervalNumber--;
384  i--;
385  }
386  }
387  // G4cout<<"after check of close borders"<<G4endl;
388 
389  for( i = 1; i <= fIntervalNumber; i++ )
390  {
391  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
392  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
393  }
394 
395  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
396 
397  ComputeLowEnergyCof();
398  G4double betaGammaSqRef =
399  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
400 
401  NormShift(betaGammaSqRef);
402  SplainPAI(betaGammaSqRef);
403 
404  // Preparation of integral PAI cross section for input betaGammaSq
405 
406  for(i = 1; i <= fSplineNumber; i++)
407  {
408  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
409  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
410  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
411  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
412  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
413 
414  // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
415  // <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
416  }
417  IntegralCerenkov();
418  IntegralMM();
419  IntegralPlasmon();
420  IntegralResonance();
421  IntegralPAIxSection();
422  /*
423  delete[] fEnergyInterval;
424  delete[] fA1;
425  delete[] fA2;
426  delete[] fA3;
427  delete[] fA4;
428  */
429 }
430 
432 //
433 // Test Constructor with beta*gamma square value
434 
435 G4PAIxSection::G4PAIxSection( G4int materialIndex,
436  G4double maxEnergyTransfer,
437  G4double betaGammaSq )
438 {
439  fSandia = 0;
440  fMatSandiaMatrix = 0;
441  fVerbose = 0;
442  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
443 
444  G4int i, j, numberOfElements;
445 
446  fMaterialIndex = materialIndex;
447  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
448  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
449  numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements();
450 
451  G4int* thisMaterialZ = new G4int[numberOfElements];
452 
453  for( i = 0; i < numberOfElements; i++ )
454  {
455  thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
456  GetElement(i)->GetZ();
457  }
458  // fSandia = new G4SandiaTable(materialIndex);
459  fSandia = (*theMaterialTable)[materialIndex]->GetSandiaTable();
460  G4SandiaTable thisMaterialSandiaTable(materialIndex);
461  fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals(thisMaterialZ,
462  numberOfElements);
463  fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
464  ( thisMaterialZ ,
465  (*theMaterialTable)[materialIndex]->GetFractionVector() ,
466  numberOfElements,fIntervalNumber);
467 
468  fIntervalNumber--;
469 
470  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
471  fA1 = G4DataVector(fIntervalNumber+2,0.0);
472  fA2 = G4DataVector(fIntervalNumber+2,0.0);
473  fA3 = G4DataVector(fIntervalNumber+2,0.0);
474  fA4 = G4DataVector(fIntervalNumber+2,0.0);
475 
476  /*
477  fEnergyInterval = new G4double[fIntervalNumber+2];
478  fA1 = new G4double[fIntervalNumber+2];
479  fA2 = new G4double[fIntervalNumber+2];
480  fA3 = new G4double[fIntervalNumber+2];
481  fA4 = new G4double[fIntervalNumber+2];
482  */
483  for( i = 1; i <= fIntervalNumber; i++ )
484  {
485  if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
486  i > fIntervalNumber)
487  {
488  fEnergyInterval[i] = maxEnergyTransfer;
489  fIntervalNumber = i;
490  break;
491  }
492  fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0);
493  fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity;
494  fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity;
495  fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity;
496  fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity;
497 
498  }
499  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
500  {
501  fIntervalNumber++;
502  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
503  fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
504  fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
505  fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
506  fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
507  }
508  for(i=1;i<=fIntervalNumber;i++)
509  {
510  // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
511  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
512  }
513  // Now checking, if two borders are too close together
514 
515  for( i = 1; i < fIntervalNumber; i++ )
516  {
517  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
518  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
519  {
520  continue;
521  }
522  else
523  {
524  for( j = i; j < fIntervalNumber; j++ )
525  {
526  fEnergyInterval[j] = fEnergyInterval[j+1];
527  fA1[j] = fA1[j+1];
528  fA2[j] = fA2[j+1];
529  fA3[j] = fA3[j+1];
530  fA4[j] = fA4[j+1];
531  }
532  fIntervalNumber--;
533  i--;
534  }
535  }
536 
537  /* *********************************
538  fSplineEnergy = new G4double[fMaxSplineSize];
539  fRePartDielectricConst = new G4double[fMaxSplineSize];
540  fImPartDielectricConst = new G4double[fMaxSplineSize];
541  fIntegralTerm = new G4double[fMaxSplineSize];
542  fDifPAIxSection = new G4double[fMaxSplineSize];
543  fIntegralPAIxSection = new G4double[fMaxSplineSize];
544 
545  for(i=0;i<fMaxSplineSize;i++)
546  {
547  fSplineEnergy[i] = 0.0;
548  fRePartDielectricConst[i] = 0.0;
549  fImPartDielectricConst[i] = 0.0;
550  fIntegralTerm[i] = 0.0;
551  fDifPAIxSection[i] = 0.0;
552  fIntegralPAIxSection[i] = 0.0;
553  }
554  */
555 
556  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
557 
558  ComputeLowEnergyCof();
559  G4double betaGammaSqRef =
560  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
561 
562  NormShift(betaGammaSqRef);
563  SplainPAI(betaGammaSqRef);
564 
565  // Preparation of integral PAI cross section for input betaGammaSq
566 
567  for(i = 1; i <= fSplineNumber; i++)
568  {
569  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
570  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
571  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
572  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
573  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
574  }
575  IntegralPAIxSection();
576  IntegralCerenkov();
577  IntegralMM();
578  IntegralPlasmon();
579  IntegralResonance();
580 }
581 
583 //
584 // Destructor
585 
586 G4PAIxSection::~G4PAIxSection()
587 {
588  /* ************************
589  delete[] fSplineEnergy ;
590  delete[] fRePartDielectricConst;
591  delete[] fImPartDielectricConst;
592  delete[] fIntegralTerm ;
593  delete[] fDifPAIxSection ;
594  delete[] fIntegralPAIxSection ;
595  */
596  delete fMatSandiaMatrix;
597 }
598 
599 
600 
601 
603 //
604 // Constructor with beta*gamma square value called from G4PAIPhotModel/Data
605 
606 void G4PAIxSection::Initialize( const G4Material* material,
607  G4double maxEnergyTransfer,
608  G4double betaGammaSq,
609  G4SandiaTable* sandia)
610 {
611  if(fVerbose > 0)
612  {
613  G4cout<<G4endl;
614  G4cout<<"G4PAIxSection::Initialize(...,G4SandiaTable* sandia)"<<G4endl;
615  G4cout<<G4endl;
616  }
617  G4int i, j;
618 
619  fSandia = sandia;
620  fIntervalNumber = sandia->GetMaxInterval();
621  fDensity = material->GetDensity();
622  fElectronDensity = material->GetElectronDensity();
623 
624  // fIntervalNumber--;
625 
626  if( fVerbose > 0 )
627  {
628  G4cout<<"fDensity = "<<fDensity<<"\t"<<fElectronDensity<<"\t fIntervalNumber = "<<fIntervalNumber<<G4endl;
629  }
630  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
631  fA1 = G4DataVector(fIntervalNumber+2,0.0);
632  fA2 = G4DataVector(fIntervalNumber+2,0.0);
633  fA3 = G4DataVector(fIntervalNumber+2,0.0);
634  fA4 = G4DataVector(fIntervalNumber+2,0.0);
635 
636  for( i = 1; i <= fIntervalNumber; i++ )
637  {
638  if ( sandia->GetSandiaMatTablePAI(i-1,0) < 1.*eV && sandia->GetLowerI1() == false)
639  {
640  fIntervalNumber--;
641  continue;
642  }
643  if( ( sandia->GetSandiaMatTablePAI(i-1,0) >= maxEnergyTransfer ) || i >= fIntervalNumber )
644  {
645  fEnergyInterval[i] = maxEnergyTransfer;
646  fIntervalNumber = i;
647  break;
648  }
649  fEnergyInterval[i] = sandia->GetSandiaMatTablePAI(i-1,0);
650  fA1[i] = sandia->GetSandiaMatTablePAI(i-1,1);
651  fA2[i] = sandia->GetSandiaMatTablePAI(i-1,2);
652  fA3[i] = sandia->GetSandiaMatTablePAI(i-1,3);
653  fA4[i] = sandia->GetSandiaMatTablePAI(i-1,4);
654 
655  if( fVerbose > 0 )
656  {
657  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
658  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
659  }
660  }
661  if( fVerbose > 0 ) G4cout<<"last i = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
662 
663  if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer )
664  {
665  fIntervalNumber++;
666  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
667  }
668  if( fVerbose > 0 )
669  {
670  for( i = 1; i <= fIntervalNumber; i++ )
671  {
672  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
673  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
674  }
675  }
676  if( fVerbose > 0 ) G4cout<<"Now checking, if two borders are too close together"<<G4endl;
677 
678  for( i = 1; i < fIntervalNumber; i++ )
679  {
680  if( fEnergyInterval[i+1]-fEnergyInterval[i] >
681  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]) && fEnergyInterval[i] > 0.) continue;
682  else
683  {
684  if( fVerbose > 0 ) G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fEnergyInterval[i+1]/keV;
685 
686  for( j = i; j < fIntervalNumber; j++ )
687  {
688  fEnergyInterval[j] = fEnergyInterval[j+1];
689  fA1[j] = fA1[j+1];
690  fA2[j] = fA2[j+1];
691  fA3[j] = fA3[j+1];
692  fA4[j] = fA4[j+1];
693  }
694  fIntervalNumber--;
695  i--;
696  }
697  }
698  if( fVerbose > 0 )
699  {
700  for( i = 1; i <= fIntervalNumber; i++ )
701  {
702  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
703  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
704  }
705  }
706  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
707 
708  ComputeLowEnergyCof(material);
709 
710  G4double betaGammaSqRef =
711  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
712 
713  NormShift(betaGammaSqRef);
714  SplainPAI(betaGammaSqRef);
715 
716  // Preparation of integral PAI cross section for input betaGammaSq
717 
718  for( i = 1; i <= fSplineNumber; i++ )
719  {
720  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
721 
722 
723  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
724  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
725  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
726  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
727  }
728  IntegralPAIxSection();
729  IntegralCerenkov();
730  IntegralMM();
731  IntegralPlasmon();
732  IntegralResonance();
733 
734  for( i = 1; i <= fSplineNumber; i++ )
735  {
736  if(fVerbose>0) G4cout<<i<<"; w = "<<fSplineEnergy[i]/keV<<" keV; dN/dx_>w = "<<fIntegralPAIxSection[i]<<" 1/mm"<<G4endl;
737  }
738 }
739 
740 
742 //
743 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
744 //
745 
746 void G4PAIxSection::ComputeLowEnergyCof(const G4Material* material)
747 {
748  G4int i, numberOfElements = material->GetNumberOfElements();
749  G4double sumZ = 0., sumCof = 0.;
750 
751  static const G4double p0 = 1.20923e+00;
752  static const G4double p1 = 3.53256e-01;
753  static const G4double p2 = -1.45052e-03;
754 
755  G4double* thisMaterialZ = new G4double[numberOfElements];
756  G4double* thisMaterialCof = new G4double[numberOfElements];
757 
758  for( i = 0; i < numberOfElements; i++ )
759  {
760  thisMaterialZ[i] = material->GetElement(i)->GetZ();
761  sumZ += thisMaterialZ[i];
762  thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
763  }
764  for( i = 0; i < numberOfElements; i++ )
765  {
766  sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
767  }
768  fLowEnergyCof = sumCof;
769  delete [] thisMaterialZ;
770  delete [] thisMaterialCof;
771  // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
772 }
773 
775 //
776 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
777 //
778 
779 void G4PAIxSection::ComputeLowEnergyCof()
780 {
781  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
782  G4int i, numberOfElements = (*theMaterialTable)[fMaterialIndex]->GetNumberOfElements();
783  G4double sumZ = 0., sumCof = 0.;
784 
785  const G4double p0 = 1.20923e+00;
786  const G4double p1 = 3.53256e-01;
787  const G4double p2 = -1.45052e-03;
788 
789  G4double* thisMaterialZ = new G4double[numberOfElements];
790  G4double* thisMaterialCof = new G4double[numberOfElements];
791 
792  for( i = 0; i < numberOfElements; i++ )
793  {
794  thisMaterialZ[i] = (*theMaterialTable)[fMaterialIndex]->GetElement(i)->GetZ();
795  sumZ += thisMaterialZ[i];
796  thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
797  }
798  for( i = 0; i < numberOfElements; i++ )
799  {
800  sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
801  }
802  fLowEnergyCof = sumCof;
803  // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
804  delete [] thisMaterialZ;
805  delete [] thisMaterialCof;
806 }
807 
809 //
810 // General control function for class G4PAIxSection
811 //
812 
813 void G4PAIxSection::InitPAI()
814 {
815  G4int i;
816  G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
817  fLorentzFactor[fRefGammaNumber] - 1;
818 
819  // Preparation of integral PAI cross section for reference gamma
820 
821  NormShift(betaGammaSq);
822  SplainPAI(betaGammaSq);
823 
824  IntegralPAIxSection();
825  IntegralCerenkov();
826  IntegralMM();
827  IntegralPlasmon();
828  IntegralResonance();
829 
830  for(i = 0; i<= fSplineNumber; i++)
831  {
832  fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i];
833  if(i != 0)
834  {
835  fPAItable[i][0] = fSplineEnergy[i];
836  }
837  }
838  fPAItable[0][0] = fSplineNumber;
839 
840  for(G4int j = 1; j < 112; j++) // for other gammas
841  {
842  if( j == fRefGammaNumber ) continue;
843 
844  betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
845 
846  for(i = 1; i <= fSplineNumber; i++)
847  {
848  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
849  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
850  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
851  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
852  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
853  }
854  IntegralPAIxSection();
855  IntegralCerenkov();
856  IntegralMM();
857  IntegralPlasmon();
858  IntegralResonance();
859 
860  for(i = 0; i <= fSplineNumber; i++)
861  {
862  fPAItable[i][j] = fIntegralPAIxSection[i];
863  }
864  }
865 
866 }
867 
869 //
870 // Shifting from borders to intervals Creation of first energy points
871 //
872 
873 void G4PAIxSection::NormShift(G4double betaGammaSq)
874 {
875  G4int i, j;
876 
877  if(fVerbose>0) G4cout<<" G4PAIxSection::NormShift call "<<G4endl;
878 
879 
880  for( i = 1; i <= fIntervalNumber-1; i++ )
881  {
882  for( j = 1; j <= 2; j++ )
883  {
884  fSplineNumber = (i-1)*2 + j;
885 
886  if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
887  else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
888  if(fVerbose>0) G4cout<<"cn = "<<fSplineNumber<<"; "<<"w = "<<fSplineEnergy[fSplineNumber]/keV<<" keV"<<G4endl;
889  }
890  }
891  fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
892 
893  j = 1;
894 
895  for( i = 2; i <= fSplineNumber; i++ )
896  {
897  if( fSplineEnergy[i]<fEnergyInterval[j+1] )
898  {
899  fIntegralTerm[i] = fIntegralTerm[i-1] +
900  RutherfordIntegral(j,fSplineEnergy[i-1],
901  fSplineEnergy[i] );
902  }
903  else
904  {
905  G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
906  fEnergyInterval[j+1] );
907  j++;
908  fIntegralTerm[i] = fIntegralTerm[i-1] + x +
909  RutherfordIntegral(j,fEnergyInterval[j],
910  fSplineEnergy[i] );
911  }
912  if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV \t"<<fIntegralTerm[i]<<"\n"<<G4endl;
913  }
914  fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
915  fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
916 
917  // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl;
918 
919  // Calculation of PAI differrential cross-section (1/(keV*cm))
920  // in the energy points near borders of energy intervals
921 
922  for(G4int k = 1; k <= fIntervalNumber-1; k++ )
923  {
924  for( j = 1; j <= 2; j++ )
925  {
926  i = (k-1)*2 + j;
927  fImPartDielectricConst[i] = fNormalizationCof*
928  ImPartDielectricConst(k,fSplineEnergy[i]);
929  fRePartDielectricConst[i] = fNormalizationCof*
930  RePartDielectricConst(fSplineEnergy[i]);
931  fIntegralTerm[i] *= fNormalizationCof;
932 
933  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
934  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
935  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
936  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
937  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
938  if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV, xsc = "<<fDifPAIxSection[i]<<"\n"<<G4endl;
939  }
940  }
941 
942 } // end of NormShift
943 
945 //
946 // Creation of new energy points as geometrical mean of existing
947 // one, calculation PAI_cs for them, while the error of logarithmic
948 // linear approximation would be smaller than 'fError'
949 
950 void G4PAIxSection::SplainPAI(G4double betaGammaSq)
951 {
952  G4int j, k = 1, i = 1;
953 
954  if(fVerbose>0) G4cout<<" G4PAIxSection::SplainPAI call "<<G4endl;
955 
956  while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
957  {
958  // if( std::abs(fSplineEnergy[i+1]-fEnergyInterval[k+1]) > (fSplineEnergy[i+1]+fEnergyInterval[k+1])*5.e-7 )
959  if( fSplineEnergy[i+1] > fEnergyInterval[k+1] )
960  {
961  k++; // Here next energy point is in next energy interval
962  i++;
963  if(fVerbose>0) G4cout<<" in if: i = "<<i<<"; k = "<<k<<G4endl;
964  continue;
965  }
966  if(fVerbose>0) G4cout<<" out if: i = "<<i<<"; k = "<<k<<G4endl;
967 
968  // Shifting of arrayes for inserting the geometrical
969  // average of 'i' and 'i+1' energy points to 'i+1' place
970  fSplineNumber++;
971 
972  for( j = fSplineNumber; j >= i+2; j-- )
973  {
974  fSplineEnergy[j] = fSplineEnergy[j-1];
975  fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
976  fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
977  fIntegralTerm[j] = fIntegralTerm[j-1];
978 
979  fDifPAIxSection[j] = fDifPAIxSection[j-1];
980  fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
981  fdNdxMM[j] = fdNdxMM[j-1];
982  fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
983  fdNdxResonance[j] = fdNdxResonance[j-1];
984  }
985  G4double x1 = fSplineEnergy[i];
986  G4double x2 = fSplineEnergy[i+1];
987  G4double yy1 = fDifPAIxSection[i];
988  G4double y2 = fDifPAIxSection[i+1];
989 
990  if(fVerbose>0) G4cout<<"Spline: x1 = "<<x1<<"; x2 = "<<x2<<", yy1 = "<<yy1<<"; y2 = "<<y2<<G4endl;
991 
992 
993  G4double en1 = sqrt(x1*x2);
994  // G4double en1 = 0.5*(x1 + x2);
995 
996 
997  fSplineEnergy[i+1] = en1;
998 
999  // Calculation of logarithmic linear approximation
1000  // in this (enr) energy point, which number is 'i+1' now
1001 
1002  G4double a = log10(y2/yy1)/log10(x2/x1);
1003  G4double b = log10(yy1) - a*log10(x1);
1004  G4double y = a*log10(en1) + b;
1005 
1006  y = pow(10.,y);
1007 
1008  // Calculation of the PAI dif. cross-section at this point
1009 
1010  fImPartDielectricConst[i+1] = fNormalizationCof*
1011  ImPartDielectricConst(k,fSplineEnergy[i+1]);
1012  fRePartDielectricConst[i+1] = fNormalizationCof*
1013  RePartDielectricConst(fSplineEnergy[i+1]);
1014  fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
1015  RutherfordIntegral(k,fSplineEnergy[i],
1016  fSplineEnergy[i+1]);
1017 
1018  fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
1019  fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
1020  fdNdxMM[i+1] = PAIdNdxMM(i+1,betaGammaSq);
1021  fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
1022  fdNdxResonance[i+1] = PAIdNdxResonance(i+1,betaGammaSq);
1023 
1024  // Condition for next division of this segment or to pass
1025 
1026  if(fVerbose>0) G4cout<<"Spline, a = "<<a<<"; b = "<<b<<"; new xsc = "<<y<<"; compxsc = "<<fDifPAIxSection[i+1]<<G4endl;
1027 
1028  // to higher energies
1029 
1030  G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
1031 
1032  G4double delta = 2.*(fSplineEnergy[i+1]-fSplineEnergy[i])/(fSplineEnergy[i+1]+fSplineEnergy[i]);
1033 
1034  if( x < 0 )
1035  {
1036  x = -x;
1037  }
1038  if( x > fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*fDelta )
1039  {
1040  continue; // next division
1041  }
1042  i += 2; // pass to next segment
1043 
1044  } // 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  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 
1684 G4double G4PAIxSection::SumOverIntervaldEdx( G4int i )
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 
1718 G4double G4PAIxSection::SumOverInterCerenkov( G4int i )
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 
1754 G4double G4PAIxSection::SumOverInterMM( G4int i )
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 
1792 G4double G4PAIxSection::SumOverInterPlasmon( G4int i )
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 
1827 G4double G4PAIxSection::SumOverInterResonance( G4int i )
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 
1861 G4double G4PAIxSection::SumOverBorder( G4int i ,
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 
1934 G4double G4PAIxSection::SumOverBorderdEdx( G4int i ,
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 
1989 G4double G4PAIxSection::SumOverBordCerenkov( G4int i ,
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 
2054 G4double G4PAIxSection::SumOverBordMM( G4int i ,
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 
2119 G4double G4PAIxSection::SumOverBordPlasmon( G4int i ,
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 
2173 G4double G4PAIxSection::SumOverBordResonance( G4int i ,
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 
2226 G4double G4PAIxSection::GetStepEnergyLoss( G4double step )
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 
2252 G4double G4PAIxSection::GetEnergyTransfer()
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 
2279 G4double G4PAIxSection::GetStepCerenkovLoss( G4double step )
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 
2305 G4double G4PAIxSection::GetStepMMLoss( G4double step )
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 
2331 G4double G4PAIxSection::GetCerenkovEnergyTransfer()
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 
2358 G4double G4PAIxSection::GetMMEnergyTransfer()
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 
2385 G4double G4PAIxSection::GetStepPlasmonLoss( G4double step )
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 
2411 G4double G4PAIxSection::GetPlasmonEnergyTransfer()
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 
2438 G4double G4PAIxSection::GetStepResonanceLoss( G4double step )
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 
2465 G4double G4PAIxSection::GetResonanceEnergyTransfer()
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 
2493 G4double G4PAIxSection::GetRutherfordEnergyTransfer()
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 + "()";
2522  G4ExceptionDescription ed;
2523  ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
2524  G4Exception(head,"pai001",FatalException,ed);
2525 }
2526 
2528 //
2529 // Init array of Lorentz factors
2530 //
2531 
2532 G4int G4PAIxSection::fNumberOfGammas = 111;
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
2568 G4int G4PAIxSection::fRefGammaNumber = 29;
2569 
2570 
2571 //
2572 // end of G4PAIxSection implementation file
2573 //
2575 
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:2493
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:128
G4Poisson
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4PAIxSection::GetPlasmonEnergyTransfer
G4double GetPlasmonEnergyTransfer()
Definition: G4PAIxSection.cc:2411
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:2438
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:2119
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:262
G4PAIxSection::SumOverInterMM
G4double SumOverInterMM(G4int intervalNumber)
Definition: G4PAIxSection.cc:1754
betaBohr
static const G4double betaBohr
Definition: G4PAIySection.cc:69
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:68
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:588
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:180
G4long
long G4long
Definition: G4Types.hh:80
G4PAIxSection::NormShift
void NormShift(G4double betaGammaSq)
Definition: G4PAIxSection.cc:873
G4SandiaTable::GetMaxInterval
G4int GetMaxInterval() const
Definition: G4SandiaTable.cc:685
G4PAIxSection::PAIdNdxPlasmon
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIxSection.cc:1409
G4PAIxSection::SumOverBordResonance
G4double SumOverBordResonance(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:2173
G4PAIxSection::fNumberOfGammas
static G4int fNumberOfGammas
Definition: G4PAIxSection.hh:236
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int) const
Definition: G4SandiaTable.cc:1008
G4PAIxSection::SumOverInterPlasmon
G4double SumOverInterPlasmon(G4int intervalNumber)
Definition: G4PAIxSection.cc:1792
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:202
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:2519
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:834
betaBohr4
static const G4double betaBohr4
Definition: G4PAIySection.cc:72
G4PAIxSection.hh
G4PAIxSection::GetResonanceEnergyTransfer
G4double GetResonanceEnergyTransfer()
Definition: G4PAIxSection.cc:2465
position
#define position
Definition: xmlparse.cc:605
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:95
G4cout
G4GLOB_DLL std::ostream G4cout
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4PAIxSection::ImPartDielectricConst
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1075
G4SandiaTable::GetLowerI1
G4bool GetLowerI1()
Definition: G4SandiaTable.hh:175
G4PAIxSection::GetPhotonRange
G4double GetPhotonRange(G4double energy)
Definition: G4PAIxSection.cc:1095
G4PAIxSection::SumOverBordMM
G4double SumOverBordMM(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:2054
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:2331
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:1827
G4SandiaTable::SandiaIntervals
G4int SandiaIntervals(G4int Z[], G4int el)
Definition: G4SandiaTable.cc:738
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:2385
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:2226
G4PAIxSection::SumOverBordCerenkov
G4double SumOverBordCerenkov(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1989
G4PAIxSection::GetMMEnergyTransfer
G4double GetMMEnergyTransfer()
Definition: G4PAIxSection.cc:2358
G4PAIxSection::SumOverBorder
G4double SumOverBorder(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1861
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
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
G4PAIxSection::SumOverInterCerenkov
G4double SumOverInterCerenkov(G4int intervalNumber)
Definition: G4PAIxSection.cc:1718
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:711
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
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:2252
G4PAIxSection::GetStepMMLoss
G4double GetStepMMLoss(G4double step)
Definition: G4PAIxSection.cc:2305
G4PAIxSection::SumOverIntervaldEdx
G4double SumOverIntervaldEdx(G4int intervalNumber)
Definition: G4PAIxSection.cc:1684
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4PAIxSection::GetStepCerenkovLoss
G4double GetStepCerenkovLoss(G4double step)
Definition: G4PAIxSection.cc:2279
G4PAIxSection::SumOverBorderdEdx
G4double SumOverBorderdEdx(G4int intervalNumber, G4double energy)
Definition: G4PAIxSection.cc:1934
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