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