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