source: trunk/source/processes/electromagnetic/standard/src/G4PAIxSection.cc @ 846

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 47.6 KB
Line 
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,v 1.21 2006/06/29 19:53:20 gunter Exp $
28// GEANT4 tag $Name:  $
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 "G4ios.hh"
55#include "G4Poisson.hh"
56#include "G4Material.hh"
57#include "G4MaterialCutsCouple.hh"
58#include "G4SandiaTable.hh"
59
60using namespace std;
61
62/* ******************************************************************
63
64// Init  array of Lorentz factors
65
66const G4double G4PAIxSection::fLorentzFactor[22] =
67{
68          0.0 ,     1.1 ,   1.2 ,   1.3 ,    1.5 ,    1.8 ,  2.0 ,
69          2.5 ,     3.0 ,   4.0 ,   7.0 ,   10.0 ,   20.0 , 40.0 ,
70         70.0 ,   100.0 , 300.0 , 600.0 , 1000.0 , 3000.0 ,
71      10000.0 , 50000.0
72} ;
73
74const G4int G4PAIxSection::
75fRefGammaNumber = 29 ;         // The number of gamma for creation of
76                               // spline (9)
77
78***************************************************************** */ 
79
80// Local class constants
81
82const G4double G4PAIxSection::fDelta = 0.005 ; // energy shift from interval border
83const G4double G4PAIxSection::fError = 0.005 ; // error in lin-log approximation
84
85const G4int G4PAIxSection::fMaxSplineSize = 500 ;  // Max size of output spline
86                                                    // arrays
87
88//////////////////////////////////////////////////////////////////
89//
90// Constructor
91//
92
93G4PAIxSection::G4PAIxSection(G4MaterialCutsCouple* matCC)
94{
95  fDensity       = matCC->GetMaterial()->GetDensity();
96  G4int matIndex = matCC->GetMaterial()->GetIndex();
97  fSandia        = new G4SandiaTable(matIndex);
98
99  G4int i, j; 
100  fMatSandiaMatrix = new G4OrderedTable();
101 
102  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
103  {
104     fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
105  }                     
106  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
107  {
108    (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
109
110    for(j = 1; j < 5 ; j++)
111    {
112      (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
113    }     
114  }                     
115
116
117
118}
119
120G4PAIxSection::G4PAIxSection(G4int materialIndex,
121                             G4double maxEnergyTransfer)
122{
123   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ;
124   G4int i, j ;   
125
126      fDensity                = (*theMaterialTable)[materialIndex]->GetDensity() ;
127      fElectronDensity        = (*theMaterialTable)[materialIndex]->
128                             GetElectronDensity() ;
129      fIntervalNumber         = (*theMaterialTable)[materialIndex]->
130                             GetSandiaTable()->GetMatNbOfIntervals() ;
131      fIntervalNumber--;     
132      // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl ;
133
134      fEnergyInterval = new G4double[fIntervalNumber+2] ;
135      fA1             = new G4double[fIntervalNumber+2] ;
136      fA2             = new G4double[fIntervalNumber+2] ;
137      fA3             = new G4double[fIntervalNumber+2] ;
138      fA4             = new G4double[fIntervalNumber+2] ;
139
140      for(i = 1 ; i <= fIntervalNumber ; i++ )
141      {
142         if(((*theMaterialTable)[materialIndex]->
143    GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
144              i > fIntervalNumber               )
145         {
146            fEnergyInterval[i] = maxEnergyTransfer ;
147            fIntervalNumber = i ;
148            break;
149         }
150         fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
151                              GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
152         fA1[i]             = (*theMaterialTable)[materialIndex]->
153                              GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
154         fA2[i]             = (*theMaterialTable)[materialIndex]->
155                              GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
156         fA3[i]             = (*theMaterialTable)[materialIndex]->
157                              GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
158         fA4[i]             = (*theMaterialTable)[materialIndex]->
159                              GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
160         // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
161         //                               <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
162      }   
163      if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
164      {
165         fIntervalNumber++;
166         fEnergyInterval[fIntervalNumber] = maxEnergyTransfer ;
167      }
168
169      // Now checking, if two borders are too close together
170
171      for(i=1;i<fIntervalNumber;i++)
172      {
173        if(fEnergyInterval[i+1]-fEnergyInterval[i] >
174           1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
175        {
176          continue ;
177        }
178        else
179        {
180          for(j=i;j<fIntervalNumber;j++)
181          {
182            fEnergyInterval[j] = fEnergyInterval[j+1] ;
183                        fA1[j] = fA1[j+1] ;
184                        fA2[j] = fA2[j+1] ;
185                        fA3[j] = fA3[j+1] ;
186                        fA4[j] = fA4[j+1] ;
187          }
188          fIntervalNumber-- ;
189          i-- ;
190        }
191      }
192
193
194      /* *********************************
195
196      fSplineEnergy          = new G4double[fMaxSplineSize] ;   
197      fRePartDielectricConst = new G4double[fMaxSplineSize] ;   
198      fImPartDielectricConst = new G4double[fMaxSplineSize] ;   
199      fIntegralTerm          = new G4double[fMaxSplineSize] ;   
200      fDifPAIxSection        = new G4double[fMaxSplineSize] ;   
201      fIntegralPAIxSection   = new G4double[fMaxSplineSize] ;   
202     
203      for(i=0;i<fMaxSplineSize;i++)
204      {
205         fSplineEnergy[i]          = 0.0 ;   
206         fRePartDielectricConst[i] = 0.0 ;   
207         fImPartDielectricConst[i] = 0.0 ;   
208         fIntegralTerm[i]          = 0.0 ;   
209         fDifPAIxSection[i]        = 0.0 ;   
210         fIntegralPAIxSection[i]   = 0.0 ;   
211      }
212      **************************************************  */   
213
214      InitPAI() ;  // create arrays allocated above
215     
216      delete[] fEnergyInterval ;
217      delete[] fA1 ;
218      delete[] fA2 ;
219      delete[] fA3 ;
220      delete[] fA4 ;   
221}
222
223////////////////////////////////////////////////////////////////////////
224//
225// Constructor with beta*gamma square value
226
227G4PAIxSection::G4PAIxSection( G4int materialIndex,
228                              G4double maxEnergyTransfer,
229                              G4double betaGammaSq,
230                              G4double** photoAbsCof, 
231                              G4int intNumber                   )
232{
233   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
234   G4int i, j ;   
235   
236      fDensity                = (*theMaterialTable)[materialIndex]->GetDensity();
237      fElectronDensity        = (*theMaterialTable)[materialIndex]->
238                             GetElectronDensity() ;
239
240      fIntervalNumber         = intNumber ;
241      fIntervalNumber--;
242      //   G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl ;
243 
244      fEnergyInterval = new G4double[fIntervalNumber+2] ;
245      fA1             = new G4double[fIntervalNumber+2] ;
246      fA2             = new G4double[fIntervalNumber+2] ;
247      fA3             = new G4double[fIntervalNumber+2] ;
248      fA4             = new G4double[fIntervalNumber+2] ;
249
250      for( i = 1 ; i <= fIntervalNumber ; i++ )
251      {
252         if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
253             i > fIntervalNumber )
254         {
255            fEnergyInterval[i] = maxEnergyTransfer ;
256            fIntervalNumber = i ;
257            break;
258         }
259         fEnergyInterval[i] = photoAbsCof[i-1][0] ;
260         fA1[i]             = photoAbsCof[i-1][1] ;
261         fA2[i]             = photoAbsCof[i-1][2] ;
262         fA3[i]             = photoAbsCof[i-1][3] ;
263         fA4[i]             = photoAbsCof[i-1][4] ;
264         // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
265         //        <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
266      }
267  // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;   
268      if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
269      {
270         fIntervalNumber++;
271         fEnergyInterval[fIntervalNumber] = maxEnergyTransfer ;
272      }
273      for(i=1;i<=fIntervalNumber;i++)
274      {
275        //  G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
276        //    <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
277      }
278      // Now checking, if two borders are too close together
279
280      for( i = 1 ; i < fIntervalNumber ; i++ )
281      {
282        if(fEnergyInterval[i+1]-fEnergyInterval[i] >
283           1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
284        {
285          continue ;
286        }
287        else
288        {
289          for(j=i;j<fIntervalNumber;j++)
290          {
291            fEnergyInterval[j] = fEnergyInterval[j+1] ;
292                        fA1[j] = fA1[j+1] ;
293                        fA2[j] = fA2[j+1] ;
294                        fA3[j] = fA3[j+1] ;
295                        fA4[j] = fA4[j+1] ;
296          }
297          fIntervalNumber-- ;
298          i-- ;
299        }
300      }
301
302      // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
303     
304      G4double   betaGammaSqRef = 
305      fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
306
307      NormShift(betaGammaSqRef) ;             
308      SplainPAI(betaGammaSqRef) ;
309     
310      // Preparation of integral PAI cross section for input betaGammaSq
311   
312      for(i = 1 ; i <= fSplineNumber ; i++)
313      {
314         fdNdxCerenkov[i]   = PAIdNdxCerenkov(i,betaGammaSq);
315         fdNdxPlasmon[i]    = PAIdNdxPlasmon(i,betaGammaSq);
316         fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
317         // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
318         //    <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
319      }
320      IntegralCerenkov() ;
321      IntegralPlasmon() ;
322      IntegralPAIxSection() ;
323     
324      delete[] fEnergyInterval ;
325      delete[] fA1 ;
326      delete[] fA2 ;
327      delete[] fA3 ;
328      delete[] fA4 ;   
329}
330
331////////////////////////////////////////////////////////////////////////
332//
333// Test Constructor with beta*gamma square value
334
335G4PAIxSection::G4PAIxSection( G4int materialIndex,
336                              G4double maxEnergyTransfer,
337                              G4double betaGammaSq          )
338{
339   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
340
341   G4int i, j, numberOfElements ;   
342   
343   fDensity         = (*theMaterialTable)[materialIndex]->GetDensity();
344   fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity() ;
345   numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements() ;
346
347   G4int* thisMaterialZ = new G4int[numberOfElements] ;
348   
349   for( i = 0 ; i < numberOfElements ; i++ )
350   {
351         thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
352                                      GetElement(i)->GetZ() ;
353   }
354   G4SandiaTable     thisMaterialSandiaTable(materialIndex) ;
355   fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
356                           (thisMaterialZ,numberOfElements) ;   
357   fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
358                           ( thisMaterialZ ,
359                      (*theMaterialTable)[materialIndex]->GetFractionVector() ,
360                             numberOfElements,fIntervalNumber) ;
361
362   fIntervalNumber--;
363
364      fEnergyInterval = new G4double[fIntervalNumber+2] ;
365      fA1             = new G4double[fIntervalNumber+2] ;
366      fA2             = new G4double[fIntervalNumber+2] ;
367      fA3             = new G4double[fIntervalNumber+2] ;
368      fA4             = new G4double[fIntervalNumber+2] ;
369
370      for(i=1;i<=fIntervalNumber;i++)
371      {
372  if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
373          i > fIntervalNumber)
374         {
375            fEnergyInterval[i] = maxEnergyTransfer ;
376            fIntervalNumber = i ;
377            break;
378         }
379   fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) ;
380   fA1[i]             = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity ;
381   fA2[i]             = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity ;
382   fA3[i]             = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity ;
383   fA4[i]             = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity ;
384
385      }   
386      if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
387      {
388         fIntervalNumber++;
389         fEnergyInterval[fIntervalNumber] = maxEnergyTransfer ;
390         fA1[fIntervalNumber] = fA1[fIntervalNumber-1] ;
391         fA2[fIntervalNumber] = fA2[fIntervalNumber-1] ;
392         fA3[fIntervalNumber] = fA3[fIntervalNumber-1] ;
393         fA4[fIntervalNumber] = fA4[fIntervalNumber-1] ;
394      }
395      for(i=1;i<=fIntervalNumber;i++)
396      {
397        // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
398        //   <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
399      }
400      // Now checking, if two borders are too close together
401
402      for(i=1;i<fIntervalNumber;i++)
403      {
404        if(fEnergyInterval[i+1]-fEnergyInterval[i] >
405           1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
406        {
407          continue ;
408        }
409        else
410        {
411          for(j=i;j<fIntervalNumber;j++)
412          {
413            fEnergyInterval[j] = fEnergyInterval[j+1] ;
414                        fA1[j] = fA1[j+1] ;
415                        fA2[j] = fA2[j+1] ;
416                        fA3[j] = fA3[j+1] ;
417                        fA4[j] = fA4[j+1] ;
418          }
419          fIntervalNumber-- ;
420          i-- ;
421        }
422      }
423
424      /* *********************************
425      fSplineEnergy          = new G4double[fMaxSplineSize] ;   
426      fRePartDielectricConst = new G4double[fMaxSplineSize] ;   
427      fImPartDielectricConst = new G4double[fMaxSplineSize] ;   
428      fIntegralTerm          = new G4double[fMaxSplineSize] ;   
429      fDifPAIxSection        = new G4double[fMaxSplineSize] ;   
430      fIntegralPAIxSection   = new G4double[fMaxSplineSize] ;   
431     
432      for(i=0;i<fMaxSplineSize;i++)
433      {
434         fSplineEnergy[i]          = 0.0 ;   
435         fRePartDielectricConst[i] = 0.0 ;   
436         fImPartDielectricConst[i] = 0.0 ;   
437         fIntegralTerm[i]          = 0.0 ;   
438         fDifPAIxSection[i]        = 0.0 ;   
439         fIntegralPAIxSection[i]   = 0.0 ;   
440      }
441      */ ////////////////////////
442
443      // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
444     
445      G4double   betaGammaSqRef = 
446      fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
447
448      NormShift(betaGammaSqRef) ;             
449      SplainPAI(betaGammaSqRef) ;
450     
451      // Preparation of integral PAI cross section for input betaGammaSq
452   
453      for(i = 1 ; i <= fSplineNumber ; i++)
454      {
455         fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
456         fdNdxCerenkov[i]   = PAIdNdxCerenkov(i,betaGammaSq);
457         fdNdxPlasmon[i]    = PAIdNdxPlasmon(i,betaGammaSq);
458      }
459      IntegralPAIxSection() ;
460      IntegralCerenkov() ;
461      IntegralPlasmon() ;
462     
463      //   delete[] fEnergyInterval ;
464      delete[] fA1 ;
465      delete[] fA2 ;
466      delete[] fA3 ;
467      delete[] fA4 ;   
468}
469
470
471////////////////////////////////////////////////////////////////////////////
472//
473// Destructor
474
475G4PAIxSection::~G4PAIxSection()
476{
477   /* ************************
478   delete[] fSplineEnergy          ;   
479   delete[] fRePartDielectricConst ;   
480   delete[] fImPartDielectricConst ;   
481   delete[] fIntegralTerm          ;   
482   delete[] fDifPAIxSection        ;   
483   delete[] fIntegralPAIxSection   ;
484   */ ////////////////////////
485}
486
487/////////////////////////////////////////////////////////////////////////
488//
489// General control function for class G4PAIxSection
490//
491
492void G4PAIxSection::InitPAI()
493{   
494   G4int i ;
495   G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
496                          fLorentzFactor[fRefGammaNumber] - 1;
497
498   // Preparation of integral PAI cross section for reference gamma
499   
500   NormShift(betaGammaSq) ;             
501   SplainPAI(betaGammaSq) ;
502
503   IntegralPAIxSection() ;
504   IntegralCerenkov() ;
505   IntegralPlasmon() ;
506
507   for(i = 0 ; i<=fSplineNumber ; i++)
508   {
509      fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i] ;
510      if(i != 0) 
511      {
512         fPAItable[i][0] = fSplineEnergy[i] ;
513      }
514   }
515   fPAItable[0][0] = fSplineNumber ;
516   
517   for(G4int j = 1 ; j < 112 ; j++)       // for other gammas
518   {
519      if( j == fRefGammaNumber ) continue ;
520     
521      betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1 ;
522     
523      for(i = 1 ; i <= fSplineNumber ; i++)
524      {
525         fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
526         fdNdxCerenkov[i]   = PAIdNdxCerenkov(i,betaGammaSq);
527         fdNdxPlasmon[i]    = PAIdNdxPlasmon(i,betaGammaSq);
528      }
529      IntegralPAIxSection() ;
530      IntegralCerenkov() ;
531      IntegralPlasmon() ;
532     
533      for(i = 0 ; i <= fSplineNumber ; i++)
534      {
535         fPAItable[i][j] = fIntegralPAIxSection[i] ;
536      }
537   } 
538
539} 
540
541///////////////////////////////////////////////////////////////////////
542//
543// Shifting from borders to intervals Creation of first energy points
544//
545
546void G4PAIxSection::NormShift(G4double betaGammaSq)
547{
548  G4int i, j ;
549
550  for( i = 1 ; i <= fIntervalNumber-1 ; i++ )
551  {
552    for( j = 1 ; j <= 2 ; j++ )
553    {
554      fSplineNumber = (i-1)*2 + j ;
555
556      if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i  ]*(1+fDelta);
557      else         fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta); 
558      //    G4cout<<"cn = "<<fSplineNumber<<"; "<<"energy = "
559      //  <<fSplineEnergy[fSplineNumber]<<G4endl;
560    }
561  }
562  fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
563
564  j = 1 ;
565
566  for(i=2;i<=fSplineNumber;i++)
567  {
568    if(fSplineEnergy[i]<fEnergyInterval[j+1])
569    {
570         fIntegralTerm[i] = fIntegralTerm[i-1] + 
571                            RutherfordIntegral(j,fSplineEnergy[i-1],
572                                                 fSplineEnergy[i]   ) ;
573    }
574    else
575    {
576       G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
577                                           fEnergyInterval[j+1]   ) ;
578         j++;
579         fIntegralTerm[i] = fIntegralTerm[i-1] + x + 
580                            RutherfordIntegral(j,fEnergyInterval[j],
581                                                 fSplineEnergy[i]    ) ;
582    }
583    // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t"<<fIntegralTerm[i]<<"\n"<<G4endl;
584  } 
585  fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
586  fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber] ;
587
588  // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl ;
589
590          // Calculation of PAI differrential cross-section (1/(keV*cm))
591          // in the energy points near borders of energy intervals
592
593   for(G4int k=1;k<=fIntervalNumber-1;k++)
594   {
595      for(j=1;j<=2;j++)
596      {
597         i = (k-1)*2 + j ;
598         fImPartDielectricConst[i] = fNormalizationCof*
599                                     ImPartDielectricConst(k,fSplineEnergy[i]);
600         fRePartDielectricConst[i] = fNormalizationCof*
601                                     RePartDielectricConst(fSplineEnergy[i]);
602         fIntegralTerm[i] *= fNormalizationCof;
603
604         fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
605         fdNdxCerenkov[i]   = PAIdNdxCerenkov(i,betaGammaSq);
606         fdNdxPlasmon[i]    = PAIdNdxPlasmon(i,betaGammaSq);
607      }
608   }
609
610}  // end of NormShift
611
612/////////////////////////////////////////////////////////////////////////
613//
614// Creation of new energy points as geometrical mean of existing
615// one, calculation PAI_cs for them, while the error of logarithmic
616// linear approximation would be smaller than 'fError'
617
618void
619   G4PAIxSection::SplainPAI(G4double betaGammaSq)
620{
621   G4int k = 1 ;
622   G4int i = 1 ;
623
624   while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
625   {
626      if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
627      {
628          k++ ;   // Here next energy point is in next energy interval
629          i++;
630          continue;
631      }
632                       // Shifting of arrayes for inserting the geometrical
633                       // average of 'i' and 'i+1' energy points to 'i+1' place
634      fSplineNumber++;
635
636      for(G4int j = fSplineNumber; j >= i+2 ; j-- )
637      {
638         fSplineEnergy[j]          = fSplineEnergy[j-1];
639         fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
640         fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
641         fIntegralTerm[j]          = fIntegralTerm[j-1];
642
643         fDifPAIxSection[j] = fDifPAIxSection[j-1];
644         fdNdxCerenkov[j]   = fdNdxCerenkov[j-1];
645         fdNdxPlasmon[j]    = fdNdxPlasmon[j-1];
646      }
647      G4double x1  = fSplineEnergy[i];
648      G4double x2  = fSplineEnergy[i+1];
649      G4double yy1 = fDifPAIxSection[i];
650      G4double y2  = fDifPAIxSection[i+1];
651
652      G4double en1 = sqrt(x1*x2);
653      fSplineEnergy[i+1] = en1;
654
655                 // Calculation of logarithmic linear approximation
656                 // in this (enr) energy point, which number is 'i+1' now
657
658      G4double a = log10(y2/yy1)/log10(x2/x1);
659      G4double b = log10(yy1) - a*log10(x1);
660      G4double y = a*log10(en1) + b ;
661      y = pow(10.,y);
662
663                 // Calculation of the PAI dif. cross-section at this point
664
665      fImPartDielectricConst[i+1] = fNormalizationCof*
666                                    ImPartDielectricConst(k,fSplineEnergy[i+1]);
667      fRePartDielectricConst[i+1] = fNormalizationCof*
668                                    RePartDielectricConst(fSplineEnergy[i+1]);
669      fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
670                           RutherfordIntegral(k,fSplineEnergy[i],
671                                                fSplineEnergy[i+1]);
672
673      fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
674      fdNdxCerenkov[i+1]   = PAIdNdxCerenkov(i+1,betaGammaSq);
675      fdNdxPlasmon[i+1]    = PAIdNdxPlasmon(i+1,betaGammaSq);
676
677                  // Condition for next division of this segment or to pass
678                  // to higher energies
679
680      G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
681
682      if( x < 0 ) 
683      {
684         x = -x ;
685      }
686      if( x > fError && fSplineNumber < fMaxSplineSize-1 )
687      {
688         continue;  // next division
689      }
690      i += 2;  // pass to next segment
691
692   }   // close 'while'
693
694}  // end of SplainPAI
695
696
697////////////////////////////////////////////////////////////////////
698//
699// Integration over electrons that could be considered
700// quasi-free at energy transfer of interest
701
702G4double G4PAIxSection::RutherfordIntegral( G4int k,
703                                            G4double x1,
704                                            G4double x2   )
705{
706   G4double  c1, c2, c3 ;
707   // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;   
708   c1 = (x2 - x1)/x1/x2 ;
709   c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
710   c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
711   // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;   
712   
713   return  fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3 ;
714
715}   // end of RutherfordIntegral
716
717
718/////////////////////////////////////////////////////////////////
719//
720// Imaginary part of dielectric constant
721// (G4int k - interval number, G4double en1 - energy point)
722
723G4double G4PAIxSection::ImPartDielectricConst( G4int    k ,
724                                               G4double energy1 )
725{
726   G4double energy2,energy3,energy4,result;
727
728   energy2 = energy1*energy1;
729   energy3 = energy2*energy1;
730   energy4 = energy3*energy1;
731   
732   result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4 ; 
733   result *=hbarc/energy1 ;
734   
735   return result ;
736
737}  // end of ImPartDielectricConst
738
739
740//////////////////////////////////////////////////////////////////////////////
741//
742// Real part of dielectric constant minus unit: epsilon_1 - 1
743// (G4double enb - energy point)
744//
745
746G4double G4PAIxSection::RePartDielectricConst(G4double enb)
747{       
748   G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
749            c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
750
751   x0 = enb ;
752   result = 0 ;
753   
754   for(G4int i=1;i<=fIntervalNumber-1;i++)
755   {
756      x1 = fEnergyInterval[i] ;
757      x2 = fEnergyInterval[i+1] ;
758      xx1 = x1 - x0 ;
759      xx2 = x2 - x0 ;
760      xx12 = xx2/xx1 ;
761     
762      if(xx12<0)
763      {
764         xx12 = -xx12;
765      }
766      xln1 = log(x2/x1) ;
767      xln2 = log(xx12) ;
768      xln3 = log((x2 + x0)/(x1 + x0)) ;
769      x02 = x0*x0 ;
770      x03 = x02*x0 ;
771      x04 = x03*x0 ;
772      x05 = x04*x0;
773      c1  = (x2 - x1)/x1/x2 ;
774      c2  = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
775      c3  = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
776
777      result -= (fA1[i]/x02 + fA3[i]/x04)*xln1 ;
778      result -= (fA2[i]/x02 + fA4[i]/x04)*c1 ;
779      result -= fA3[i]*c2/2/x02 ;
780      result -= fA4[i]*c3/3/x02 ;
781
782      cof1 = fA1[i]/x02 + fA3[i]/x04 ;
783      cof2 = fA2[i]/x03 + fA4[i]/x05 ;
784
785      result += 0.5*(cof1 +cof2)*xln2 ;
786      result += 0.5*(cof1 - cof2)*xln3 ;
787   } 
788   result *= 2*hbarc/pi ;
789   
790   return result ;
791
792}   // end of RePartDielectricConst
793
794//////////////////////////////////////////////////////////////////////
795//
796// PAI differential cross-section in terms of
797// simplified Allison's equation
798//
799
800G4double G4PAIxSection::DifPAIxSection( G4int              i ,
801                                        G4double betaGammaSq  )
802{       
803   G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
804   //G4double beta, be4 ;
805   G4double be4 ;
806   G4double betaBohr2 = fine_structure_const*fine_structure_const ;
807   G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
808   be2 = betaGammaSq/(1 + betaGammaSq) ;
809   be4 = be2*be2 ;
810   //  beta = sqrt(be2) ;
811   cof = 1 ;
812   x1 = log(2*electron_mass_c2/fSplineEnergy[i]) ;
813
814   if( betaGammaSq < 0.01 ) x2 = log(be2) ;
815   else
816   {
817     x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
818                (1/betaGammaSq - fRePartDielectricConst[i]) + 
819                fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2 ;
820   }
821   if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
822   {
823     x6=0 ;
824   }
825   else
826   {
827     x3 = -fRePartDielectricConst[i] + 1/betaGammaSq ;
828     x5 = -1 - fRePartDielectricConst[i] +
829          be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
830          fImPartDielectricConst[i]*fImPartDielectricConst[i]) ;
831
832     x7 = atan2(fImPartDielectricConst[i],x3) ;
833     x6 = x5 * x7 ;
834   }
835    // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
836   
837   x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc ;
838   //   if( x4 < 0.0 ) x4 = 0.0 ;
839   x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 
840        fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
841
842   result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]) ;
843   if(result < 1.0e-8) result = 1.0e-8 ;
844   result *= fine_structure_const/be2/pi ;
845   //   result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
846   //  result *= (1-exp(-be2/betaBohr2)) ;
847   result *= (1-exp(-be4/betaBohr4)) ;
848   if(fDensity >= 0.1)
849   { 
850      result /= x8 ;
851   }
852   return result ;
853
854} // end of DifPAIxSection
855
856//////////////////////////////////////////////////////////////////////////
857//
858// Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
859
860G4double G4PAIxSection::PAIdNdxCerenkov( G4int    i ,
861                                         G4double betaGammaSq  )
862{       
863   G4double cof, logarithm, x3, x5, argument, modul2, dNdxC ; 
864   G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr ;
865
866   cof         = 1.0 ;
867   cofBetaBohr = 4.0 ;
868   betaBohr2   = fine_structure_const*fine_structure_const ;
869   betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr ;
870
871   be2 = betaGammaSq/(1 + betaGammaSq) ;
872   be4 = be2*be2 ;
873
874   if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
875   else
876   {
877     logarithm  = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
878                        (1/betaGammaSq - fRePartDielectricConst[i]) + 
879                        fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5 ;
880     logarithm += log(1+1.0/betaGammaSq) ;
881   }
882
883   if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
884   {
885     argument = 0.0 ;
886   }
887   else
888   {
889     x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq ;
890     x5 = -1.0 - fRePartDielectricConst[i] +
891          be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
892          fImPartDielectricConst[i]*fImPartDielectricConst[i]) ;
893     if( x3 == 0.0 ) argument = 0.5*pi;
894     else            argument = atan2(fImPartDielectricConst[i],x3) ;
895     argument *= x5  ;
896   }   
897   dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc ;
898 
899   if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
900
901   dNdxC *= fine_structure_const/be2/pi ;
902
903   dNdxC *= (1-exp(-be4/betaBohr4)) ;
904
905   if(fDensity >= 0.1)
906   { 
907      modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) + 
908                    fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
909      dNdxC /= modul2 ;
910   }
911   return dNdxC ;
912
913} // end of PAIdNdxCerenkov
914
915//////////////////////////////////////////////////////////////////////////
916//
917// Calculation od dN/dx of collisions with creation of longitudinal EM
918// excitations (plasmons, delta-electrons)
919
920G4double G4PAIxSection::PAIdNdxPlasmon( G4int    i ,
921                                        G4double betaGammaSq  )
922{       
923   G4double cof, resonance, modul2, dNdxP ;
924   G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr ;
925
926   cof = 1 ;
927   cofBetaBohr = 4.0 ;
928   betaBohr2   = fine_structure_const*fine_structure_const ;
929   betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr ;
930
931   be2 = betaGammaSq/(1 + betaGammaSq) ;
932   be4 = be2*be2 ;
933 
934   resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]) ; 
935   resonance *= fImPartDielectricConst[i]/hbarc ;
936
937
938   dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] ) ;
939
940   if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
941
942   dNdxP *= fine_structure_const/be2/pi ;
943   dNdxP *= (1-exp(-be4/betaBohr4)) ;
944
945   if( fDensity >= 0.1 )
946   { 
947     modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 
948        fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
949     dNdxP /= modul2 ;
950   }
951   return dNdxP ;
952
953} // end of PAIdNdxPlasmon
954
955////////////////////////////////////////////////////////////////////////
956//
957// Calculation of the PAI integral cross-section
958// fIntegralPAIxSection[1] = specific primary ionisation, 1/cm
959// and fIntegralPAIxSection[0] = mean energy loss per cm  in keV/cm
960
961void G4PAIxSection::IntegralPAIxSection()
962{
963  fIntegralPAIxSection[fSplineNumber] = 0 ;
964  fIntegralPAIdEdx[fSplineNumber]     = 0 ;
965  fIntegralPAIxSection[0]             = 0 ;
966  G4int k = fIntervalNumber -1 ;
967
968  for(G4int i = fSplineNumber-1 ; i >= 1 ; i--)
969  {
970    if(fSplineEnergy[i] >= fEnergyInterval[k])
971    {
972      fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i) ;
973      fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i) ;
974    }
975    else
976    {
977      fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + 
978                                   SumOverBorder(i+1,fEnergyInterval[k]) ;
979      fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + 
980                                   SumOverBorderdEdx(i+1,fEnergyInterval[k]) ;
981      k-- ;
982    }
983  }
984}   // end of IntegralPAIxSection
985
986////////////////////////////////////////////////////////////////////////
987//
988// Calculation of the PAI Cerenkov integral cross-section
989// fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
990// and fIntegralCerenkov[0] = mean Cerenkov loss per cm  in keV/cm
991
992void G4PAIxSection::IntegralCerenkov()
993{
994  G4int i, k ;
995   fIntegralCerenkov[fSplineNumber] = 0 ;
996   fIntegralCerenkov[0] = 0 ;
997   k = fIntervalNumber -1 ;
998
999   for( i = fSplineNumber-1 ; i >= 1 ; i-- )
1000   {
1001      if(fSplineEnergy[i] >= fEnergyInterval[k])
1002      {
1003        fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i) ;
1004        // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1005      }
1006      else
1007      {
1008        fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + 
1009                                   SumOverBordCerenkov(i+1,fEnergyInterval[k]) ;
1010        k-- ;
1011        // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1012      }
1013   }
1014
1015}   // end of IntegralCerenkov
1016
1017////////////////////////////////////////////////////////////////////////
1018//
1019// Calculation of the PAI Plasmon integral cross-section
1020// fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
1021// and fIntegralPlasmon[0] = mean plasmon loss per cm  in keV/cm
1022
1023void G4PAIxSection::IntegralPlasmon()
1024{
1025   fIntegralPlasmon[fSplineNumber] = 0 ;
1026   fIntegralPlasmon[0] = 0 ;
1027   G4int k = fIntervalNumber -1 ;
1028   for(G4int i=fSplineNumber-1;i>=1;i--)
1029   {
1030      if(fSplineEnergy[i] >= fEnergyInterval[k])
1031      {
1032        fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i) ;
1033      }
1034      else
1035      {
1036        fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + 
1037                                   SumOverBordPlasmon(i+1,fEnergyInterval[k]) ;
1038        k-- ;
1039      }
1040   }
1041
1042}   // end of IntegralPlasmon
1043
1044//////////////////////////////////////////////////////////////////////
1045//
1046// Calculation the PAI integral cross-section inside
1047// of interval of continuous values of photo-ionisation
1048// cross-section. Parameter  'i' is the number of interval.
1049
1050G4double G4PAIxSection::SumOverInterval( G4int i )
1051{         
1052   G4double x0,x1,y0,yy1,a,b,c,result ;
1053
1054   x0 = fSplineEnergy[i] ;
1055   x1 = fSplineEnergy[i+1] ;
1056   y0 = fDifPAIxSection[i] ;
1057   yy1 = fDifPAIxSection[i+1];
1058   c = x1/x0;
1059   a = log10(yy1/y0)/log10(c) ;
1060   // b = log10(y0) - a*log10(x0) ;
1061   b = y0/pow(x0,a) ;
1062   a += 1 ;
1063   if(a == 0) 
1064   {
1065      result = b*log(x1/x0) ;
1066   }
1067   else
1068   {
1069      result = y0*(x1*pow(c,a-1) - x0)/a ;
1070   }
1071   a++;
1072   if(a == 0) 
1073   {
1074      fIntegralPAIxSection[0] += b*log(x1/x0) ;
1075   }
1076   else
1077   {
1078      fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
1079   }
1080   return result ;
1081
1082} //  end of SumOverInterval
1083
1084/////////////////////////////////
1085
1086G4double G4PAIxSection::SumOverIntervaldEdx( G4int i )
1087{         
1088   G4double x0,x1,y0,yy1,a,b,c,result ;
1089
1090   x0 = fSplineEnergy[i] ;
1091   x1 = fSplineEnergy[i+1] ;
1092   y0 = fDifPAIxSection[i] ;
1093   yy1 = fDifPAIxSection[i+1];
1094   c = x1/x0;
1095   a = log10(yy1/y0)/log10(c) ;
1096   // b = log10(y0) - a*log10(x0) ;
1097   b = y0/pow(x0,a) ;
1098   a += 2 ;
1099   if(a == 0) 
1100   {
1101     result = b*log(x1/x0) ;
1102   }
1103   else
1104   {
1105     result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
1106   }
1107   return result ;
1108
1109} //  end of SumOverInterval
1110
1111//////////////////////////////////////////////////////////////////////
1112//
1113// Calculation the PAI Cerenkov integral cross-section inside
1114// of interval of continuous values of photo-ionisation Cerenkov
1115// cross-section. Parameter  'i' is the number of interval.
1116
1117G4double G4PAIxSection::SumOverInterCerenkov( G4int i )
1118{         
1119   G4double x0,x1,y0,yy1,a,b,c,result ;
1120
1121   x0  = fSplineEnergy[i] ;
1122   x1  = fSplineEnergy[i+1] ;
1123   y0  = fdNdxCerenkov[i] ;
1124   yy1 = fdNdxCerenkov[i+1];
1125   // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1126   //   <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1127
1128   c = x1/x0;
1129   a = log10(yy1/y0)/log10(c) ;
1130   b = y0/pow(x0,a) ;
1131
1132   a += 1.0 ;
1133   if(a == 0) result = b*log(c) ;
1134   else       result = y0*(x1*pow(c,a-1) - x0)/a ;   
1135   a += 1.0 ;
1136
1137   if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0) ;
1138   else         fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
1139   //  G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;   
1140   return result ;
1141
1142} //  end of SumOverInterCerenkov
1143
1144//////////////////////////////////////////////////////////////////////
1145//
1146// Calculation the PAI Plasmon integral cross-section inside
1147// of interval of continuous values of photo-ionisation Plasmon
1148// cross-section. Parameter  'i' is the number of interval.
1149
1150G4double G4PAIxSection::SumOverInterPlasmon( G4int i )
1151{         
1152   G4double x0,x1,y0,yy1,a,b,c,result ;
1153
1154   x0  = fSplineEnergy[i] ;
1155   x1  = fSplineEnergy[i+1] ;
1156   y0  = fdNdxPlasmon[i] ;
1157   yy1 = fdNdxPlasmon[i+1];
1158   c =x1/x0;
1159   a = log10(yy1/y0)/log10(c) ;
1160   // b = log10(y0) - a*log10(x0) ;
1161   b = y0/pow(x0,a) ;
1162
1163   a += 1.0 ;
1164   if(a == 0) result = b*log(x1/x0) ;
1165   else       result = y0*(x1*pow(c,a-1) - x0)/a ;   
1166   a += 1.0 ;
1167
1168   if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0) ;
1169   else         fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
1170   
1171   return result ;
1172
1173} //  end of SumOverInterPlasmon
1174
1175///////////////////////////////////////////////////////////////////////////////
1176//
1177// Integration of PAI cross-section for the case of
1178// passing across border between intervals
1179
1180G4double G4PAIxSection::SumOverBorder( G4int      i , 
1181                                       G4double en0    )
1182{               
1183   G4double x0,x1,y0,yy1,a,b,c,d,e0,result ;
1184
1185   e0 = en0 ;
1186   x0 = fSplineEnergy[i] ;
1187   x1 = fSplineEnergy[i+1] ;
1188   y0 = fDifPAIxSection[i] ;
1189   yy1 = fDifPAIxSection[i+1] ;
1190
1191   c = x1/x0;
1192   d = e0/x0;   
1193   a = log10(yy1/y0)/log10(x1/x0) ;
1194   // b0 = log10(y0) - a*log10(x0) ;
1195   b = y0/pow(x0,a);  // pow(10.,b) ;
1196   
1197   a += 1 ;
1198   if(a == 0)
1199   {
1200      result = b*log(x0/e0) ;
1201   }
1202   else
1203   {
1204      result = y0*(x0 - e0*pow(d,a-1))/a ;
1205   }
1206   a++ ;
1207   if(a == 0)
1208   {
1209      fIntegralPAIxSection[0] += b*log(x0/e0) ;
1210   }
1211   else 
1212   {
1213      fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1214   }
1215   x0 = fSplineEnergy[i - 1] ;
1216   x1 = fSplineEnergy[i - 2] ;
1217   y0 = fDifPAIxSection[i - 1] ;
1218   yy1 = fDifPAIxSection[i - 2] ;
1219
1220   c = x1/x0;
1221   d = e0/x0;   
1222   a = log10(yy1/y0)/log10(x1/x0) ;
1223   //  b0 = log10(y0) - a*log10(x0) ;
1224   b = y0/pow(x0,a) ;
1225   a += 1 ;
1226   if(a == 0)
1227   {
1228      result += b*log(e0/x0) ;
1229   }
1230   else
1231   {
1232      result += y0*(e0*pow(d,a-1) - x0)/a ;
1233   }
1234   a++ ;
1235   if(a == 0) 
1236   {
1237      fIntegralPAIxSection[0] += b*log(e0/x0) ;
1238   }
1239   else
1240   {
1241      fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1242   }
1243   return result ;
1244
1245} 
1246
1247///////////////////////////////////////////////////////////////////////
1248
1249G4double G4PAIxSection::SumOverBorderdEdx( G4int      i , 
1250                                       G4double en0    )
1251{               
1252   G4double x0,x1,y0,yy1,a,b,c,d,e0,result ;
1253
1254   e0 = en0 ;
1255   x0 = fSplineEnergy[i] ;
1256   x1 = fSplineEnergy[i+1] ;
1257   y0 = fDifPAIxSection[i] ;
1258   yy1 = fDifPAIxSection[i+1] ;
1259
1260   c = x1/x0;
1261   d = e0/x0;   
1262   a = log10(yy1/y0)/log10(x1/x0) ;
1263   // b0 = log10(y0) - a*log10(x0) ;
1264   b = y0/pow(x0,a);  // pow(10.,b) ;
1265   
1266   a += 2 ;
1267   if(a == 0)
1268   {
1269      result = b*log(x0/e0) ;
1270   }
1271   else 
1272   {
1273      result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1274   }
1275   x0 = fSplineEnergy[i - 1] ;
1276   x1 = fSplineEnergy[i - 2] ;
1277   y0 = fDifPAIxSection[i - 1] ;
1278   yy1 = fDifPAIxSection[i - 2] ;
1279
1280   c = x1/x0;
1281   d = e0/x0;   
1282   a = log10(yy1/y0)/log10(x1/x0) ;
1283   //  b0 = log10(y0) - a*log10(x0) ;
1284   b = y0/pow(x0,a) ;
1285   a += 2 ;
1286   if(a == 0) 
1287   {
1288      result += b*log(e0/x0) ;
1289   }
1290   else
1291   {
1292      result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1293   }
1294   return result ;
1295
1296} 
1297
1298///////////////////////////////////////////////////////////////////////////////
1299//
1300// Integration of Cerenkov cross-section for the case of
1301// passing across border between intervals
1302
1303G4double G4PAIxSection::SumOverBordCerenkov( G4int      i , 
1304                                             G4double en0    )
1305{               
1306   G4double x0,x1,y0,yy1,a,b,e0,c,d,result ;
1307
1308   e0 = en0 ;
1309   x0 = fSplineEnergy[i] ;
1310   x1 = fSplineEnergy[i+1] ;
1311   y0 = fdNdxCerenkov[i] ;
1312   yy1 = fdNdxCerenkov[i+1] ;
1313
1314   //  G4cout<<G4endl;
1315   //  G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1316   //     <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1317   c = x1/x0 ;
1318   d = e0/x0 ;
1319   a = log10(yy1/y0)/log10(c) ;
1320   // b0 = log10(y0) - a*log10(x0) ;
1321   b = y0/pow(x0,a); // pow(10.,b0) ;   
1322   
1323   a += 1.0 ;
1324   if( a == 0 ) result = b*log(x0/e0) ;
1325   else         result = y0*(x0 - e0*pow(d,a-1))/a ;   
1326   a += 1.0 ;
1327
1328   if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0) ;
1329   else         fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1330
1331// G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
1332   
1333   x0  = fSplineEnergy[i - 1] ;
1334   x1  = fSplineEnergy[i - 2] ;
1335   y0  = fdNdxCerenkov[i - 1] ;
1336   yy1 = fdNdxCerenkov[i - 2] ;
1337
1338   // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1339   //    <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1340
1341   c = x1/x0 ;
1342   d = e0/x0 ;
1343   a  = log10(yy1/y0)/log10(x1/x0) ;
1344   // b0 = log10(y0) - a*log10(x0) ;
1345   b  =  y0/pow(x0,a);  // pow(10.,b0) ;
1346
1347   a += 1.0 ;
1348   if( a == 0 ) result += b*log(e0/x0) ;
1349   else         result += y0*(e0*pow(d,a-1) - x0 )/a ;
1350   a += 1.0 ;
1351
1352   if( a == 0 )   fIntegralCerenkov[0] += b*log(e0/x0) ;
1353   else           fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1354
1355   // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
1356   // <<b<<"; result = "<<result<<G4endl;   
1357
1358   return result ;
1359
1360} 
1361
1362///////////////////////////////////////////////////////////////////////////////
1363//
1364// Integration of Plasmon cross-section for the case of
1365// passing across border between intervals
1366
1367G4double G4PAIxSection::SumOverBordPlasmon( G4int      i , 
1368                                             G4double en0    )
1369{               
1370   G4double x0,x1,y0,yy1,a,b,c,d,e0,result ;
1371
1372   e0 = en0 ;
1373   x0 = fSplineEnergy[i] ;
1374   x1 = fSplineEnergy[i+1] ;
1375   y0 = fdNdxPlasmon[i] ;
1376   yy1 = fdNdxPlasmon[i+1] ;
1377
1378   c = x1/x0 ;
1379   d = e0/x0 ;   
1380   a = log10(yy1/y0)/log10(c) ;
1381   //  b0 = log10(y0) - a*log10(x0) ;
1382   b = y0/pow(x0,a); //pow(10.,b) ;
1383   
1384   a += 1.0 ;
1385   if( a == 0 ) result = b*log(x0/e0) ;
1386   else         result = y0*(x0 - e0*pow(d,a-1))/a ;   
1387   a += 1.0 ;
1388
1389   if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0) ;
1390   else         fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1391   
1392   x0 = fSplineEnergy[i - 1] ;
1393   x1 = fSplineEnergy[i - 2] ;
1394   y0 = fdNdxPlasmon[i - 1] ;
1395   yy1 = fdNdxPlasmon[i - 2] ;
1396
1397   c = x1/x0 ;
1398   d = e0/x0 ;
1399   a = log10(yy1/y0)/log10(c) ;
1400   // b0 = log10(y0) - a*log10(x0) ;
1401   b = y0/pow(x0,a);// pow(10.,b0) ;
1402
1403   a += 1.0 ;
1404   if( a == 0 ) result += b*log(e0/x0) ;
1405   else         result += y0*(e0*pow(d,a-1) - x0)/a ;
1406   a += 1.0 ;
1407
1408   if( a == 0 )   fIntegralPlasmon[0] += b*log(e0/x0) ;
1409   else           fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1410   
1411   return result ;
1412
1413} 
1414
1415/////////////////////////////////////////////////////////////////////////
1416//
1417//
1418
1419G4double G4PAIxSection::GetStepEnergyLoss( G4double step )
1420{ 
1421  G4int iTransfer  ;
1422  G4long numOfCollisions ;
1423  G4double loss = 0.0 ;
1424  G4double meanNumber, position ;
1425
1426  // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl ;
1427
1428
1429
1430  meanNumber = fIntegralPAIxSection[1]*step ;
1431  numOfCollisions = G4Poisson(meanNumber) ;
1432
1433  //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1434
1435  while(numOfCollisions)
1436  {
1437    position = fIntegralPAIxSection[1]*G4UniformRand() ;
1438
1439    for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1440    {
1441        if( position >= fIntegralPAIxSection[iTransfer] ) break ;
1442    }
1443    loss += fSplineEnergy[iTransfer]  ;
1444    numOfCollisions-- ;
1445  }
1446  // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl ;
1447
1448  return loss ;
1449}
1450
1451/////////////////////////////////////////////////////////////////////////
1452//
1453//
1454
1455G4double G4PAIxSection::GetStepCerenkovLoss( G4double step )
1456{ 
1457  G4int iTransfer  ;
1458  G4long numOfCollisions ;
1459  G4double loss = 0.0 ;
1460  G4double meanNumber, position ;
1461
1462  // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl ;
1463
1464
1465
1466  meanNumber = fIntegralCerenkov[1]*step ;
1467  numOfCollisions = G4Poisson(meanNumber) ;
1468
1469  //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1470
1471  while(numOfCollisions)
1472  {
1473    position = fIntegralCerenkov[1]*G4UniformRand() ;
1474
1475    for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1476    {
1477        if( position >= fIntegralCerenkov[iTransfer] ) break ;
1478    }
1479    loss += fSplineEnergy[iTransfer]  ;
1480    numOfCollisions-- ;
1481  }
1482  // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl ;
1483
1484  return loss ;
1485}
1486
1487/////////////////////////////////////////////////////////////////////////
1488//
1489//
1490
1491G4double G4PAIxSection::GetStepPlasmonLoss( G4double step )
1492{ 
1493  G4int iTransfer  ;
1494  G4long numOfCollisions ;
1495  G4double loss = 0.0 ;
1496  G4double meanNumber, position ;
1497
1498  // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl ;
1499
1500
1501
1502  meanNumber = fIntegralPlasmon[1]*step ;
1503  numOfCollisions = G4Poisson(meanNumber) ;
1504
1505  //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1506
1507  while(numOfCollisions)
1508  {
1509    position = fIntegralPlasmon[1]*G4UniformRand() ;
1510
1511    for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1512    {
1513        if( position >= fIntegralPlasmon[iTransfer] ) break ;
1514    }
1515    loss += fSplineEnergy[iTransfer]  ;
1516    numOfCollisions-- ;
1517  }
1518  // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl ;
1519
1520  return loss ;
1521}
1522
1523
1524
1525/////////////////////////////////////////////////////////////////////////////
1526//
1527// Init  array of Lorentz factors
1528//
1529
1530G4int G4PAIxSection::fNumberOfGammas = 111 ;
1531
1532const G4double G4PAIxSection::fLorentzFactor[112] =     // fNumberOfGammas+1
1533{
15340.0,
15351.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
15361.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
15371.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
15381.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
15392.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
15403.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
15415.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
15428.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
15431.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
15442.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
15455.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
15461.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
15471.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
15483.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
15496.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
15501.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
15512.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
15524.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
15538.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
15541.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
15553.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
15565.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
15571.065799e+05
1558} ;
1559
1560///////////////////////////////////////////////////////////////////////
1561//
1562// The number of gamma for creation of  spline (near ion-min , G ~ 4 )
1563//
1564
1565const
1566G4int G4PAIxSection::fRefGammaNumber = 29 ; 
1567
1568   
1569//   
1570// end of G4PAIxSection implementation file
1571//
1572////////////////////////////////////////////////////////////////////////////
1573
Note: See TracBrowser for help on using the repository browser.