Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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$
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 = 500; // Max size of output spline
88  // arrays
89 
91 //
92 // Constructor
93 //
94 
95 G4PAIxSection::G4PAIxSection(G4MaterialCutsCouple* matCC)
96 {
97  fDensity = matCC->GetMaterial()->GetDensity();
98  G4int matIndex = matCC->GetMaterial()->GetIndex();
99  fMaterialIndex = matIndex;
100  fSandia = new G4SandiaTable(matIndex);
101 
102  G4int i, j;
103  fMatSandiaMatrix = new G4OrderedTable();
104 
105  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
106  {
107  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
108  }
109  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
110  {
111  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
112 
113  for(j = 1; j < 5; j++)
114  {
115  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
116  }
117  }
118  fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0;
119 }
120 
122 
123 G4PAIxSection::G4PAIxSection(G4int materialIndex,
124  G4double maxEnergyTransfer)
125 {
126  fSandia = 0;
127  fMatSandiaMatrix = 0;
128  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
129  G4int i, j;
130 
131  fMaterialIndex = materialIndex;
132  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
133  fElectronDensity = (*theMaterialTable)[materialIndex]->
134  GetElectronDensity();
135  fIntervalNumber = (*theMaterialTable)[materialIndex]->
136  GetSandiaTable()->GetMatNbOfIntervals();
137  fIntervalNumber--;
138  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
139 
140  fEnergyInterval = new G4double[fIntervalNumber+2];
141  fA1 = new G4double[fIntervalNumber+2];
142  fA2 = new G4double[fIntervalNumber+2];
143  fA3 = new G4double[fIntervalNumber+2];
144  fA4 = new G4double[fIntervalNumber+2];
145 
146  for(i = 1; i <= fIntervalNumber; i++ )
147  {
148  if(((*theMaterialTable)[materialIndex]->
149  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
150  i > fIntervalNumber )
151  {
152  fEnergyInterval[i] = maxEnergyTransfer;
153  fIntervalNumber = i;
154  break;
155  }
156  fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
157  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
158  fA1[i] = (*theMaterialTable)[materialIndex]->
159  GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
160  fA2[i] = (*theMaterialTable)[materialIndex]->
161  GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
162  fA3[i] = (*theMaterialTable)[materialIndex]->
163  GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
164  fA4[i] = (*theMaterialTable)[materialIndex]->
165  GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
166  // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
167  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
168  }
169  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
170  {
171  fIntervalNumber++;
172  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
173  }
174 
175  // Now checking, if two borders are too close together
176 
177  for(i=1;i<fIntervalNumber;i++)
178  {
179  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
180  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
181  {
182  continue;
183  }
184  else
185  {
186  for(j=i;j<fIntervalNumber;j++)
187  {
188  fEnergyInterval[j] = fEnergyInterval[j+1];
189  fA1[j] = fA1[j+1];
190  fA2[j] = fA2[j+1];
191  fA3[j] = fA3[j+1];
192  fA4[j] = fA4[j+1];
193  }
194  fIntervalNumber--;
195  i--;
196  }
197  }
198 
199 
200  /* *********************************
201 
202  fSplineEnergy = new G4double[fMaxSplineSize];
203  fRePartDielectricConst = new G4double[fMaxSplineSize];
204  fImPartDielectricConst = new G4double[fMaxSplineSize];
205  fIntegralTerm = new G4double[fMaxSplineSize];
206  fDifPAIxSection = new G4double[fMaxSplineSize];
207  fIntegralPAIxSection = new G4double[fMaxSplineSize];
208 
209  for(i=0;i<fMaxSplineSize;i++)
210  {
211  fSplineEnergy[i] = 0.0;
212  fRePartDielectricConst[i] = 0.0;
213  fImPartDielectricConst[i] = 0.0;
214  fIntegralTerm[i] = 0.0;
215  fDifPAIxSection[i] = 0.0;
216  fIntegralPAIxSection[i] = 0.0;
217  }
218  ************************************************** */
219 
220  InitPAI(); // create arrays allocated above
221 
222  delete[] fEnergyInterval;
223  delete[] fA1;
224  delete[] fA2;
225  delete[] fA3;
226  delete[] fA4;
227 }
228 
230 //
231 // Constructor with beta*gamma square value
232 
233 G4PAIxSection::G4PAIxSection( G4int materialIndex,
234  G4double maxEnergyTransfer,
235  G4double betaGammaSq,
236  G4double** photoAbsCof,
237  G4int intNumber )
238 {
239  fSandia = 0;
240  fMatSandiaMatrix = 0;
241  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
242  G4int i, j;
243 
244  fMaterialIndex = materialIndex;
245  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
246  fElectronDensity = (*theMaterialTable)[materialIndex]->
247  GetElectronDensity();
248 
249  fIntervalNumber = intNumber;
250  fIntervalNumber--;
251  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
252 
253  fEnergyInterval = new G4double[fIntervalNumber+2];
254  fA1 = new G4double[fIntervalNumber+2];
255  fA2 = new G4double[fIntervalNumber+2];
256  fA3 = new G4double[fIntervalNumber+2];
257  fA4 = new G4double[fIntervalNumber+2];
258 
259  for( i = 1; i <= fIntervalNumber; i++ )
260  {
261  if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
262  i > fIntervalNumber )
263  {
264  fEnergyInterval[i] = maxEnergyTransfer;
265  fIntervalNumber = i;
266  break;
267  }
268  fEnergyInterval[i] = photoAbsCof[i-1][0];
269  fA1[i] = photoAbsCof[i-1][1];
270  fA2[i] = photoAbsCof[i-1][2];
271  fA3[i] = photoAbsCof[i-1][3];
272  fA4[i] = photoAbsCof[i-1][4];
273  // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
274  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
275  }
276  // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
277  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
278  {
279  fIntervalNumber++;
280  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
281  }
282  for(i=1;i<=fIntervalNumber;i++)
283  {
284  // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
285  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
286  }
287  // Now checking, if two borders are too close together
288 
289  for( i = 1; i < fIntervalNumber; i++ )
290  {
291  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
292  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
293  {
294  continue;
295  }
296  else
297  {
298  for(j=i;j<fIntervalNumber;j++)
299  {
300  fEnergyInterval[j] = fEnergyInterval[j+1];
301  fA1[j] = fA1[j+1];
302  fA2[j] = fA2[j+1];
303  fA3[j] = fA3[j+1];
304  fA4[j] = fA4[j+1];
305  }
306  fIntervalNumber--;
307  i--;
308  }
309  }
310 
311  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
312 
313  G4double betaGammaSqRef =
314  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
315 
316  NormShift(betaGammaSqRef);
317  SplainPAI(betaGammaSqRef);
318 
319  // Preparation of integral PAI cross section for input betaGammaSq
320 
321  for(i = 1; i <= fSplineNumber; i++)
322  {
323  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
324  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
325  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
326  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
327  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
328 
329  // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
330  // <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
331  }
332  IntegralCerenkov();
333  IntegralMM();
334  IntegralPlasmon();
335  IntegralResonance();
336  IntegralPAIxSection();
337 
338  delete[] fEnergyInterval;
339  delete[] fA1;
340  delete[] fA2;
341  delete[] fA3;
342  delete[] fA4;
343 }
344 
346 //
347 // Test Constructor with beta*gamma square value
348 
349 G4PAIxSection::G4PAIxSection( G4int materialIndex,
350  G4double maxEnergyTransfer,
351  G4double betaGammaSq )
352 {
353  fSandia = 0;
354  fMatSandiaMatrix = 0;
355  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
356 
357  G4int i, j, numberOfElements;
358 
359  fMaterialIndex = materialIndex;
360  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
361  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
362  numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements();
363 
364  G4int* thisMaterialZ = new G4int[numberOfElements];
365 
366  for( i = 0; i < numberOfElements; i++ )
367  {
368  thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
369  GetElement(i)->GetZ();
370  }
371  // fSandia = new G4SandiaTable(materialIndex);
372  fSandia = (*theMaterialTable)[materialIndex]->
373  GetSandiaTable();
374  G4SandiaTable thisMaterialSandiaTable(materialIndex);
375  fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
376  (thisMaterialZ,numberOfElements);
377  fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
378  ( thisMaterialZ ,
379  (*theMaterialTable)[materialIndex]->GetFractionVector() ,
380  numberOfElements,fIntervalNumber);
381 
382  fIntervalNumber--;
383 
384  fEnergyInterval = new G4double[fIntervalNumber+2];
385  fA1 = new G4double[fIntervalNumber+2];
386  fA2 = new G4double[fIntervalNumber+2];
387  fA3 = new G4double[fIntervalNumber+2];
388  fA4 = new G4double[fIntervalNumber+2];
389 
390  for( i = 1; i <= fIntervalNumber; i++ )
391  {
392  if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
393  i > fIntervalNumber)
394  {
395  fEnergyInterval[i] = maxEnergyTransfer;
396  fIntervalNumber = i;
397  break;
398  }
399  fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0);
400  fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity;
401  fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity;
402  fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity;
403  fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity;
404 
405  }
406  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
407  {
408  fIntervalNumber++;
409  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
410  fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
411  fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
412  fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
413  fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
414  }
415  for(i=1;i<=fIntervalNumber;i++)
416  {
417  // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
418  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
419  }
420  // Now checking, if two borders are too close together
421 
422  for( i = 1; i < fIntervalNumber; i++ )
423  {
424  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
425  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
426  {
427  continue;
428  }
429  else
430  {
431  for( j = i; j < fIntervalNumber; j++ )
432  {
433  fEnergyInterval[j] = fEnergyInterval[j+1];
434  fA1[j] = fA1[j+1];
435  fA2[j] = fA2[j+1];
436  fA3[j] = fA3[j+1];
437  fA4[j] = fA4[j+1];
438  }
439  fIntervalNumber--;
440  i--;
441  }
442  }
443 
444  /* *********************************
445  fSplineEnergy = new G4double[fMaxSplineSize];
446  fRePartDielectricConst = new G4double[fMaxSplineSize];
447  fImPartDielectricConst = new G4double[fMaxSplineSize];
448  fIntegralTerm = new G4double[fMaxSplineSize];
449  fDifPAIxSection = new G4double[fMaxSplineSize];
450  fIntegralPAIxSection = new G4double[fMaxSplineSize];
451 
452  for(i=0;i<fMaxSplineSize;i++)
453  {
454  fSplineEnergy[i] = 0.0;
455  fRePartDielectricConst[i] = 0.0;
456  fImPartDielectricConst[i] = 0.0;
457  fIntegralTerm[i] = 0.0;
458  fDifPAIxSection[i] = 0.0;
459  fIntegralPAIxSection[i] = 0.0;
460  }
461  */
462 
463  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
464 
465  G4double betaGammaSqRef =
466  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
467 
468  NormShift(betaGammaSqRef);
469  SplainPAI(betaGammaSqRef);
470 
471  // Preparation of integral PAI cross section for input betaGammaSq
472 
473  for(i = 1; i <= fSplineNumber; i++)
474  {
475  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
476  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
477  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
478  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
479  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
480  }
481  IntegralPAIxSection();
482  IntegralCerenkov();
483  IntegralMM();
484  IntegralPlasmon();
485  IntegralResonance();
486 
487  delete[] fEnergyInterval;
488  delete[] fA1;
489  delete[] fA2;
490  delete[] fA3;
491  delete[] fA4;
492  delete[] thisMaterialZ;
493 }
494 
495 
497 //
498 // Destructor
499 
500 G4PAIxSection::~G4PAIxSection()
501 {
502  /* ************************
503  delete[] fSplineEnergy ;
504  delete[] fRePartDielectricConst;
505  delete[] fImPartDielectricConst;
506  delete[] fIntegralTerm ;
507  delete[] fDifPAIxSection ;
508  delete[] fIntegralPAIxSection ;
509  */
510  delete fSandia;
511  delete fMatSandiaMatrix;
512 }
513 
515 //
516 // General control function for class G4PAIxSection
517 //
518 
519 void G4PAIxSection::InitPAI()
520 {
521  G4int i;
522  G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
523  fLorentzFactor[fRefGammaNumber] - 1;
524 
525  // Preparation of integral PAI cross section for reference gamma
526 
527  NormShift(betaGammaSq);
528  SplainPAI(betaGammaSq);
529 
530  IntegralPAIxSection();
531  IntegralCerenkov();
532  IntegralMM();
533  IntegralPlasmon();
534  IntegralResonance();
535 
536  for(i = 0; i<= fSplineNumber; i++)
537  {
538  fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i];
539  if(i != 0)
540  {
541  fPAItable[i][0] = fSplineEnergy[i];
542  }
543  }
544  fPAItable[0][0] = fSplineNumber;
545 
546  for(G4int j = 1; j < 112; j++) // for other gammas
547  {
548  if( j == fRefGammaNumber ) continue;
549 
550  betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
551 
552  for(i = 1; i <= fSplineNumber; i++)
553  {
554  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
555  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
556  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
557  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
558  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
559  }
560  IntegralPAIxSection();
561  IntegralCerenkov();
562  IntegralMM();
563  IntegralPlasmon();
564  IntegralResonance();
565 
566  for(i = 0; i <= fSplineNumber; i++)
567  {
568  fPAItable[i][j] = fIntegralPAIxSection[i];
569  }
570  }
571 
572 }
573 
575 //
576 // Shifting from borders to intervals Creation of first energy points
577 //
578 
579 void G4PAIxSection::NormShift(G4double betaGammaSq)
580 {
581  G4int i, j;
582 
583  for( i = 1; i <= fIntervalNumber-1; i++ )
584  {
585  for( j = 1; j <= 2; j++ )
586  {
587  fSplineNumber = (i-1)*2 + j;
588 
589  if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
590  else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
591  // G4cout<<"cn = "<<fSplineNumber<<"; "<<"energy = "
592  // <<fSplineEnergy[fSplineNumber]<<G4endl;
593  }
594  }
595  fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
596 
597  j = 1;
598 
599  for( i = 2; i <= fSplineNumber; i++ )
600  {
601  if(fSplineEnergy[i]<fEnergyInterval[j+1])
602  {
603  fIntegralTerm[i] = fIntegralTerm[i-1] +
604  RutherfordIntegral(j,fSplineEnergy[i-1],
605  fSplineEnergy[i] );
606  }
607  else
608  {
609  G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
610  fEnergyInterval[j+1] );
611  j++;
612  fIntegralTerm[i] = fIntegralTerm[i-1] + x +
613  RutherfordIntegral(j,fEnergyInterval[j],
614  fSplineEnergy[i] );
615  }
616  // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t"<<fIntegralTerm[i]<<"\n"<<G4endl;
617  }
618  fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
619  fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
620 
621  // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl;
622 
623  // Calculation of PAI differrential cross-section (1/(keV*cm))
624  // in the energy points near borders of energy intervals
625 
626  for(G4int k = 1; k <= fIntervalNumber-1; k++ )
627  {
628  for( j = 1; j <= 2; j++ )
629  {
630  i = (k-1)*2 + j;
631  fImPartDielectricConst[i] = fNormalizationCof*
632  ImPartDielectricConst(k,fSplineEnergy[i]);
633  fRePartDielectricConst[i] = fNormalizationCof*
634  RePartDielectricConst(fSplineEnergy[i]);
635  fIntegralTerm[i] *= fNormalizationCof;
636 
637  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
638  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
639  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
640  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
641  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
642  }
643  }
644 
645 } // end of NormShift
646 
648 //
649 // Creation of new energy points as geometrical mean of existing
650 // one, calculation PAI_cs for them, while the error of logarithmic
651 // linear approximation would be smaller than 'fError'
652 
653 void G4PAIxSection::SplainPAI(G4double betaGammaSq)
654 {
655  G4int k = 1;
656  G4int i = 1;
657 
658  while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
659  {
660  if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
661  {
662  k++; // Here next energy point is in next energy interval
663  i++;
664  continue;
665  }
666  // Shifting of arrayes for inserting the geometrical
667  // average of 'i' and 'i+1' energy points to 'i+1' place
668  fSplineNumber++;
669 
670  for(G4int j = fSplineNumber; j >= i+2; j-- )
671  {
672  fSplineEnergy[j] = fSplineEnergy[j-1];
673  fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
674  fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
675  fIntegralTerm[j] = fIntegralTerm[j-1];
676 
677  fDifPAIxSection[j] = fDifPAIxSection[j-1];
678  fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
679  fdNdxMM[j] = fdNdxMM[j-1];
680  fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
681  fdNdxResonance[j] = fdNdxResonance[j-1];
682  }
683  G4double x1 = fSplineEnergy[i];
684  G4double x2 = fSplineEnergy[i+1];
685  G4double yy1 = fDifPAIxSection[i];
686  G4double y2 = fDifPAIxSection[i+1];
687 
688  G4double en1 = sqrt(x1*x2);
689  fSplineEnergy[i+1] = en1;
690 
691  // Calculation of logarithmic linear approximation
692  // in this (enr) energy point, which number is 'i+1' now
693 
694  G4double a = log10(y2/yy1)/log10(x2/x1);
695  G4double b = log10(yy1) - a*log10(x1);
696  G4double y = a*log10(en1) + b;
697  y = pow(10.,y);
698 
699  // Calculation of the PAI dif. cross-section at this point
700 
701  fImPartDielectricConst[i+1] = fNormalizationCof*
702  ImPartDielectricConst(k,fSplineEnergy[i+1]);
703  fRePartDielectricConst[i+1] = fNormalizationCof*
704  RePartDielectricConst(fSplineEnergy[i+1]);
705  fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
706  RutherfordIntegral(k,fSplineEnergy[i],
707  fSplineEnergy[i+1]);
708 
709  fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
710  fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
711  fdNdxMM[i+1] = PAIdNdxMM(i+1,betaGammaSq);
712  fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
713  fdNdxResonance[i+1] = PAIdNdxResonance(i+1,betaGammaSq);
714 
715  // Condition for next division of this segment or to pass
716  // to higher energies
717 
718  G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
719 
720  if( x < 0 )
721  {
722  x = -x;
723  }
724  if( x > fError && fSplineNumber < fMaxSplineSize-1 )
725  {
726  continue; // next division
727  }
728  i += 2; // pass to next segment
729 
730  } // close 'while'
731 
732 } // end of SplainPAI
733 
734 
736 //
737 // Integration over electrons that could be considered
738 // quasi-free at energy transfer of interest
739 
740 G4double G4PAIxSection::RutherfordIntegral( G4int k,
741  G4double x1,
742  G4double x2 )
743 {
744  G4double c1, c2, c3;
745  // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
746  c1 = (x2 - x1)/x1/x2;
747  c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2;
748  c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
749  // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
750 
751  return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
752 
753 } // end of RutherfordIntegral
754 
755 
757 //
758 // Imaginary part of dielectric constant
759 // (G4int k - interval number, G4double en1 - energy point)
760 
761 G4double G4PAIxSection::ImPartDielectricConst( G4int k ,
762  G4double energy1 )
763 {
764  G4double energy2,energy3,energy4,result;
765 
766  energy2 = energy1*energy1;
767  energy3 = energy2*energy1;
768  energy4 = energy3*energy1;
769 
770  result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
771  result *=hbarc/energy1;
772 
773  return result;
774 
775 } // end of ImPartDielectricConst
776 
778 //
779 // Returns lambda of photon with energy1 in current material
780 
781 G4double G4PAIxSection::GetPhotonRange( G4double energy1 )
782 {
783  G4int i;
784  G4double energy2, energy3, energy4, result, lambda;
785 
786  energy2 = energy1*energy1;
787  energy3 = energy2*energy1;
788  energy4 = energy3*energy1;
789 
790  // G4double* SandiaCof = fSandia->GetSandiaCofForMaterialPAI(energy1);
791  // result = SandiaCof[0]/energy1+SandiaCof[1]/energy2+SandiaCof[2]/energy3+SandiaCof[3]/energy4;
792  // result *= fDensity;
793 
794  for( i = 1; i <= fIntervalNumber; i++ )
795  {
796  if( energy1 < fEnergyInterval[i]) break;
797  }
798  i--;
799  if(i == 0) i = 1;
800 
801  result = fA1[i]/energy1+fA2[i]/energy2+fA3[i]/energy3+fA4[i]/energy4;
802 
803  if( result > DBL_MIN ) lambda = 1./result;
804  else lambda = DBL_MAX;
805 
806  return lambda;
807 }
808 
810 //
811 // Return lambda of electron with energy1 in current material
812 // parametrisation from NIM A554(2005)474-493
813 
814 G4double G4PAIxSection::GetElectronRange( G4double energy )
815 {
816  G4double range;
817  /*
818  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
819 
820  G4double Z = (*theMaterialTable)[fMaterialIndex]->GetIonisation()->GetZeffective();
821  G4double A = (*theMaterialTable)[fMaterialIndex]->GetA();
822 
823  energy /= keV; // energy in keV in parametrised formula
824 
825  if (energy < 10.)
826  {
827  range = 3.872e-3*A/Z;
828  range *= pow(energy,1.492);
829  }
830  else
831  {
832  range = 6.97e-3*pow(energy,1.6);
833  }
834  */
835  // Blum&Rolandi Particle Detection with Drift Chambers, p. 7
836 
837  G4double cofA = 5.37e-4*g/cm2/keV;
838  G4double cofB = 0.9815;
839  G4double cofC = 3.123e-3/keV;
840  // energy /= keV;
841 
842  range = cofA*energy*( 1 - cofB/(1 + cofC*energy) );
843 
844  // range *= g/cm2;
845  range /= fDensity;
846 
847  return range;
848 }
849 
851 //
852 // Real part of dielectric constant minus unit: epsilon_1 - 1
853 // (G4double enb - energy point)
854 //
855 
856 G4double G4PAIxSection::RePartDielectricConst(G4double enb)
857 {
858  G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
859  c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
860 
861  x0 = enb;
862  result = 0;
863 
864  for(G4int i=1;i<=fIntervalNumber-1;i++)
865  {
866  x1 = fEnergyInterval[i];
867  x2 = fEnergyInterval[i+1];
868  xx1 = x1 - x0;
869  xx2 = x2 - x0;
870  xx12 = xx2/xx1;
871 
872  if(xx12<0)
873  {
874  xx12 = -xx12;
875  }
876  xln1 = log(x2/x1);
877  xln2 = log(xx12);
878  xln3 = log((x2 + x0)/(x1 + x0));
879  x02 = x0*x0;
880  x03 = x02*x0;
881  x04 = x03*x0;
882  x05 = x04*x0;
883  c1 = (x2 - x1)/x1/x2;
884  c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2;
885  c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
886 
887  result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
888  result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
889  result -= fA3[i]*c2/2/x02;
890  result -= fA4[i]*c3/3/x02;
891 
892  cof1 = fA1[i]/x02 + fA3[i]/x04;
893  cof2 = fA2[i]/x03 + fA4[i]/x05;
894 
895  result += 0.5*(cof1 +cof2)*xln2;
896  result += 0.5*(cof1 - cof2)*xln3;
897  }
898  result *= 2*hbarc/pi;
899 
900  return result;
901 
902 } // end of RePartDielectricConst
903 
905 //
906 // PAI differential cross-section in terms of
907 // simplified Allison's equation
908 //
909 
910 G4double G4PAIxSection::DifPAIxSection( G4int i ,
911  G4double betaGammaSq )
912 {
913  G4double cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
914 
915  G4double betaBohr = fine_structure_const;
916  // G4double betaBohr2 = fine_structure_const*fine_structure_const;
917  // G4double betaBohr3 = betaBohr*betaBohr2; // *4.0;
918 
919  G4double be2 = betaGammaSq/(1 + betaGammaSq);
920  G4double beta = sqrt(be2);
921  // G4double be3 = beta*be2;
922 
923  cof = 1.;
924  x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
925 
926  if( betaGammaSq < 0.01 ) x2 = log(be2);
927  else
928  {
929  x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
930  (1/betaGammaSq - fRePartDielectricConst[i]) +
931  fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
932  }
933  if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
934  {
935  x6 = 0.;
936  }
937  else
938  {
939  x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
940  x5 = -1 - fRePartDielectricConst[i] +
941  be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
942  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
943 
944  x7 = atan2(fImPartDielectricConst[i],x3);
945  x6 = x5 * x7;
946  }
947  // if(fImPartDielectricConst[i] == 0) x6 = 0.;
948 
949  x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
950 
951  // if( x4 < 0.0 ) x4 = 0.0;
952 
953  x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
954  fImPartDielectricConst[i]*fImPartDielectricConst[i];
955 
956  result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
957 
958  if( result < 1.0e-8 ) result = 1.0e-8;
959 
960  result *= fine_structure_const/be2/pi;
961 
962  // low energy correction
963 
964  G4double lowCof = 4.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)?
965 
966  result *= (1 - exp(-beta/betaBohr/lowCof));
967 
968  // result *= (1 - exp(-be2/betaBohr2/lowCof));
969 
970  // result *= (1 - exp(-be3/betaBohr3/lowCof)); // ~ be for be<<betaBohr
971 
972  // result *= (1 - exp(-be4/betaBohr4/lowCof));
973 
974  if(fDensity >= 0.1)
975  {
976  result /= x8;
977  }
978  return result;
979 
980 } // end of DifPAIxSection
981 
983 //
984 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
985 
986 G4double G4PAIxSection::PAIdNdxCerenkov( G4int i ,
987  G4double betaGammaSq )
988 {
989  G4double logarithm, x3, x5, argument, modul2, dNdxC;
990  G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
991 
992  cofBetaBohr = 4.0;
993  betaBohr2 = fine_structure_const*fine_structure_const;
994  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
995 
996  be2 = betaGammaSq/(1 + betaGammaSq);
997  be4 = be2*be2;
998 
999  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1000  else
1001  {
1002  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1003  (1/betaGammaSq - fRePartDielectricConst[i]) +
1004  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1005  logarithm += log(1+1.0/betaGammaSq);
1006  }
1007 
1008  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1009  {
1010  argument = 0.0;
1011  }
1012  else
1013  {
1014  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1015  x5 = -1.0 - fRePartDielectricConst[i] +
1016  be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1017  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1018  if( x3 == 0.0 ) argument = 0.5*pi;
1019  else argument = atan2(fImPartDielectricConst[i],x3);
1020  argument *= x5 ;
1021  }
1022  dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
1023 
1024  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1025 
1026  dNdxC *= fine_structure_const/be2/pi;
1027 
1028  dNdxC *= (1-exp(-be4/betaBohr4));
1029 
1030  if(fDensity >= 0.1)
1031  {
1032  modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1033  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1034  dNdxC /= modul2;
1035  }
1036  return dNdxC;
1037 
1038 } // end of PAIdNdxCerenkov
1039 
1041 //
1042 // Calculation od dN/dx of collisions of MM with creation of Cerenkov pseudo-photons
1043 
1044 G4double G4PAIxSection::PAIdNdxMM( G4int i ,
1045  G4double betaGammaSq )
1046 {
1047  G4double logarithm, x3, x5, argument, dNdxC;
1048  G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
1049 
1050  cofBetaBohr = 4.0;
1051  betaBohr2 = fine_structure_const*fine_structure_const;
1052  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1053 
1054  be2 = betaGammaSq/(1 + betaGammaSq);
1055  be4 = be2*be2;
1056 
1057  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1058  else
1059  {
1060  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1061  (1/betaGammaSq - fRePartDielectricConst[i]) +
1062  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1063  logarithm += log(1+1.0/betaGammaSq);
1064  }
1065 
1066  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1067  {
1068  argument = 0.0;
1069  }
1070  else
1071  {
1072  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1073  x5 = be2*( 1.0 + fRePartDielectricConst[i] ) - 1.0;
1074  if( x3 == 0.0 ) argument = 0.5*pi;
1075  else argument = atan2(fImPartDielectricConst[i],x3);
1076  argument *= x5 ;
1077  }
1078  dNdxC = ( logarithm*fImPartDielectricConst[i]*be2 + argument )/hbarc;
1079 
1080  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1081 
1082  dNdxC *= fine_structure_const/be2/pi;
1083 
1084  dNdxC *= (1-exp(-be4/betaBohr4));
1085  return dNdxC;
1086 
1087 } // end of PAIdNdxMM
1088 
1090 //
1091 // Calculation od dN/dx of collisions with creation of longitudinal EM
1092 // excitations (plasmons, delta-electrons)
1093 
1094 G4double G4PAIxSection::PAIdNdxPlasmon( G4int i ,
1095  G4double betaGammaSq )
1096 {
1097  G4double resonance, modul2, dNdxP, cof = 1.;
1098  G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
1099 
1100 
1101  cofBetaBohr = 4.0;
1102  betaBohr2 = fine_structure_const*fine_structure_const;
1103  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1104 
1105  be2 = betaGammaSq/(1 + betaGammaSq);
1106  be4 = be2*be2;
1107 
1108  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1109  resonance *= fImPartDielectricConst[i]/hbarc;
1110 
1111 
1112  dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
1113 
1114  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1115 
1116  dNdxP *= fine_structure_const/be2/pi;
1117  dNdxP *= (1-exp(-be4/betaBohr4));
1118 
1119  if( fDensity >= 0.1 )
1120  {
1121  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1122  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1123  dNdxP /= modul2;
1124  }
1125  return dNdxP;
1126 
1127 } // end of PAIdNdxPlasmon
1128 
1130 //
1131 // Calculation od dN/dx of collisions with creation of longitudinal EM
1132 // resonance excitations (plasmons, delta-electrons)
1133 
1134 G4double G4PAIxSection::PAIdNdxResonance( G4int i ,
1135  G4double betaGammaSq )
1136 {
1137  G4double resonance, modul2, dNdxP;
1138  G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
1139 
1140  cofBetaBohr = 4.0;
1141  betaBohr2 = fine_structure_const*fine_structure_const;
1142  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1143 
1144  be2 = betaGammaSq/(1 + betaGammaSq);
1145  be4 = be2*be2;
1146 
1147  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1148  resonance *= fImPartDielectricConst[i]/hbarc;
1149 
1150 
1151  dNdxP = resonance;
1152 
1153  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1154 
1155  dNdxP *= fine_structure_const/be2/pi;
1156  dNdxP *= (1-exp(-be4/betaBohr4));
1157 
1158  if( fDensity >= 0.1 )
1159  {
1160  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1161  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1162  dNdxP /= modul2;
1163  }
1164  return dNdxP;
1165 
1166 } // end of PAIdNdxResonance
1167 
1169 //
1170 // Calculation of the PAI integral cross-section
1171 // fIntegralPAIxSection[1] = specific primary ionisation, 1/cm
1172 // and fIntegralPAIxSection[0] = mean energy loss per cm in keV/cm
1173 
1174 void G4PAIxSection::IntegralPAIxSection()
1175 {
1176  fIntegralPAIxSection[fSplineNumber] = 0;
1177  fIntegralPAIdEdx[fSplineNumber] = 0;
1178  fIntegralPAIxSection[0] = 0;
1179  G4int k = fIntervalNumber -1;
1180 
1181  for(G4int i = fSplineNumber-1; i >= 1; i--)
1182  {
1183  if(fSplineEnergy[i] >= fEnergyInterval[k])
1184  {
1185  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i);
1186  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
1187  }
1188  else
1189  {
1190  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] +
1191  SumOverBorder(i+1,fEnergyInterval[k]);
1192  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
1193  SumOverBorderdEdx(i+1,fEnergyInterval[k]);
1194  k--;
1195  }
1196  }
1197 } // end of IntegralPAIxSection
1198 
1200 //
1201 // Calculation of the PAI Cerenkov integral cross-section
1202 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
1203 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
1204 
1205 void G4PAIxSection::IntegralCerenkov()
1206 {
1207  G4int i, k;
1208  fIntegralCerenkov[fSplineNumber] = 0;
1209  fIntegralCerenkov[0] = 0;
1210  k = fIntervalNumber -1;
1211 
1212  for( i = fSplineNumber-1; i >= 1; i-- )
1213  {
1214  if(fSplineEnergy[i] >= fEnergyInterval[k])
1215  {
1216  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
1217  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1218  }
1219  else
1220  {
1221  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
1222  SumOverBordCerenkov(i+1,fEnergyInterval[k]);
1223  k--;
1224  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1225  }
1226  }
1227 
1228 } // end of IntegralCerenkov
1229 
1231 //
1232 // Calculation of the PAI MM-Cerenkov integral cross-section
1233 // fIntegralMM[1] = specific MM-Cerenkov ionisation, 1/cm
1234 // and fIntegralMM[0] = mean MM-Cerenkov loss per cm in keV/cm
1235 
1236 void G4PAIxSection::IntegralMM()
1237 {
1238  G4int i, k;
1239  fIntegralMM[fSplineNumber] = 0;
1240  fIntegralMM[0] = 0;
1241  k = fIntervalNumber -1;
1242 
1243  for( i = fSplineNumber-1; i >= 1; i-- )
1244  {
1245  if(fSplineEnergy[i] >= fEnergyInterval[k])
1246  {
1247  fIntegralMM[i] = fIntegralMM[i+1] + SumOverInterMM(i);
1248  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1249  }
1250  else
1251  {
1252  fIntegralMM[i] = fIntegralMM[i+1] +
1253  SumOverBordMM(i+1,fEnergyInterval[k]);
1254  k--;
1255  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1256  }
1257  }
1258 
1259 } // end of IntegralMM
1260 
1262 //
1263 // Calculation of the PAI Plasmon integral cross-section
1264 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
1265 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
1266 
1267 void G4PAIxSection::IntegralPlasmon()
1268 {
1269  fIntegralPlasmon[fSplineNumber] = 0;
1270  fIntegralPlasmon[0] = 0;
1271  G4int k = fIntervalNumber -1;
1272  for(G4int i=fSplineNumber-1;i>=1;i--)
1273  {
1274  if(fSplineEnergy[i] >= fEnergyInterval[k])
1275  {
1276  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
1277  }
1278  else
1279  {
1280  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
1281  SumOverBordPlasmon(i+1,fEnergyInterval[k]);
1282  k--;
1283  }
1284  }
1285 
1286 } // end of IntegralPlasmon
1287 
1289 //
1290 // Calculation of the PAI resonance integral cross-section
1291 // fIntegralResonance[1] = resonance primary ionisation, 1/cm
1292 // and fIntegralResonance[0] = mean resonance loss per cm in keV/cm
1293 
1294 void G4PAIxSection::IntegralResonance()
1295 {
1296  fIntegralResonance[fSplineNumber] = 0;
1297  fIntegralResonance[0] = 0;
1298  G4int k = fIntervalNumber -1;
1299  for(G4int i=fSplineNumber-1;i>=1;i--)
1300  {
1301  if(fSplineEnergy[i] >= fEnergyInterval[k])
1302  {
1303  fIntegralResonance[i] = fIntegralResonance[i+1] + SumOverInterResonance(i);
1304  }
1305  else
1306  {
1307  fIntegralResonance[i] = fIntegralResonance[i+1] +
1308  SumOverBordResonance(i+1,fEnergyInterval[k]);
1309  k--;
1310  }
1311  }
1312 
1313 } // end of IntegralResonance
1314 
1316 //
1317 // Calculation the PAI integral cross-section inside
1318 // of interval of continuous values of photo-ionisation
1319 // cross-section. Parameter 'i' is the number of interval.
1320 
1321 G4double G4PAIxSection::SumOverInterval( G4int i )
1322 {
1323  G4double x0,x1,y0,yy1,a,b,c,result;
1324 
1325  x0 = fSplineEnergy[i];
1326  x1 = fSplineEnergy[i+1];
1327  y0 = fDifPAIxSection[i];
1328  yy1 = fDifPAIxSection[i+1];
1329  c = x1/x0;
1330  a = log10(yy1/y0)/log10(c);
1331  // b = log10(y0) - a*log10(x0);
1332  b = y0/pow(x0,a);
1333  a += 1.;
1334  if( std::fabs(a) < 1.e-6 )
1335  {
1336  result = b*log(x1/x0);
1337  }
1338  else
1339  {
1340  result = y0*(x1*pow(c,a-1) - x0)/a;
1341  }
1342  a += 1.;
1343  if( std::fabs(a) < 1.e-6 )
1344  {
1345  fIntegralPAIxSection[0] += b*log(x1/x0);
1346  }
1347  else
1348  {
1349  fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1350  }
1351  return result;
1352 
1353 } // end of SumOverInterval
1354 
1356 
1357 G4double G4PAIxSection::SumOverIntervaldEdx( G4int i )
1358 {
1359  G4double x0,x1,y0,yy1,a,b,c,result;
1360 
1361  x0 = fSplineEnergy[i];
1362  x1 = fSplineEnergy[i+1];
1363  y0 = fDifPAIxSection[i];
1364  yy1 = fDifPAIxSection[i+1];
1365  c = x1/x0;
1366  a = log10(yy1/y0)/log10(c);
1367  // b = log10(y0) - a*log10(x0);
1368  b = y0/pow(x0,a);
1369  a += 2;
1370  if(a == 0)
1371  {
1372  result = b*log(x1/x0);
1373  }
1374  else
1375  {
1376  result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1377  }
1378  return result;
1379 
1380 } // end of SumOverInterval
1381 
1383 //
1384 // Calculation the PAI Cerenkov integral cross-section inside
1385 // of interval of continuous values of photo-ionisation Cerenkov
1386 // cross-section. Parameter 'i' is the number of interval.
1387 
1388 G4double G4PAIxSection::SumOverInterCerenkov( G4int i )
1389 {
1390  G4double x0,x1,y0,yy1,a,b,c,result;
1391 
1392  x0 = fSplineEnergy[i];
1393  x1 = fSplineEnergy[i+1];
1394  y0 = fdNdxCerenkov[i];
1395  yy1 = fdNdxCerenkov[i+1];
1396  // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1397  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1398 
1399  c = x1/x0;
1400  a = log10(yy1/y0)/log10(c);
1401  b = y0/pow(x0,a);
1402 
1403  a += 1.0;
1404  if(a == 0) result = b*log(c);
1405  else result = y0*(x1*pow(c,a-1) - x0)/a;
1406  a += 1.0;
1407 
1408  if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
1409  else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1410  // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1411  return result;
1412 
1413 } // end of SumOverInterCerenkov
1414 
1416 //
1417 // Calculation the PAI MM-Cerenkov integral cross-section inside
1418 // of interval of continuous values of photo-ionisation Cerenkov
1419 // cross-section. Parameter 'i' is the number of interval.
1420 
1421 G4double G4PAIxSection::SumOverInterMM( G4int i )
1422 {
1423  G4double x0,x1,y0,yy1,a,b,c,result;
1424 
1425  x0 = fSplineEnergy[i];
1426  x1 = fSplineEnergy[i+1];
1427  y0 = fdNdxMM[i];
1428  yy1 = fdNdxMM[i+1];
1429  // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1430  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1431 
1432  c = x1/x0;
1433  a = log10(yy1/y0)/log10(c);
1434  b = y0/pow(x0,a);
1435 
1436  a += 1.0;
1437  if(a == 0) result = b*log(c);
1438  else result = y0*(x1*pow(c,a-1) - x0)/a;
1439  a += 1.0;
1440 
1441  if( a == 0 ) fIntegralMM[0] += b*log(x1/x0);
1442  else fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1443  // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1444  return result;
1445 
1446 } // end of SumOverInterMM
1447 
1449 //
1450 // Calculation the PAI Plasmon integral cross-section inside
1451 // of interval of continuous values of photo-ionisation Plasmon
1452 // cross-section. Parameter 'i' is the number of interval.
1453 
1454 G4double G4PAIxSection::SumOverInterPlasmon( G4int i )
1455 {
1456  G4double x0,x1,y0,yy1,a,b,c,result;
1457 
1458  x0 = fSplineEnergy[i];
1459  x1 = fSplineEnergy[i+1];
1460  y0 = fdNdxPlasmon[i];
1461  yy1 = fdNdxPlasmon[i+1];
1462  c =x1/x0;
1463  a = log10(yy1/y0)/log10(c);
1464  // b = log10(y0) - a*log10(x0);
1465  b = y0/pow(x0,a);
1466 
1467  a += 1.0;
1468  if(a == 0) result = b*log(x1/x0);
1469  else result = y0*(x1*pow(c,a-1) - x0)/a;
1470  a += 1.0;
1471 
1472  if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
1473  else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1474 
1475  return result;
1476 
1477 } // end of SumOverInterPlasmon
1478 
1480 //
1481 // Calculation the PAI resonance integral cross-section inside
1482 // of interval of continuous values of photo-ionisation resonance
1483 // cross-section. Parameter 'i' is the number of interval.
1484 
1485 G4double G4PAIxSection::SumOverInterResonance( G4int i )
1486 {
1487  G4double x0,x1,y0,yy1,a,b,c,result;
1488 
1489  x0 = fSplineEnergy[i];
1490  x1 = fSplineEnergy[i+1];
1491  y0 = fdNdxResonance[i];
1492  yy1 = fdNdxResonance[i+1];
1493  c =x1/x0;
1494  a = log10(yy1/y0)/log10(c);
1495  // b = log10(y0) - a*log10(x0);
1496  b = y0/pow(x0,a);
1497 
1498  a += 1.0;
1499  if(a == 0) result = b*log(x1/x0);
1500  else result = y0*(x1*pow(c,a-1) - x0)/a;
1501  a += 1.0;
1502 
1503  if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0);
1504  else fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1505 
1506  return result;
1507 
1508 } // end of SumOverInterResonance
1509 
1511 //
1512 // Integration of PAI cross-section for the case of
1513 // passing across border between intervals
1514 
1515 G4double G4PAIxSection::SumOverBorder( G4int i ,
1516  G4double en0 )
1517 {
1518  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1519 
1520  e0 = en0;
1521  x0 = fSplineEnergy[i];
1522  x1 = fSplineEnergy[i+1];
1523  y0 = fDifPAIxSection[i];
1524  yy1 = fDifPAIxSection[i+1];
1525 
1526  //c = x1/x0;
1527  d = e0/x0;
1528  a = log10(yy1/y0)/log10(x1/x0);
1529  // b0 = log10(y0) - a*log10(x0);
1530  b = y0/pow(x0,a); // pow(10.,b);
1531 
1532  a += 1.;
1533  if( std::fabs(a) < 1.e-6 )
1534  {
1535  result = b*log(x0/e0);
1536  }
1537  else
1538  {
1539  result = y0*(x0 - e0*pow(d,a-1))/a;
1540  }
1541  a += 1.;
1542  if( std::fabs(a) < 1.e-6 )
1543  {
1544  fIntegralPAIxSection[0] += b*log(x0/e0);
1545  }
1546  else
1547  {
1548  fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1549  }
1550  x0 = fSplineEnergy[i - 1];
1551  x1 = fSplineEnergy[i - 2];
1552  y0 = fDifPAIxSection[i - 1];
1553  yy1 = fDifPAIxSection[i - 2];
1554 
1555  //c = x1/x0;
1556  d = e0/x0;
1557  a = log10(yy1/y0)/log10(x1/x0);
1558  // b0 = log10(y0) - a*log10(x0);
1559  b = y0/pow(x0,a);
1560  a += 1.;
1561  if( std::fabs(a) < 1.e-6 )
1562  {
1563  result += b*log(e0/x0);
1564  }
1565  else
1566  {
1567  result += y0*(e0*pow(d,a-1) - x0)/a;
1568  }
1569  a += 1.;
1570  if( std::fabs(a) < 1.e-6 )
1571  {
1572  fIntegralPAIxSection[0] += b*log(e0/x0);
1573  }
1574  else
1575  {
1576  fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1577  }
1578  return result;
1579 
1580 }
1581 
1583 
1584 G4double G4PAIxSection::SumOverBorderdEdx( G4int i ,
1585  G4double en0 )
1586 {
1587  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1588 
1589  e0 = en0;
1590  x0 = fSplineEnergy[i];
1591  x1 = fSplineEnergy[i+1];
1592  y0 = fDifPAIxSection[i];
1593  yy1 = fDifPAIxSection[i+1];
1594 
1595  //c = x1/x0;
1596  d = e0/x0;
1597  a = log10(yy1/y0)/log10(x1/x0);
1598  // b0 = log10(y0) - a*log10(x0);
1599  b = y0/pow(x0,a); // pow(10.,b);
1600 
1601  a += 2;
1602  if(a == 0)
1603  {
1604  result = b*log(x0/e0);
1605  }
1606  else
1607  {
1608  result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1609  }
1610  x0 = fSplineEnergy[i - 1];
1611  x1 = fSplineEnergy[i - 2];
1612  y0 = fDifPAIxSection[i - 1];
1613  yy1 = fDifPAIxSection[i - 2];
1614 
1615  // c = x1/x0;
1616  d = e0/x0;
1617  a = log10(yy1/y0)/log10(x1/x0);
1618  // b0 = log10(y0) - a*log10(x0);
1619  b = y0/pow(x0,a);
1620  a += 2;
1621  if(a == 0)
1622  {
1623  result += b*log(e0/x0);
1624  }
1625  else
1626  {
1627  result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1628  }
1629  return result;
1630 
1631 }
1632 
1634 //
1635 // Integration of Cerenkov cross-section for the case of
1636 // passing across border between intervals
1637 
1638 G4double G4PAIxSection::SumOverBordCerenkov( G4int i ,
1639  G4double en0 )
1640 {
1641  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1642 
1643  e0 = en0;
1644  x0 = fSplineEnergy[i];
1645  x1 = fSplineEnergy[i+1];
1646  y0 = fdNdxCerenkov[i];
1647  yy1 = fdNdxCerenkov[i+1];
1648 
1649  // G4cout<<G4endl;
1650  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1651  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1652  c = x1/x0;
1653  d = e0/x0;
1654  a = log10(yy1/y0)/log10(c);
1655  // b0 = log10(y0) - a*log10(x0);
1656  b = y0/pow(x0,a); // pow(10.,b0);
1657 
1658  a += 1.0;
1659  if( a == 0 ) result = b*log(x0/e0);
1660  else result = y0*(x0 - e0*pow(d,a-1))/a;
1661  a += 1.0;
1662 
1663  if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
1664  else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1665 
1666 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
1667 
1668  x0 = fSplineEnergy[i - 1];
1669  x1 = fSplineEnergy[i - 2];
1670  y0 = fdNdxCerenkov[i - 1];
1671  yy1 = fdNdxCerenkov[i - 2];
1672 
1673  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1674  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1675 
1676  c = x1/x0;
1677  d = e0/x0;
1678  a = log10(yy1/y0)/log10(x1/x0);
1679  // b0 = log10(y0) - a*log10(x0);
1680  b = y0/pow(x0,a); // pow(10.,b0);
1681 
1682  a += 1.0;
1683  if( a == 0 ) result += b*log(e0/x0);
1684  else result += y0*(e0*pow(d,a-1) - x0 )/a;
1685  a += 1.0;
1686 
1687  if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
1688  else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1689 
1690  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
1691  // <<b<<"; result = "<<result<<G4endl;
1692 
1693  return result;
1694 
1695 }
1696 
1698 //
1699 // Integration of MM-Cerenkov cross-section for the case of
1700 // passing across border between intervals
1701 
1702 G4double G4PAIxSection::SumOverBordMM( G4int i ,
1703  G4double en0 )
1704 {
1705  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1706 
1707  e0 = en0;
1708  x0 = fSplineEnergy[i];
1709  x1 = fSplineEnergy[i+1];
1710  y0 = fdNdxMM[i];
1711  yy1 = fdNdxMM[i+1];
1712 
1713  // G4cout<<G4endl;
1714  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1715  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1716  c = x1/x0;
1717  d = e0/x0;
1718  a = log10(yy1/y0)/log10(c);
1719  // b0 = log10(y0) - a*log10(x0);
1720  b = y0/pow(x0,a); // pow(10.,b0);
1721 
1722  a += 1.0;
1723  if( a == 0 ) result = b*log(x0/e0);
1724  else result = y0*(x0 - e0*pow(d,a-1))/a;
1725  a += 1.0;
1726 
1727  if( a == 0 ) fIntegralMM[0] += b*log(x0/e0);
1728  else fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1729 
1730 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
1731 
1732  x0 = fSplineEnergy[i - 1];
1733  x1 = fSplineEnergy[i - 2];
1734  y0 = fdNdxMM[i - 1];
1735  yy1 = fdNdxMM[i - 2];
1736 
1737  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1738  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1739 
1740  c = x1/x0;
1741  d = e0/x0;
1742  a = log10(yy1/y0)/log10(x1/x0);
1743  // b0 = log10(y0) - a*log10(x0);
1744  b = y0/pow(x0,a); // pow(10.,b0);
1745 
1746  a += 1.0;
1747  if( a == 0 ) result += b*log(e0/x0);
1748  else result += y0*(e0*pow(d,a-1) - x0 )/a;
1749  a += 1.0;
1750 
1751  if( a == 0 ) fIntegralMM[0] += b*log(e0/x0);
1752  else fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1753 
1754  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
1755  // <<b<<"; result = "<<result<<G4endl;
1756 
1757  return result;
1758 
1759 }
1760 
1762 //
1763 // Integration of Plasmon cross-section for the case of
1764 // passing across border between intervals
1765 
1766 G4double G4PAIxSection::SumOverBordPlasmon( G4int i ,
1767  G4double en0 )
1768 {
1769  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
1770 
1771  e0 = en0;
1772  x0 = fSplineEnergy[i];
1773  x1 = fSplineEnergy[i+1];
1774  y0 = fdNdxPlasmon[i];
1775  yy1 = fdNdxPlasmon[i+1];
1776 
1777  c = x1/x0;
1778  d = e0/x0;
1779  a = log10(yy1/y0)/log10(c);
1780  // b0 = log10(y0) - a*log10(x0);
1781  b = y0/pow(x0,a); //pow(10.,b);
1782 
1783  a += 1.0;
1784  if( a == 0 ) result = b*log(x0/e0);
1785  else result = y0*(x0 - e0*pow(d,a-1))/a;
1786  a += 1.0;
1787 
1788  if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
1789  else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1790 
1791  x0 = fSplineEnergy[i - 1];
1792  x1 = fSplineEnergy[i - 2];
1793  y0 = fdNdxPlasmon[i - 1];
1794  yy1 = fdNdxPlasmon[i - 2];
1795 
1796  c = x1/x0;
1797  d = e0/x0;
1798  a = log10(yy1/y0)/log10(c);
1799  // b0 = log10(y0) - a*log10(x0);
1800  b = y0/pow(x0,a);// pow(10.,b0);
1801 
1802  a += 1.0;
1803  if( a == 0 ) result += b*log(e0/x0);
1804  else result += y0*(e0*pow(d,a-1) - x0)/a;
1805  a += 1.0;
1806 
1807  if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
1808  else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1809 
1810  return result;
1811 
1812 }
1813 
1815 //
1816 // Integration of resonance cross-section for the case of
1817 // passing across border between intervals
1818 
1819 G4double G4PAIxSection::SumOverBordResonance( G4int i ,
1820  G4double en0 )
1821 {
1822  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
1823 
1824  e0 = en0;
1825  x0 = fSplineEnergy[i];
1826  x1 = fSplineEnergy[i+1];
1827  y0 = fdNdxResonance[i];
1828  yy1 = fdNdxResonance[i+1];
1829 
1830  c = x1/x0;
1831  d = e0/x0;
1832  a = log10(yy1/y0)/log10(c);
1833  // b0 = log10(y0) - a*log10(x0);
1834  b = y0/pow(x0,a); //pow(10.,b);
1835 
1836  a += 1.0;
1837  if( a == 0 ) result = b*log(x0/e0);
1838  else result = y0*(x0 - e0*pow(d,a-1))/a;
1839  a += 1.0;
1840 
1841  if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0);
1842  else fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1843 
1844  x0 = fSplineEnergy[i - 1];
1845  x1 = fSplineEnergy[i - 2];
1846  y0 = fdNdxResonance[i - 1];
1847  yy1 = fdNdxResonance[i - 2];
1848 
1849  c = x1/x0;
1850  d = e0/x0;
1851  a = log10(yy1/y0)/log10(c);
1852  // b0 = log10(y0) - a*log10(x0);
1853  b = y0/pow(x0,a);// pow(10.,b0);
1854 
1855  a += 1.0;
1856  if( a == 0 ) result += b*log(e0/x0);
1857  else result += y0*(e0*pow(d,a-1) - x0)/a;
1858  a += 1.0;
1859 
1860  if( a == 0 ) fIntegralResonance[0] += b*log(e0/x0);
1861  else fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1862 
1863  return result;
1864 
1865 }
1866 
1868 //
1869 // Returns random PAI-total energy loss over step
1870 
1871 G4double G4PAIxSection::GetStepEnergyLoss( G4double step )
1872 {
1873  G4long numOfCollisions;
1874  G4double meanNumber, loss = 0.0;
1875 
1876  // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl;
1877 
1878  meanNumber = fIntegralPAIxSection[1]*step;
1879  numOfCollisions = G4Poisson(meanNumber);
1880 
1881  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1882 
1883  while(numOfCollisions)
1884  {
1885  loss += GetEnergyTransfer();
1886  numOfCollisions--;
1887  }
1888  // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
1889 
1890  return loss;
1891 }
1892 
1894 //
1895 // Returns random PAI-total energy transfer in one collision
1896 
1897 G4double G4PAIxSection::GetEnergyTransfer()
1898 {
1899  G4int iTransfer ;
1900 
1901  G4double energyTransfer, position;
1902 
1903  position = fIntegralPAIxSection[1]*G4UniformRand();
1904 
1905  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
1906  {
1907  if( position >= fIntegralPAIxSection[iTransfer] ) break;
1908  }
1909  if(iTransfer > fSplineNumber) iTransfer--;
1910 
1911  energyTransfer = fSplineEnergy[iTransfer];
1912 
1913  if(iTransfer > 1)
1914  {
1915  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
1916  }
1917  return energyTransfer;
1918 }
1919 
1921 //
1922 // Returns random Cerenkov energy loss over step
1923 
1924 G4double G4PAIxSection::GetStepCerenkovLoss( G4double step )
1925 {
1926  G4long numOfCollisions;
1927  G4double meanNumber, loss = 0.0;
1928 
1929  // G4cout<<" G4PAIxSection::GetStepCerenkovLoss "<<G4endl;
1930 
1931  meanNumber = fIntegralCerenkov[1]*step;
1932  numOfCollisions = G4Poisson(meanNumber);
1933 
1934  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1935 
1936  while(numOfCollisions)
1937  {
1938  loss += GetCerenkovEnergyTransfer();
1939  numOfCollisions--;
1940  }
1941  // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
1942 
1943  return loss;
1944 }
1945 
1947 //
1948 // Returns random MM-Cerenkov energy loss over step
1949 
1950 G4double G4PAIxSection::GetStepMMLoss( G4double step )
1951 {
1952  G4long numOfCollisions;
1953  G4double meanNumber, loss = 0.0;
1954 
1955  // G4cout<<" G4PAIxSection::GetStepMMLoss "<<G4endl;
1956 
1957  meanNumber = fIntegralMM[1]*step;
1958  numOfCollisions = G4Poisson(meanNumber);
1959 
1960  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1961 
1962  while(numOfCollisions)
1963  {
1964  loss += GetMMEnergyTransfer();
1965  numOfCollisions--;
1966  }
1967  // G4cout<<"PAI MM-Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
1968 
1969  return loss;
1970 }
1971 
1973 //
1974 // Returns Cerenkov energy transfer in one collision
1975 
1976 G4double G4PAIxSection::GetCerenkovEnergyTransfer()
1977 {
1978  G4int iTransfer ;
1979 
1980  G4double energyTransfer, position;
1981 
1982  position = fIntegralCerenkov[1]*G4UniformRand();
1983 
1984  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
1985  {
1986  if( position >= fIntegralCerenkov[iTransfer] ) break;
1987  }
1988  if(iTransfer > fSplineNumber) iTransfer--;
1989 
1990  energyTransfer = fSplineEnergy[iTransfer];
1991 
1992  if(iTransfer > 1)
1993  {
1994  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
1995  }
1996  return energyTransfer;
1997 }
1998 
2000 //
2001 // Returns MM-Cerenkov energy transfer in one collision
2002 
2003 G4double G4PAIxSection::GetMMEnergyTransfer()
2004 {
2005  G4int iTransfer ;
2006 
2007  G4double energyTransfer, position;
2008 
2009  position = fIntegralMM[1]*G4UniformRand();
2010 
2011  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2012  {
2013  if( position >= fIntegralMM[iTransfer] ) break;
2014  }
2015  if(iTransfer > fSplineNumber) iTransfer--;
2016 
2017  energyTransfer = fSplineEnergy[iTransfer];
2018 
2019  if(iTransfer > 1)
2020  {
2021  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2022  }
2023  return energyTransfer;
2024 }
2025 
2027 //
2028 // Returns random plasmon energy loss over step
2029 
2030 G4double G4PAIxSection::GetStepPlasmonLoss( G4double step )
2031 {
2032  G4long numOfCollisions;
2033  G4double meanNumber, loss = 0.0;
2034 
2035  // G4cout<<" G4PAIxSection::GetStepPlasmonLoss "<<G4endl;
2036 
2037  meanNumber = fIntegralPlasmon[1]*step;
2038  numOfCollisions = G4Poisson(meanNumber);
2039 
2040  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2041 
2042  while(numOfCollisions)
2043  {
2044  loss += GetPlasmonEnergyTransfer();
2045  numOfCollisions--;
2046  }
2047  // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
2048 
2049  return loss;
2050 }
2051 
2053 //
2054 // Returns plasmon energy transfer in one collision
2055 
2056 G4double G4PAIxSection::GetPlasmonEnergyTransfer()
2057 {
2058  G4int iTransfer ;
2059 
2060  G4double energyTransfer, position;
2061 
2062  position = fIntegralPlasmon[1]*G4UniformRand();
2063 
2064  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2065  {
2066  if( position >= fIntegralPlasmon[iTransfer] ) break;
2067  }
2068  if(iTransfer > fSplineNumber) iTransfer--;
2069 
2070  energyTransfer = fSplineEnergy[iTransfer];
2071 
2072  if(iTransfer > 1)
2073  {
2074  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2075  }
2076  return energyTransfer;
2077 }
2078 
2080 //
2081 // Returns random resonance energy loss over step
2082 
2083 G4double G4PAIxSection::GetStepResonanceLoss( G4double step )
2084 {
2085  G4long numOfCollisions;
2086  G4double meanNumber, loss = 0.0;
2087 
2088  // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl;
2089 
2090  meanNumber = fIntegralResonance[1]*step;
2091  numOfCollisions = G4Poisson(meanNumber);
2092 
2093  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2094 
2095  while(numOfCollisions)
2096  {
2097  loss += GetResonanceEnergyTransfer();
2098  numOfCollisions--;
2099  }
2100  // G4cout<<"PAI resonance loss = "<<loss/keV<<" keV"<<G4endl;
2101 
2102  return loss;
2103 }
2104 
2105 
2107 //
2108 // Returns resonance energy transfer in one collision
2109 
2110 G4double G4PAIxSection::GetResonanceEnergyTransfer()
2111 {
2112  G4int iTransfer ;
2113 
2114  G4double energyTransfer, position;
2115 
2116  position = fIntegralResonance[1]*G4UniformRand();
2117 
2118  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2119  {
2120  if( position >= fIntegralResonance[iTransfer] ) break;
2121  }
2122  if(iTransfer > fSplineNumber) iTransfer--;
2123 
2124  energyTransfer = fSplineEnergy[iTransfer];
2125 
2126  if(iTransfer > 1)
2127  {
2128  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2129  }
2130  return energyTransfer;
2131 }
2132 
2133 
2135 //
2136 // Returns Rutherford energy transfer in one collision
2137 
2138 G4double G4PAIxSection::GetRutherfordEnergyTransfer()
2139 {
2140  G4int iTransfer ;
2141 
2142  G4double energyTransfer, position;
2143 
2144  position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand();
2145 
2146  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2147  {
2148  if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break;
2149  }
2150  if(iTransfer > fSplineNumber) iTransfer--;
2151 
2152  energyTransfer = fSplineEnergy[iTransfer];
2153 
2154  if(iTransfer > 1)
2155  {
2156  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2157  }
2158  return energyTransfer;
2159 }
2160 
2162 //
2163 
2164 void G4PAIxSection::CallError(G4int i, const G4String& methodName) const
2165 {
2166  G4String head = "G4PAIxSection::" + methodName + "()";
2167  G4ExceptionDescription ed;
2168  ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber << G4endl;
2169  G4Exception(head,"pai001",FatalException,ed);
2170 }
2171 
2173 //
2174 // Init array of Lorentz factors
2175 //
2176 
2177 G4int G4PAIxSection::fNumberOfGammas = 111;
2178 
2179 const G4double G4PAIxSection::fLorentzFactor[112] = // fNumberOfGammas+1
2180 {
2181 0.0,
2182 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
2183 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
2184 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
2185 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
2186 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
2187 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
2188 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
2189 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
2190 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
2191 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
2192 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
2193 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
2194 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
2195 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
2196 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
2197 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
2198 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
2199 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
2200 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
2201 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
2202 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
2203 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
2204 1.065799e+05
2205 };
2206 
2208 //
2209 // The number of gamma for creation of spline (near ion-min , G ~ 4 )
2210 //
2211 
2212 const
2213 G4int G4PAIxSection::fRefGammaNumber = 29;
2214 
2215 
2216 //
2217 // end of G4PAIxSection implementation file
2218 //
2220 
Generated on Sat May 25 2013 14:33:37 for Geant4 by   doxygen 1.8.4