Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4PAIySection.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/standard/src/G4PAIySection.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4PAIySection.cc (Version 11.0)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 //                                                 27 // 
 28 // G4PAIySection.cc -- class implementation fi     28 // G4PAIySection.cc -- class implementation file
 29 //                                                 29 //
 30 // GEANT 4 class implementation file               30 // GEANT 4 class implementation file
 31 //                                                 31 //
 32 // For information related to this code, pleas     32 // For information related to this code, please, contact
 33 // the Geant4 Collaboration.                       33 // the Geant4 Collaboration.
 34 //                                                 34 //
 35 // R&D: Vladimir.Grichine@cern.ch                  35 // R&D: Vladimir.Grichine@cern.ch
 36 //                                                 36 //
 37 // History:                                        37 // History:
 38 //                                                 38 //
 39 // 01.10.07, V.Ivanchenko create using V.Grich     39 // 01.10.07, V.Ivanchenko create using V.Grichine G4PAIxSection class
 40 // 26.07.09, V.Ivanchenko added protection for     40 // 26.07.09, V.Ivanchenko added protection for mumerical exceptions for 
 41 //              low-density materials              41 //              low-density materials
 42 // 21.11.10 V. Grichine bug fixed in Initialis     42 // 21.11.10 V. Grichine bug fixed in Initialise for reading sandia table from
 43 //            material. Warning: the table is      43 //            material. Warning: the table is tuned for photo-effect not PAI model.
 44 // 23.06.13 V.Grichine arrays->G4DataVectors       44 // 23.06.13 V.Grichine arrays->G4DataVectors
 45 //                                                 45 //
 46                                                    46 
 47 #include "G4PAIySection.hh"                        47 #include "G4PAIySection.hh"
 48                                                    48 
 49 #include "globals.hh"                              49 #include "globals.hh"
 50 #include "G4PhysicalConstants.hh"                  50 #include "G4PhysicalConstants.hh"
 51 #include "G4SystemOfUnits.hh"                      51 #include "G4SystemOfUnits.hh"
 52 #include "G4ios.hh"                                52 #include "G4ios.hh"
 53 #include "G4Poisson.hh"                            53 #include "G4Poisson.hh"
 54 #include "G4Material.hh"                           54 #include "G4Material.hh"
 55 #include "G4MaterialCutsCouple.hh"                 55 #include "G4MaterialCutsCouple.hh"
 56 #include "G4SandiaTable.hh"                        56 #include "G4SandiaTable.hh"
 57 #include "G4Exp.hh"                                57 #include "G4Exp.hh"
 58 #include "G4Log.hh"                                58 #include "G4Log.hh"
 59                                                    59 
 60 using namespace std;                               60 using namespace std;
 61                                                    61 
 62 // Local class constants                           62 // Local class constants
 63                                                    63 
 64 const G4double G4PAIySection::fDelta = 0.005;      64 const G4double G4PAIySection::fDelta = 0.005; // energy shift from interval border
 65 const G4double G4PAIySection::fError = 0.005;      65 const G4double G4PAIySection::fError = 0.005; // error in lin-log approximation
 66                                                    66 
 67 const G4int G4PAIySection::fMaxSplineSize = 50     67 const G4int G4PAIySection::fMaxSplineSize = 500;  // Max size of output spline
 68                                                    68                                                   // arrays
 69                                                    69 
 70 //////////////////////////////////////////////     70 //////////////////////////////////////////////////////////////////
 71 //                                                 71 //
 72 // Constructor                                     72 // Constructor
 73 //                                                 73 //
 74                                                    74 
 75 G4PAIySection::G4PAIySection()                     75 G4PAIySection::G4PAIySection()
 76 {                                                  76 {
 77   fSandia = nullptr;                           <<  77   fSandia = 0;
 78   fDensity = fElectronDensity = fNormalization     78   fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
 79   fIntervalNumber = fSplineNumber = 0;             79   fIntervalNumber = fSplineNumber = 0;
 80   fVerbose = 0;                                    80   fVerbose = 0;
 81                                                    81 
 82   betaBohr = fine_structure_const;                 82   betaBohr = fine_structure_const;
 83   G4double cofBetaBohr = 4.0;                      83   G4double cofBetaBohr = 4.0;
 84   G4double betaBohr2 = fine_structure_const*fi <<  84   G4double betaBohr2   = fine_structure_const*fine_structure_const;
 85   betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr; <<  85   betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr;
 86                                                    86   
 87   fSplineEnergy          = G4DataVector(fMaxSp     87   fSplineEnergy          = G4DataVector(fMaxSplineSize,0.0);
 88   fRePartDielectricConst = G4DataVector(fMaxSp     88   fRePartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
 89   fImPartDielectricConst = G4DataVector(fMaxSp     89   fImPartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
 90   fIntegralTerm          = G4DataVector(fMaxSp     90   fIntegralTerm          = G4DataVector(fMaxSplineSize,0.0);
 91   fDifPAIySection        = G4DataVector(fMaxSp     91   fDifPAIySection        = G4DataVector(fMaxSplineSize,0.0);
 92   fdNdxCerenkov          = G4DataVector(fMaxSp     92   fdNdxCerenkov          = G4DataVector(fMaxSplineSize,0.0);
 93   fdNdxPlasmon           = G4DataVector(fMaxSp     93   fdNdxPlasmon           = G4DataVector(fMaxSplineSize,0.0);
 94   fIntegralPAIySection   = G4DataVector(fMaxSp     94   fIntegralPAIySection   = G4DataVector(fMaxSplineSize,0.0);
 95   fIntegralPAIdEdx       = G4DataVector(fMaxSp     95   fIntegralPAIdEdx       = G4DataVector(fMaxSplineSize,0.0);
 96   fIntegralCerenkov      = G4DataVector(fMaxSp     96   fIntegralCerenkov      = G4DataVector(fMaxSplineSize,0.0);
 97   fIntegralPlasmon       = G4DataVector(fMaxSp     97   fIntegralPlasmon       = G4DataVector(fMaxSplineSize,0.0);
 98                                                    98 
 99   for( G4int i = 0; i < 500; ++i )                 99   for( G4int i = 0; i < 500; ++i ) 
100   {                                               100   {
101     for( G4int j = 0; j < 112; ++j ) { fPAItab    101     for( G4int j = 0; j < 112; ++j ) { fPAItable[i][j] = 0.0; }
102   }                                               102   }
103 }                                                 103 }
104                                                   104 
105 //////////////////////////////////////////////    105 ////////////////////////////////////////////////////////////////////////////
106 //                                                106 //
107 //                                                107 // 
108                                                   108 
109 G4double G4PAIySection::GetLorentzFactor(G4int    109 G4double G4PAIySection::GetLorentzFactor(G4int j) const
110 {                                                 110 {
111    return fLorentzFactor[j];                      111    return fLorentzFactor[j];
112 }                                                 112 }
113                                                   113 
114 //////////////////////////////////////////////    114 ////////////////////////////////////////////////////////////////////////
115 //                                                115 //
116 // Constructor with beta*gamma square value ca    116 // Constructor with beta*gamma square value called from G4PAIModel
117                                                   117 
118 void G4PAIySection::Initialize( const G4Materi    118 void G4PAIySection::Initialize( const G4Material* material,
119                                 G4double maxEn    119                                 G4double maxEnergyTransfer,
120                                 G4double betaG    120                                 G4double betaGammaSq, 
121                                 G4SandiaTable*    121                                 G4SandiaTable* sandia)
122 {                                                 122 {
123   if(fVerbose > 0)                                123   if(fVerbose > 0)
124   {                                               124   {
125     G4cout<<G4endl;                               125     G4cout<<G4endl;
126     G4cout<<"G4PAIySection::Initialize(...,G4S    126     G4cout<<"G4PAIySection::Initialize(...,G4SandiaTable* sandia)"<<G4endl;
127     G4cout<<G4endl;                               127     G4cout<<G4endl;
128   }                                               128   }
129   G4int i, j;                                     129   G4int i, j;
130                                                   130 
131   fSandia          = sandia;                      131   fSandia          = sandia;
132   fIntervalNumber  = sandia->GetMaxInterval();    132   fIntervalNumber  = sandia->GetMaxInterval();
133   fDensity         = material->GetDensity();      133   fDensity         = material->GetDensity();
134   fElectronDensity = material->GetElectronDens    134   fElectronDensity = material->GetElectronDensity();
135                                                   135 
136   // fIntervalNumber--;                           136   // fIntervalNumber--;
137                                                   137 
138   if( fVerbose > 0 )                              138   if( fVerbose > 0 )
139   {                                               139   {
140     G4cout<<"fDensity = "<<fDensity<<"\t"<<fEl    140     G4cout<<"fDensity = "<<fDensity<<"\t"<<fElectronDensity<<"\t fIntervalNumber = "
141           <<fIntervalNumber<< " (beta*gamma)^2    141           <<fIntervalNumber<< " (beta*gamma)^2= " << betaGammaSq << G4endl;
142   }                                               142   }  
143   fEnergyInterval = G4DataVector(fIntervalNumb    143   fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
144   fA1             = G4DataVector(fIntervalNumb    144   fA1             = G4DataVector(fIntervalNumber+2,0.0);
145   fA2             = G4DataVector(fIntervalNumb    145   fA2             = G4DataVector(fIntervalNumber+2,0.0);
146   fA3             = G4DataVector(fIntervalNumb    146   fA3             = G4DataVector(fIntervalNumber+2,0.0);
147   fA4             = G4DataVector(fIntervalNumb    147   fA4             = G4DataVector(fIntervalNumber+2,0.0);
148                                                   148 
149   for( i = 1; i <= fIntervalNumber; ++i )         149   for( i = 1; i <= fIntervalNumber; ++i ) 
150   {                                               150   {
151     if ( sandia->GetSandiaMatTablePAI(i-1,0) <    151     if ( sandia->GetSandiaMatTablePAI(i-1,0) < 1.*eV ) 
152     {                                             152     { 
153       fIntervalNumber--;                          153       fIntervalNumber--;
154       continue;                                   154       continue;
155     }                                             155     }
156     if( ( sandia->GetSandiaMatTablePAI(i-1,0)     156     if( ( sandia->GetSandiaMatTablePAI(i-1,0) >= maxEnergyTransfer ) 
157         || i >= fIntervalNumber )                 157         || i >= fIntervalNumber ) 
158     {                                             158     {
159       fEnergyInterval[i] = maxEnergyTransfer;     159       fEnergyInterval[i] = maxEnergyTransfer;
160       fIntervalNumber = i;                        160       fIntervalNumber = i;
161       break;                                      161       break;
162     }                                             162     }
163     fEnergyInterval[i] = sandia->GetSandiaMatT    163     fEnergyInterval[i] = sandia->GetSandiaMatTablePAI(i-1,0);
164     fA1[i]             = sandia->GetSandiaMatT    164     fA1[i]             = sandia->GetSandiaMatTablePAI(i-1,1);
165     fA2[i]             = sandia->GetSandiaMatT    165     fA2[i]             = sandia->GetSandiaMatTablePAI(i-1,2);
166     fA3[i]             = sandia->GetSandiaMatT    166     fA3[i]             = sandia->GetSandiaMatTablePAI(i-1,3);
167     fA4[i]             = sandia->GetSandiaMatT    167     fA4[i]             = sandia->GetSandiaMatTablePAI(i-1,4);
168                                                   168 
169     if( fVerbose > 0 ) {                          169     if( fVerbose > 0 ) {
170       G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<    170       G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
171             <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4en    171             <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
172     }                                             172     }
173   }                                               173   }
174   if( fVerbose > 0 ) {                            174   if( fVerbose > 0 ) { 
175     G4cout<<"last i = "<<i<<"; "<<"fIntervalNu    175     G4cout<<"last i = "<<i<<"; "<<"fIntervalNumber = "
176           <<fIntervalNumber<<G4endl;              176           <<fIntervalNumber<<G4endl;   
177   }                                               177   }
178   if( fEnergyInterval[fIntervalNumber] != maxE    178   if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer )
179   {                                               179   {
180       fIntervalNumber++;                          180       fIntervalNumber++;
181       fEnergyInterval[fIntervalNumber] = maxEn    181       fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
182   }                                               182   }
183   if( fVerbose > 0 )                              183   if( fVerbose > 0 )
184   {                                               184   {  
185     for( i = 1; i <= fIntervalNumber; ++i )       185     for( i = 1; i <= fIntervalNumber; ++i )
186     {                                             186     {
187       G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<    187       G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
188         <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;     188         <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
189     }                                             189     }
190   }                                               190   }  
191   if( fVerbose > 0 ) {                            191   if( fVerbose > 0 ) {
192     G4cout<<"Now checking, if two borders are     192     G4cout<<"Now checking, if two borders are too close together"<<G4endl;
193   }                                               193   }
194   for( i = 1; i < fIntervalNumber; ++i )          194   for( i = 1; i < fIntervalNumber; ++i )
195   {                                               195   {
196     if( fEnergyInterval[i+1]-fEnergyInterval[i    196     if( fEnergyInterval[i+1]-fEnergyInterval[i] >
197          1.5*fDelta*(fEnergyInterval[i+1]+fEne    197          1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]) ) continue;
198     else                                          198     else
199     {                                             199     {
200       for( j = i; j < fIntervalNumber; j++ )      200       for( j = i; j < fIntervalNumber; j++ )
201       {                                           201       {
202               fEnergyInterval[j] = fEnergyInte    202               fEnergyInterval[j] = fEnergyInterval[j+1];
203               fA1[j]             = fA1[j+1];      203               fA1[j]             = fA1[j+1];
204               fA2[j]             = fA2[j+1];      204               fA2[j]             = fA2[j+1];
205               fA3[j]             = fA3[j+1];      205               fA3[j]             = fA3[j+1];
206               fA4[j]             = fA4[j+1];      206               fA4[j]             = fA4[j+1];
207       }                                           207       }
208       fIntervalNumber--;                          208       fIntervalNumber--;
209     }                                             209     }
210   }                                               210   }
211   if( fVerbose > 0 )                              211   if( fVerbose > 0 )
212   {                                               212   {
213     for( i = 1; i <= fIntervalNumber; ++i )       213     for( i = 1; i <= fIntervalNumber; ++i )
214     {                                             214     {
215       G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<    215       G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
216         <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;     216         <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
217     }                                             217     }
218   }                                               218   }
219   // Preparation of fSplineEnergy array corres    219   // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
220                                                   220 
221   ComputeLowEnergyCof(material);                  221   ComputeLowEnergyCof(material);
222                                                   222             
223   G4double   betaGammaSqRef =                     223   G4double   betaGammaSqRef = 
224     fLorentzFactor[fRefGammaNumber]*fLorentzFa    224     fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
225                                                   225 
226   NormShift(betaGammaSqRef);                      226   NormShift(betaGammaSqRef);             
227   SplainPAI(betaGammaSqRef);                      227   SplainPAI(betaGammaSqRef);
228                                                   228       
229   // Preparation of integral PAI cross section    229   // Preparation of integral PAI cross section for input betaGammaSq
230                                                   230    
231   for( i = 1; i <= fSplineNumber; ++i )           231   for( i = 1; i <= fSplineNumber; ++i )
232   {                                               232   {
233      fDifPAIySection[i] = DifPAIySection(i,bet    233      fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
234                                                   234 
235      if( fVerbose > 0 ) G4cout<<i<<"; dNdxPAI     235      if( fVerbose > 0 ) G4cout<<i<<"; dNdxPAI = "<<fDifPAIySection[i]<<G4endl;
236   }                                               236   }
237   IntegralPAIySection();                          237   IntegralPAIySection();   
238 }                                                 238 }
239                                                   239 
240 //////////////////////////////////////////////    240 /////////////////////////////////////////////////////////////////////////
241 //                                                241 //
242 // Compute low energy cof. It reduces PAI xsc     242 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
243 //                                                243 //
244                                                   244 
245 void G4PAIySection::ComputeLowEnergyCof(const     245 void G4PAIySection::ComputeLowEnergyCof(const G4Material* material)
246 {                                                 246 {    
247   G4int i, numberOfElements = (G4int)material- << 247   G4int i, numberOfElements = material->GetNumberOfElements();
248   G4double sumZ = 0., sumCof = 0.;                248   G4double sumZ = 0., sumCof = 0.; 
249                                                   249 
250   static const G4double p0 =  1.20923e+00;        250   static const G4double p0 =  1.20923e+00; 
251   static const G4double p1 =  3.53256e-01;        251   static const G4double p1 =  3.53256e-01; 
252   static const G4double p2 = -1.45052e-03;        252   static const G4double p2 = -1.45052e-03; 
253                                                   253   
254   G4double* thisMaterialZ   = new G4double[num    254   G4double* thisMaterialZ   = new G4double[numberOfElements];
255   G4double* thisMaterialCof = new G4double[num    255   G4double* thisMaterialCof = new G4double[numberOfElements];
256                                                   256    
257   for( i = 0; i < numberOfElements; ++i )         257   for( i = 0; i < numberOfElements; ++i )
258   {                                               258   {
259     thisMaterialZ[i] = material->GetElement(i)    259     thisMaterialZ[i] = material->GetElement(i)->GetZ();
260     sumZ += thisMaterialZ[i];                     260     sumZ += thisMaterialZ[i];
261     thisMaterialCof[i] = p0+p1*thisMaterialZ[i    261     thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];   
262   }                                               262   }
263   for( i = 0; i < numberOfElements; ++i )         263   for( i = 0; i < numberOfElements; ++i )
264   {                                               264   {
265     sumCof += thisMaterialCof[i]*thisMaterialZ    265     sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
266   }                                               266   }
267   fLowEnergyCof = sumCof;                         267   fLowEnergyCof = sumCof;
268   delete [] thisMaterialZ;                        268   delete [] thisMaterialZ;
269   delete [] thisMaterialCof;                      269   delete [] thisMaterialCof;
270   // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof    270   // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
271 }                                                 271 }
272                                                   272 
273 //////////////////////////////////////////////    273 /////////////////////////////////////////////////////////////////////////
274 //                                                274 //
275 // General control function for class G4PAIySe    275 // General control function for class G4PAIySection
276 //                                                276 //
277                                                   277 
278 void G4PAIySection::InitPAI()                     278 void G4PAIySection::InitPAI()
279 {                                                 279 {    
280    G4int i;                                       280    G4int i;
281    G4double betaGammaSq = fLorentzFactor[fRefG    281    G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
282                           fLorentzFactor[fRefG    282                           fLorentzFactor[fRefGammaNumber] - 1;
283                                                   283 
284    // Preparation of integral PAI cross sectio    284    // Preparation of integral PAI cross section for reference gamma
285                                                   285    
286    NormShift(betaGammaSq);                        286    NormShift(betaGammaSq);             
287    SplainPAI(betaGammaSq);                        287    SplainPAI(betaGammaSq);
288                                                   288 
289    IntegralPAIySection();                         289    IntegralPAIySection();
290    IntegralCerenkov();                            290    IntegralCerenkov();
291    IntegralPlasmon();                             291    IntegralPlasmon();
292                                                   292 
293    for( i = 0; i<= fSplineNumber; ++i)            293    for( i = 0; i<= fSplineNumber; ++i)
294    {                                              294    {
295      fPAItable[i][fRefGammaNumber] = fIntegral    295      fPAItable[i][fRefGammaNumber] = fIntegralPAIySection[i];
296                                                   296       
297      if(i != 0)  fPAItable[i][0] = fSplineEner    297      if(i != 0)  fPAItable[i][0] = fSplineEnergy[i];     
298    }                                              298    }
299    fPAItable[0][0] = fSplineNumber;               299    fPAItable[0][0] = fSplineNumber;
300                                                   300    
301    for( G4int j = 1; j < 112; ++j)       // fo    301    for( G4int j = 1; j < 112; ++j)       // for other gammas
302    {                                              302    {
303       if( j == fRefGammaNumber ) continue;        303       if( j == fRefGammaNumber ) continue;
304                                                   304       
305       betaGammaSq = fLorentzFactor[j]*fLorentz    305       betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
306                                                   306       
307       for(i = 1; i <= fSplineNumber; ++i)         307       for(i = 1; i <= fSplineNumber; ++i)
308       {                                           308       {
309          fDifPAIySection[i] = DifPAIySection(i    309          fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
310          fdNdxCerenkov[i]   = PAIdNdxCerenkov(    310          fdNdxCerenkov[i]   = PAIdNdxCerenkov(i,betaGammaSq);
311          fdNdxPlasmon[i]    = PAIdNdxPlasmon(i    311          fdNdxPlasmon[i]    = PAIdNdxPlasmon(i,betaGammaSq);
312       }                                           312       }
313       IntegralPAIySection();                      313       IntegralPAIySection();
314       IntegralCerenkov();                         314       IntegralCerenkov();
315       IntegralPlasmon();                          315       IntegralPlasmon();
316                                                   316       
317       for(i = 0; i <= fSplineNumber; ++i)         317       for(i = 0; i <= fSplineNumber; ++i)
318       {                                           318       {
319         fPAItable[i][j] = fIntegralPAIySection    319         fPAItable[i][j] = fIntegralPAIySection[i];
320       }                                           320       }
321    }                                              321    } 
322 }                                                 322 }  
323                                                   323 
324 //////////////////////////////////////////////    324 ///////////////////////////////////////////////////////////////////////
325 //                                                325 //
326 // Shifting from borders to intervals Creation    326 // Shifting from borders to intervals Creation of first energy points
327 //                                                327 //
328                                                   328 
329 void G4PAIySection::NormShift(G4double betaGam    329 void G4PAIySection::NormShift(G4double betaGammaSq)
330 {                                                 330 {
331   G4int i, j;                                     331   G4int i, j;
332                                                   332 
333   for( i = 1; i <= fIntervalNumber-1; ++i)        333   for( i = 1; i <= fIntervalNumber-1; ++i)
334   {                                               334   {
335     for( j = 1; j <= 2; ++j)                      335     for( j = 1; j <= 2; ++j)
336     {                                             336     {
337       fSplineNumber = (i-1)*2 + j;                337       fSplineNumber = (i-1)*2 + j;
338                                                   338 
339       if( j == 1 ) fSplineEnergy[fSplineNumber    339       if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i  ]*(1+fDelta);
340       else         fSplineEnergy[fSplineNumber    340       else         fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta); 
341       //    G4cout<<"cn = "<<fSplineNumber<<";    341       //    G4cout<<"cn = "<<fSplineNumber<<"; "<<"energy = "
342       //  <<fSplineEnergy[fSplineNumber]<<G4en    342       //  <<fSplineEnergy[fSplineNumber]<<G4endl;
343     }                                             343     }
344   }                                               344   }
345   fIntegralTerm[1]=RutherfordIntegral(1,fEnerg    345   fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
346                                                   346 
347   j = 1;                                          347   j = 1;
348                                                   348 
349   for(i=2;i<=fSplineNumber;++i)                   349   for(i=2;i<=fSplineNumber;++i)
350   {                                               350   {
351     if(fSplineEnergy[i]<fEnergyInterval[j+1])     351     if(fSplineEnergy[i]<fEnergyInterval[j+1])
352     {                                             352     {
353          fIntegralTerm[i] = fIntegralTerm[i-1]    353          fIntegralTerm[i] = fIntegralTerm[i-1] + 
354                             RutherfordIntegral    354                             RutherfordIntegral(j,fSplineEnergy[i-1],
355                                                   355                                                  fSplineEnergy[i]   );
356     }                                             356     }
357     else                                          357     else
358     {                                             358     {
359        G4double x = RutherfordIntegral(j,fSpli    359        G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
360                                            fEn    360                                            fEnergyInterval[j+1]   );
361          j++;                                     361          j++;
362          fIntegralTerm[i] = fIntegralTerm[i-1]    362          fIntegralTerm[i] = fIntegralTerm[i-1] + x + 
363                             RutherfordIntegral    363                             RutherfordIntegral(j,fEnergyInterval[j],
364                                                   364                                                  fSplineEnergy[i]    );
365     }                                             365     }
366     // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t"    366     // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t"<<fIntegralTerm[i]<<"\n"<<G4endl;
367   }                                               367   } 
368   static const G4double nfactor =                 368   static const G4double nfactor =
369     2*pi*pi*hbarc*hbarc*fine_structure_const/e    369     2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
370   fNormalizationCof = nfactor*fElectronDensity    370   fNormalizationCof = nfactor*fElectronDensity/fIntegralTerm[fSplineNumber];
371                                                   371 
372   // G4cout<<"fNormalizationCof = "<<fNormaliz    372   // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl;
373                                                   373 
374   // Calculation of PAI differrential cross-se << 374           // Calculation of PAI differrential cross-section (1/(keV*cm))
375   // in the energy points near borders of ener << 375           // in the energy points near borders of energy intervals
376                                                   376 
377   for(G4int k=1; k<=fIntervalNumber-1; ++k)    << 377    for(G4int k=1; k<=fIntervalNumber-1; ++k)
378    {                                              378    {
379      for(j=1; j<=2; ++j)                       << 379       for(j=1; j<=2; ++j)
380       {                                           380       {
381          i = (k-1)*2 + j;                         381          i = (k-1)*2 + j;
382          fImPartDielectricConst[i] = fNormaliz    382          fImPartDielectricConst[i] = fNormalizationCof*
383                                      ImPartDie    383                                      ImPartDielectricConst(k,fSplineEnergy[i]);
384          fRePartDielectricConst[i] = fNormaliz    384          fRePartDielectricConst[i] = fNormalizationCof*
385                                      RePartDie    385                                      RePartDielectricConst(fSplineEnergy[i]);
386          fIntegralTerm[i] *= fNormalizationCof    386          fIntegralTerm[i] *= fNormalizationCof;
387                                                   387 
388          fDifPAIySection[i] = DifPAIySection(i    388          fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
389          fdNdxCerenkov[i]   = PAIdNdxCerenkov(    389          fdNdxCerenkov[i]   = PAIdNdxCerenkov(i,betaGammaSq);
390          fdNdxPlasmon[i]    = PAIdNdxPlasmon(i    390          fdNdxPlasmon[i]    = PAIdNdxPlasmon(i,betaGammaSq);
391       }                                           391       }
392    }                                              392    }
393                                                   393 
394 }  // end of NormShift                            394 }  // end of NormShift 
395                                                   395 
396 //////////////////////////////////////////////    396 /////////////////////////////////////////////////////////////////////////
397 //                                                397 //
398 // Creation of new energy points as geometrica    398 // Creation of new energy points as geometrical mean of existing
399 // one, calculation PAI_cs for them, while the    399 // one, calculation PAI_cs for them, while the error of logarithmic
400 // linear approximation would be smaller than     400 // linear approximation would be smaller than 'fError'
401                                                   401 
402 void G4PAIySection::SplainPAI(G4double betaGam    402 void G4PAIySection::SplainPAI(G4double betaGammaSq)
403 {                                                 403 {
404    G4int k = 1;                                   404    G4int k = 1;
405    G4int i = 1;                                   405    G4int i = 1;
406                                                   406 
407    while ( (i < fSplineNumber) && (fSplineNumb    407    while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
408    {                                              408    {
409       if(fSplineEnergy[i+1] > fEnergyInterval[    409       if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
410       {                                           410       {
411           k++;   // Here next energy point is     411           k++;   // Here next energy point is in next energy interval
412           ++i;                                    412           ++i;
413           continue;                               413           continue;
414       }                                           414       }
415       // Shifting of arrayes for inserting the << 415                         // Shifting of arrayes for inserting the geometrical 
416       // average of 'i' and 'i+1' energy point << 416                        // average of 'i' and 'i+1' energy points to 'i+1' place
417       fSplineNumber++;                            417       fSplineNumber++;
418                                                   418 
419       for(G4int j = fSplineNumber; j >= i+2; j    419       for(G4int j = fSplineNumber; j >= i+2; j-- )
420       {                                           420       {
421          fSplineEnergy[j]          = fSplineEn    421          fSplineEnergy[j]          = fSplineEnergy[j-1];
422          fImPartDielectricConst[j] = fImPartDi    422          fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
423          fRePartDielectricConst[j] = fRePartDi    423          fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
424          fIntegralTerm[j]          = fIntegral    424          fIntegralTerm[j]          = fIntegralTerm[j-1];
425                                                   425 
426          fDifPAIySection[j] = fDifPAIySection[    426          fDifPAIySection[j] = fDifPAIySection[j-1];
427          fdNdxCerenkov[j]   = fdNdxCerenkov[j-    427          fdNdxCerenkov[j]   = fdNdxCerenkov[j-1];
428          fdNdxPlasmon[j]    = fdNdxPlasmon[j-1    428          fdNdxPlasmon[j]    = fdNdxPlasmon[j-1];
429       }                                           429       }
430       G4double x1  = fSplineEnergy[i];            430       G4double x1  = fSplineEnergy[i];
431       G4double x2  = fSplineEnergy[i+1];          431       G4double x2  = fSplineEnergy[i+1];
432       G4double yy1 = fDifPAIySection[i];          432       G4double yy1 = fDifPAIySection[i];
433       G4double y2  = fDifPAIySection[i+1];        433       G4double y2  = fDifPAIySection[i+1];
434                                                   434 
435       G4double en1 = sqrt(x1*x2);                 435       G4double en1 = sqrt(x1*x2);
436       fSplineEnergy[i+1] = en1;                   436       fSplineEnergy[i+1] = en1;
437                                                   437 
438       // Calculation of logarithmic linear app << 438                  // Calculation of logarithmic linear approximation
439       // in this (enr) energy point, which num << 439                  // in this (enr) energy point, which number is 'i+1' now
440                                                   440 
441       G4double a = log10(y2/yy1)/log10(x2/x1);    441       G4double a = log10(y2/yy1)/log10(x2/x1);
442       G4double b = log10(yy1) - a*log10(x1);      442       G4double b = log10(yy1) - a*log10(x1);
443       G4double y = a*log10(en1) + b;              443       G4double y = a*log10(en1) + b;
444       y = pow(10.,y);                             444       y = pow(10.,y);
445                                                   445 
446       // Calculation of the PAI dif. cross-sec << 446                  // Calculation of the PAI dif. cross-section at this point
447                                                   447 
448       fImPartDielectricConst[i+1] = fNormaliza    448       fImPartDielectricConst[i+1] = fNormalizationCof*
449                                     ImPartDiel    449                                     ImPartDielectricConst(k,fSplineEnergy[i+1]);
450       fRePartDielectricConst[i+1] = fNormaliza    450       fRePartDielectricConst[i+1] = fNormalizationCof*
451                                     RePartDiel    451                                     RePartDielectricConst(fSplineEnergy[i+1]);
452       fIntegralTerm[i+1] = fIntegralTerm[i] +     452       fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
453                            RutherfordIntegral(    453                            RutherfordIntegral(k,fSplineEnergy[i],
454                                                   454                                                 fSplineEnergy[i+1]);
455                                                   455 
456       fDifPAIySection[i+1] = DifPAIySection(i+    456       fDifPAIySection[i+1] = DifPAIySection(i+1,betaGammaSq);
457       fdNdxCerenkov[i+1]   = PAIdNdxCerenkov(i    457       fdNdxCerenkov[i+1]   = PAIdNdxCerenkov(i+1,betaGammaSq);
458       fdNdxPlasmon[i+1]    = PAIdNdxPlasmon(i+    458       fdNdxPlasmon[i+1]    = PAIdNdxPlasmon(i+1,betaGammaSq);
459                                                   459 
460                   // Condition for next divisi    460                   // Condition for next division of this segment or to pass
461                   // to higher energies           461                   // to higher energies
462                                                   462 
463       G4double x = 2*(fDifPAIySection[i+1] - y    463       G4double x = 2*(fDifPAIySection[i+1] - y)/(fDifPAIySection[i+1] + y);
464                                                   464 
465       G4double delta = 2.*(fSplineEnergy[i+1]-    465       G4double delta = 2.*(fSplineEnergy[i+1]-fSplineEnergy[i])
466         /(fSplineEnergy[i+1]+fSplineEnergy[i])    466         /(fSplineEnergy[i+1]+fSplineEnergy[i]);
467                                                   467 
468       if( x < 0 )                                 468       if( x < 0 ) 
469       {                                           469       {
470          x = -x;                                  470          x = -x;
471       }                                           471       }
472       if( x > fError && fSplineNumber < fMaxSp    472       if( x > fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*fDelta  )
473       {                                           473       {
474          continue;  // next division              474          continue;  // next division
475       }                                           475       }
476       i += 2;  // pass to next segment            476       i += 2;  // pass to next segment
477                                                   477 
478       // Loop checking, 03-Aug-2015, Vladimir     478       // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
479    }   // close 'while'                           479    }   // close 'while'
480                                                   480 
481 }  // end of SplainPAI                            481 }  // end of SplainPAI 
482                                                   482 
483                                                   483 
484 //////////////////////////////////////////////    484 ////////////////////////////////////////////////////////////////////
485 //                                                485 //
486 // Integration over electrons that could be co    486 // Integration over electrons that could be considered
487 // quasi-free at energy transfer of interest      487 // quasi-free at energy transfer of interest
488                                                   488 
489 G4double G4PAIySection::RutherfordIntegral( G4    489 G4double G4PAIySection::RutherfordIntegral( G4int k,
490                                             G4    490                                             G4double x1,
491                                                   491                                               G4double x2   )
492 {                                                 492 {
493    G4double  c1, c2, c3;                          493    G4double  c1, c2, c3;
494    // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<    494    // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
495    G4double x12 = x1*x2;                          495    G4double x12 = x1*x2;   
496    c1 = (x2 - x1)/x12;                            496    c1 = (x2 - x1)/x12;
497    c2 = (x2 - x1)*(x2 + x1)/(x12*x12);            497    c2 = (x2 - x1)*(x2 + x1)/(x12*x12);
498    c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/(x12    498    c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/(x12*x12*x12);
499    // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<    499    // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;   
500                                                   500    
501    return  fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3    501    return  fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
502                                                   502 
503 }   // end of RutherfordIntegral                  503 }   // end of RutherfordIntegral 
504                                                   504 
505                                                   505 
506 //////////////////////////////////////////////    506 /////////////////////////////////////////////////////////////////
507 //                                                507 //
508 // Imaginary part of dielectric constant          508 // Imaginary part of dielectric constant
509 // (G4int k - interval number, G4double en1 -     509 // (G4int k - interval number, G4double en1 - energy point)
510                                                   510 
511 G4double G4PAIySection::ImPartDielectricConst( << 511 G4double G4PAIySection::ImPartDielectricConst( G4int    k ,
                                                   >> 512                                                G4double energy1 )
512 {                                                 513 {
513    G4double energy2,energy3,energy4,result;       514    G4double energy2,energy3,energy4,result;
514                                                   515 
515    energy2 = energy1*energy1;                     516    energy2 = energy1*energy1;
516    energy3 = energy2*energy1;                     517    energy3 = energy2*energy1;
517    energy4 = energy3*energy1;                     518    energy4 = energy3*energy1;
518                                                   519    
519    result = fA1[k]/energy1+fA2[k]/energy2+fA3[    520    result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;  
520    result *=hbarc/energy1;                        521    result *=hbarc/energy1;
521                                                   522    
522    return result;                                 523    return result;
523                                                   524 
524 }  // end of ImPartDielectricConst                525 }  // end of ImPartDielectricConst 
525                                                   526 
526                                                   527 
527 //////////////////////////////////////////////    528 //////////////////////////////////////////////////////////////////////////////
528 //                                                529 //
529 // Real part of dielectric constant minus unit    530 // Real part of dielectric constant minus unit: epsilon_1 - 1
530 // (G4double enb - energy point)                  531 // (G4double enb - energy point)
531 //                                                532 //
532                                                   533 
533 G4double G4PAIySection::RePartDielectricConst(    534 G4double G4PAIySection::RePartDielectricConst(G4double enb)
534 {                                                 535 {       
535    G4double x0, x02, x03, x04, x05, x1, x2, xx    536    G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
536             c1, c2, c3, cof1, cof2, xln1, xln2    537             c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
537                                                   538 
538    x0 = enb;                                      539    x0 = enb;
539    result = 0;                                    540    result = 0;
540                                                   541    
541    for(G4int i=1;i<=fIntervalNumber-1;++i)        542    for(G4int i=1;i<=fIntervalNumber-1;++i)
542    {                                              543    {
543       x1 = fEnergyInterval[i];                    544       x1 = fEnergyInterval[i];
544       x2 = fEnergyInterval[i+1];                  545       x2 = fEnergyInterval[i+1];
545       xx1 = x1 - x0;                              546       xx1 = x1 - x0;
546       xx2 = x2 - x0;                              547       xx2 = x2 - x0;
547       xx12 = xx2/xx1;                             548       xx12 = xx2/xx1;
548                                                   549       
549       if(xx12<0.)                              << 550       if(xx12<0)
550       {                                           551       {
551          xx12 = -xx12;                            552          xx12 = -xx12;
552       }                                           553       }
553       xln1 = log(x2/x1);                          554       xln1 = log(x2/x1);
554       xln2 = log(xx12);                           555       xln2 = log(xx12);
555       xln3 = log((x2 + x0)/(x1 + x0));            556       xln3 = log((x2 + x0)/(x1 + x0));
556       x02 = x0*x0;                                557       x02 = x0*x0;
557       x03 = x02*x0;                               558       x03 = x02*x0;
558       x04 = x03*x0;                               559       x04 = x03*x0;
559       x05 = x04*x0;                               560       x05 = x04*x0;
560       G4double x12 = x1*x2;                       561       G4double x12 = x1*x2;
561       c1  = (x2 - x1)/x12;                        562       c1  = (x2 - x1)/x12;
562       c2  = (x2 - x1)*(x2 +x1)/(x12*x12);         563       c2  = (x2 - x1)*(x2 +x1)/(x12*x12);
563       c3  = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/(    564       c3  = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/(x12*x12*x12);
564                                                   565 
565       result -= (fA1[i]/x02 + fA3[i]/x04)*xln1    566       result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
566       result -= (fA2[i]/x02 + fA4[i]/x04)*c1;     567       result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
567       result -= fA3[i]*c2/2/x02;                  568       result -= fA3[i]*c2/2/x02;
568       result -= fA4[i]*c3/3/x02;                  569       result -= fA4[i]*c3/3/x02;
569                                                   570 
570       cof1 = fA1[i]/x02 + fA3[i]/x04;             571       cof1 = fA1[i]/x02 + fA3[i]/x04;
571       cof2 = fA2[i]/x03 + fA4[i]/x05;             572       cof2 = fA2[i]/x03 + fA4[i]/x05;
572                                                   573 
573       result += 0.5*(cof1 +cof2)*xln2;            574       result += 0.5*(cof1 +cof2)*xln2;
574       result += 0.5*(cof1 - cof2)*xln3;           575       result += 0.5*(cof1 - cof2)*xln3;
575    }                                              576    } 
576    result *= 2*hbarc/pi;                          577    result *= 2*hbarc/pi;
577                                                   578    
578    return result;                                 579    return result;
579                                                   580 
580 }   // end of RePartDielectricConst               581 }   // end of RePartDielectricConst 
581                                                   582 
582 //////////////////////////////////////////////    583 //////////////////////////////////////////////////////////////////////
583 //                                                584 //
584 // PAI differential cross-section in terms of     585 // PAI differential cross-section in terms of
585 // simplified Allison's equation                  586 // simplified Allison's equation
586 //                                                587 //
587                                                   588 
588 G4double G4PAIySection::DifPAIySection( G4int     589 G4double G4PAIySection::DifPAIySection( G4int              i ,
589                                         G4doub    590                                         G4double betaGammaSq  )
590 {                                                 591 {        
591    G4double beta, be2,cof,x1,x2,x3,x4,x5,x6,x7 << 592   G4double beta, be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
                                                   >> 593    //G4double beta, be4;
                                                   >> 594    //G4double be4;
                                                   >> 595    // G4double betaBohr2 = fine_structure_const*fine_structure_const;
                                                   >> 596    // G4double betaBohr4 = betaBohr2*betaBohr2*4.0;
592    be2 = betaGammaSq/(1 + betaGammaSq);           597    be2 = betaGammaSq/(1 + betaGammaSq);
593    beta = std::sqrt(be2);                      << 598    //be4 = be2*be2;
                                                   >> 599    beta = sqrt(be2);
594    cof = 1;                                       600    cof = 1;
595    x1 = log(2*electron_mass_c2/fSplineEnergy[i    601    x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
596                                                   602 
597    if( betaGammaSq < 0.01 ) x2 = log(be2);        603    if( betaGammaSq < 0.01 ) x2 = log(be2);
598    else                                           604    else
599    {                                              605    {
600      x2 = -log( (1/betaGammaSq - fRePartDielec    606      x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
601                 (1/betaGammaSq - fRePartDielec    607                 (1/betaGammaSq - fRePartDielectricConst[i]) + 
602                 fImPartDielectricConst[i]*fImP    608                 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
603    }                                              609    }
604    if( fImPartDielectricConst[i] == 0.0 ||beta    610    if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
605    {                                              611    {
606      x6=0;                                        612      x6=0;
607    }                                              613    }
608    else                                           614    else
609    {                                              615    {
610      x3 = -fRePartDielectricConst[i] + 1/betaG    616      x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
611      x5 = -1 - fRePartDielectricConst[i] +        617      x5 = -1 - fRePartDielectricConst[i] +
612           be2*((1 +fRePartDielectricConst[i])*    618           be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
613           fImPartDielectricConst[i]*fImPartDie    619           fImPartDielectricConst[i]*fImPartDielectricConst[i]);
614                                                   620 
615      x7 = std::atan2(fImPartDielectricConst[i] << 621      x7 = atan2(fImPartDielectricConst[i],x3);
616      x6 = x5 * x7;                                622      x6 = x5 * x7;
617    }                                              623    }
                                                   >> 624     // if(fImPartDielectricConst[i] == 0) x6 = 0;
                                                   >> 625    
618    x4 = ((x1 + x2)*fImPartDielectricConst[i] +    626    x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
                                                   >> 627    //   if( x4 < 0.0 ) x4 = 0.0;
619    x8 = (1 + fRePartDielectricConst[i])*(1 + f    628    x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 
620         fImPartDielectricConst[i]*fImPartDiele    629         fImPartDielectricConst[i]*fImPartDielectricConst[i];
621                                                   630 
622    result = (x4 + cof*fIntegralTerm[i]/fSpline    631    result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
623    result = std::max(result, 1.0e-8);             632    result = std::max(result, 1.0e-8);
624    result *= fine_structure_const/(be2*pi);    << 633    result *= fine_structure_const/be2/pi;
625    // low energy correction                       634    // low energy correction
626                                                   635 
627    G4double lowCof = fLowEnergyCof; // 6.0 ; /    636    G4double lowCof = fLowEnergyCof; // 6.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)? 
628                                                   637 
629    result *= (1 - std::exp(-beta/(betaBohr*low << 638    result *= (1 - exp(-beta/betaBohr/lowCof));
                                                   >> 639 
                                                   >> 640    //   result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr));
                                                   >> 641    //  result *= (1-exp(-be2/betaBohr2));
                                                   >> 642    // result *= (1-exp(-be4/betaBohr4));
                                                   >> 643    //   if(fDensity >= 0.1)
630    if(x8 > 0.)                                    644    if(x8 > 0.)
631    {                                              645    { 
632      result /= x8;                             << 646       result /= x8;
633    }                                              647    }
634    return result;                                 648    return result;
635                                                   649 
636 } // end of DifPAIySection                        650 } // end of DifPAIySection 
637                                                   651 
638 //////////////////////////////////////////////    652 //////////////////////////////////////////////////////////////////////////
639 //                                                653 //
640 // Calculation od dN/dx of collisions with cre    654 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
641                                                   655 
642 G4double G4PAIySection::PAIdNdxCerenkov( G4int << 656 G4double G4PAIySection::PAIdNdxCerenkov( G4int    i ,
                                                   >> 657                                          G4double betaGammaSq  )
643 {                                                 658 {        
644    G4double logarithm, x3, x5, argument, modul    659    G4double logarithm, x3, x5, argument, modul2, dNdxC; 
645    G4double be2, be4;                             660    G4double be2, be4;
646                                                   661 
                                                   >> 662    //G4double cof         = 1.0;
                                                   >> 663 
647    be2 = betaGammaSq/(1 + betaGammaSq);           664    be2 = betaGammaSq/(1 + betaGammaSq);
648    be4 = be2*be2;                                 665    be4 = be2*be2;
649                                                   666 
650    if( betaGammaSq < 0.01 ) logarithm = log(1.    667    if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
651    else                                           668    else
652    {                                              669    {
653      logarithm = -std::log( (1/betaGammaSq - f << 670      logarithm  = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
654                         (1/betaGammaSq - fRePa    671                         (1/betaGammaSq - fRePartDielectricConst[i]) + 
655                         fImPartDielectricConst    672                         fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
656      logarithm += std::log(1+1.0/betaGammaSq); << 673      logarithm += log(1+1.0/betaGammaSq);
657    }                                              674    }
658                                                   675 
659    if( fImPartDielectricConst[i] == 0.0 || bet    676    if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
660    {                                              677    {
661      argument = 0.0;                              678      argument = 0.0;
662    }                                              679    }
663    else                                           680    else
664    {                                              681    {
665      x3 = -fRePartDielectricConst[i] + 1.0/bet    682      x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
666      x5 = -1.0 - fRePartDielectricConst[i] +      683      x5 = -1.0 - fRePartDielectricConst[i] +
667           be2*((1.0 +fRePartDielectricConst[i]    684           be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
668           fImPartDielectricConst[i]*fImPartDie    685           fImPartDielectricConst[i]*fImPartDielectricConst[i]);
669      if( x3 == 0.0 ) argument = 0.5*pi;           686      if( x3 == 0.0 ) argument = 0.5*pi;
670      else            argument = std::atan2(fIm << 687      else            argument = atan2(fImPartDielectricConst[i],x3);
671      argument *= x5 ;                             688      argument *= x5 ;
672    }                                              689    }   
673    dNdxC = ( logarithm*fImPartDielectricConst[    690    dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
674                                                   691   
675    if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;             692    if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
676                                                   693 
677    dNdxC *= fine_structure_const/be2/pi;          694    dNdxC *= fine_structure_const/be2/pi;
678                                                   695 
679    dNdxC *= (1 - std::exp(-be4/betaBohr4));    << 696    dNdxC *= (1-exp(-be4/betaBohr4));
680                                                   697 
                                                   >> 698    //   if(fDensity >= 0.1)
                                                   >> 699    // { 
681    modul2 = (1.0 + fRePartDielectricConst[i])*    700    modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) + 
682                     fImPartDielectricConst[i]*    701                     fImPartDielectricConst[i]*fImPartDielectricConst[i];
683    if(modul2 > 0.)                                702    if(modul2 > 0.)
684      {                                            703      {
685        dNdxC /= modul2;                           704        dNdxC /= modul2;
686      }                                            705      }
687    return dNdxC;                                  706    return dNdxC;
688                                                   707 
689 } // end of PAIdNdxCerenkov                       708 } // end of PAIdNdxCerenkov 
690                                                   709 
691 //////////////////////////////////////////////    710 //////////////////////////////////////////////////////////////////////////
692 //                                                711 //
693 // Calculation od dN/dx of collisions with cre    712 // Calculation od dN/dx of collisions with creation of longitudinal EM
694 // excitations (plasmons, delta-electrons)        713 // excitations (plasmons, delta-electrons)
695                                                   714 
696 G4double G4PAIySection::PAIdNdxPlasmon( G4int  << 715 G4double G4PAIySection::PAIdNdxPlasmon( G4int    i ,
                                                   >> 716                                         G4double betaGammaSq  )
697 {                                                 717 {        
698    G4double cof, resonance, modul2, dNdxP;        718    G4double cof, resonance, modul2, dNdxP;
699    G4double be2, be4;                             719    G4double be2, be4;
700                                                   720 
701    cof = 1;                                       721    cof = 1;
702                                                   722 
703    be2 = betaGammaSq/(1 + betaGammaSq);           723    be2 = betaGammaSq/(1 + betaGammaSq);
704    be4 = be2*be2;                                 724    be4 = be2*be2;
705                                                   725  
706    resonance = std::log(2*electron_mass_c2*be2 << 726    resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);  
707    resonance *= fImPartDielectricConst[i]/hbar    727    resonance *= fImPartDielectricConst[i]/hbarc;
708                                                   728 
709    dNdxP = ( resonance + cof*fIntegralTerm[i]/    729    dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
710                                                   730 
711    dNdxP = std::max(dNdxP, 1.0e-8);               731    dNdxP = std::max(dNdxP, 1.0e-8);
712                                                   732 
713    dNdxP *= fine_structure_const/be2/pi;          733    dNdxP *= fine_structure_const/be2/pi;
714    dNdxP *= (1 - std::exp(-be4/betaBohr4));    << 734    dNdxP *= (1-exp(-be4/betaBohr4));
715                                                   735 
                                                   >> 736 //   if( fDensity >= 0.1 )
                                                   >> 737 //   { 
716    modul2 = (1 + fRePartDielectricConst[i])*(1    738    modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 
717      fImPartDielectricConst[i]*fImPartDielectr    739      fImPartDielectricConst[i]*fImPartDielectricConst[i];
718    if(modul2 > 0.)                                740    if(modul2 > 0.)
719      {                                            741      { 
720        dNdxP /= modul2;                           742        dNdxP /= modul2;
721      }                                            743      }
722    return dNdxP;                                  744    return dNdxP;
723                                                   745 
724 } // end of PAIdNdxPlasmon                        746 } // end of PAIdNdxPlasmon 
725                                                   747 
726 //////////////////////////////////////////////    748 ////////////////////////////////////////////////////////////////////////
727 //                                                749 //
728 // Calculation of the PAI integral cross-secti    750 // Calculation of the PAI integral cross-section
729 // fIntegralPAIySection[1] = specific primary     751 // fIntegralPAIySection[1] = specific primary ionisation, 1/cm
730 // and fIntegralPAIySection[0] = mean energy l    752 // and fIntegralPAIySection[0] = mean energy loss per cm  in keV/cm
731                                                   753 
732 void G4PAIySection::IntegralPAIySection()         754 void G4PAIySection::IntegralPAIySection()
733 {                                                 755 {
734   fIntegralPAIySection[fSplineNumber] = 0;        756   fIntegralPAIySection[fSplineNumber] = 0;
735   fIntegralPAIdEdx[fSplineNumber]     = 0;        757   fIntegralPAIdEdx[fSplineNumber]     = 0;
736   fIntegralPAIySection[0]             = 0;        758   fIntegralPAIySection[0]             = 0;
737   G4int k = fIntervalNumber -1;                   759   G4int k = fIntervalNumber -1;
738                                                   760 
739   for(G4int i = fSplineNumber-1; i >= 1; i--)     761   for(G4int i = fSplineNumber-1; i >= 1; i--)
740   {                                               762   {
741     if(fSplineEnergy[i] >= fEnergyInterval[k])    763     if(fSplineEnergy[i] >= fEnergyInterval[k])
742     {                                             764     {
743       fIntegralPAIySection[i] = fIntegralPAIyS    765       fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + SumOverInterval(i);
744       fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i    766       fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
745     }                                             767     }
746     else                                          768     else
747     {                                             769     {
748       fIntegralPAIySection[i] = fIntegralPAIyS    770       fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + 
749                                    SumOverBord    771                                    SumOverBorder(i+1,fEnergyInterval[k]);
750       fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i    772       fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + 
751                                    SumOverBord    773                                    SumOverBorderdEdx(i+1,fEnergyInterval[k]);
752       k--;                                        774       k--;
753     }                                             775     }
754   }                                               776   }
755 }   // end of IntegralPAIySection                 777 }   // end of IntegralPAIySection 
756                                                   778 
757 //////////////////////////////////////////////    779 ////////////////////////////////////////////////////////////////////////
758 //                                                780 //
759 // Calculation of the PAI Cerenkov integral cr    781 // Calculation of the PAI Cerenkov integral cross-section
760 // fIntegralCrenkov[1] = specific Crenkov ioni    782 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
761 // and fIntegralCerenkov[0] = mean Cerenkov lo    783 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm  in keV/cm
762                                                   784 
763 void G4PAIySection::IntegralCerenkov()            785 void G4PAIySection::IntegralCerenkov()
764 {                                                 786 {
765   G4int i, k;                                     787   G4int i, k;
766    fIntegralCerenkov[fSplineNumber] = 0;          788    fIntegralCerenkov[fSplineNumber] = 0;
767    fIntegralCerenkov[0] = 0;                      789    fIntegralCerenkov[0] = 0;
768    k = fIntervalNumber -1;                        790    k = fIntervalNumber -1;
769                                                   791 
770    for( i = fSplineNumber-1; i >= 1; i-- )        792    for( i = fSplineNumber-1; i >= 1; i-- )
771    {                                              793    {
772       if(fSplineEnergy[i] >= fEnergyInterval[k    794       if(fSplineEnergy[i] >= fEnergyInterval[k])
773       {                                           795       {
774         fIntegralCerenkov[i] = fIntegralCerenk    796         fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
775         // G4cout<<"int: i = "<<i<<"; sumC = "    797         // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
776       }                                           798       }
777       else                                        799       else
778       {                                           800       {
779         fIntegralCerenkov[i] = fIntegralCerenk    801         fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + 
780                                    SumOverBord    802                                    SumOverBordCerenkov(i+1,fEnergyInterval[k]);
781         k--;                                      803         k--;
782         // G4cout<<"bord: i = "<<i<<"; sumC =     804         // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
783       }                                           805       }
784    }                                              806    }
785                                                   807 
786 }   // end of IntegralCerenkov                    808 }   // end of IntegralCerenkov 
787                                                   809 
788 //////////////////////////////////////////////    810 ////////////////////////////////////////////////////////////////////////
789 //                                                811 //
790 // Calculation of the PAI Plasmon integral cro    812 // Calculation of the PAI Plasmon integral cross-section
791 // fIntegralPlasmon[1] = splasmon primary ioni    813 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
792 // and fIntegralPlasmon[0] = mean plasmon loss    814 // and fIntegralPlasmon[0] = mean plasmon loss per cm  in keV/cm
793                                                   815 
794 void G4PAIySection::IntegralPlasmon()             816 void G4PAIySection::IntegralPlasmon()
795 {                                                 817 {
796    fIntegralPlasmon[fSplineNumber] = 0;           818    fIntegralPlasmon[fSplineNumber] = 0;
797    fIntegralPlasmon[0] = 0;                       819    fIntegralPlasmon[0] = 0;
798    G4int k = fIntervalNumber -1;                  820    G4int k = fIntervalNumber -1;
799    for(G4int i=fSplineNumber-1;i>=1;i--)          821    for(G4int i=fSplineNumber-1;i>=1;i--)
800    {                                              822    {
801       if(fSplineEnergy[i] >= fEnergyInterval[k    823       if(fSplineEnergy[i] >= fEnergyInterval[k])
802       {                                           824       {
803         fIntegralPlasmon[i] = fIntegralPlasmon    825         fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
804       }                                           826       }
805       else                                        827       else
806       {                                           828       {
807         fIntegralPlasmon[i] = fIntegralPlasmon    829         fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + 
808                                    SumOverBord    830                                    SumOverBordPlasmon(i+1,fEnergyInterval[k]);
809         k--;                                      831         k--;
810       }                                           832       }
811    }                                              833    }
                                                   >> 834 
812 }   // end of IntegralPlasmon                     835 }   // end of IntegralPlasmon
813                                                   836 
814 //////////////////////////////////////////////    837 //////////////////////////////////////////////////////////////////////
815 //                                                838 //
816 // Calculation the PAI integral cross-section     839 // Calculation the PAI integral cross-section inside
817 // of interval of continuous values of photo-i    840 // of interval of continuous values of photo-ionisation
818 // cross-section. Parameter  'i' is the number    841 // cross-section. Parameter  'i' is the number of interval.
819                                                   842 
820 G4double G4PAIySection::SumOverInterval( G4int    843 G4double G4PAIySection::SumOverInterval( G4int i )
821 {                                                 844 {         
822    G4double x0,x1,y0,yy1,a,b,c,result;            845    G4double x0,x1,y0,yy1,a,b,c,result;
823                                                   846 
824    x0 = fSplineEnergy[i];                         847    x0 = fSplineEnergy[i];
825    x1 = fSplineEnergy[i+1];                       848    x1 = fSplineEnergy[i+1];
826                                                   849 
827    if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)    850    if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
828                                                   851 
829    y0 = fDifPAIySection[i];                       852    y0 = fDifPAIySection[i];
830    yy1 = fDifPAIySection[i+1];                    853    yy1 = fDifPAIySection[i+1];
831    //G4cout << "## x0= " << x0 << " x1= " << x    854    //G4cout << "## x0= " << x0 << " x1= " << x1 << G4endl;
832    c = x1/x0;                                     855    c = x1/x0;
833    //G4cout << "c= " << c << " y0= " << y0 <<     856    //G4cout << "c= " << c << " y0= " << y0 << " yy1= " << yy1 << G4endl;
834    a = log10(yy1/y0)/log10(c);                    857    a = log10(yy1/y0)/log10(c);
835    //G4cout << "a= " << a << G4endl;              858    //G4cout << "a= " << a << G4endl;
836                                                << 859    // b = log10(y0) - a*log10(x0);
837    b = 0.0;                                    << 860    b = y0/pow(x0,a);
838    if(a < 20.) b = y0/pow(x0,a);               << 
839                                                << 
840    a += 1;                                        861    a += 1;
841    if(a == 0)                                     862    if(a == 0) 
842    {                                              863    {
843       result = b*log(x1/x0);                      864       result = b*log(x1/x0);
844    }                                              865    }
845    else                                           866    else
846    {                                              867    {
847       result = y0*(x1*pow(c,a-1) - x0)/a;         868       result = y0*(x1*pow(c,a-1) - x0)/a;
848    }                                              869    }
849    a++;                                           870    a++;
850    if(a == 0)                                     871    if(a == 0) 
851    {                                              872    {
852       fIntegralPAIySection[0] += b*log(x1/x0);    873       fIntegralPAIySection[0] += b*log(x1/x0);
853    }                                              874    }
854    else                                           875    else
855    {                                              876    {
856       fIntegralPAIySection[0] += y0*(x1*x1*pow    877       fIntegralPAIySection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
857    }                                              878    }
858    return result;                                 879    return result;
859                                                   880 
860 } //  end of SumOverInterval                      881 } //  end of SumOverInterval
861                                                   882 
862 /////////////////////////////////                 883 /////////////////////////////////
863                                                   884 
864 G4double G4PAIySection::SumOverIntervaldEdx( G    885 G4double G4PAIySection::SumOverIntervaldEdx( G4int i )
865 {                                                 886 {         
866    G4double x0,x1,y0,yy1,a,b,c,result;            887    G4double x0,x1,y0,yy1,a,b,c,result;
867                                                   888 
868    x0 = fSplineEnergy[i];                         889    x0 = fSplineEnergy[i];
869    x1 = fSplineEnergy[i+1];                       890    x1 = fSplineEnergy[i+1];
870                                                   891 
871    if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)    892    if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
872                                                   893 
873    y0 = fDifPAIySection[i];                       894    y0 = fDifPAIySection[i];
874    yy1 = fDifPAIySection[i+1];                    895    yy1 = fDifPAIySection[i+1];
875    c = x1/x0;                                     896    c = x1/x0;
876    a = log10(yy1/y0)/log10(c);                    897    a = log10(yy1/y0)/log10(c);
877                                                << 898    // b = log10(y0) - a*log10(x0);
878    b = 0.0;                                    << 899    b = y0/pow(x0,a);
879    if(a < 20.) b = y0/pow(x0,a);               << 
880                                                << 
881    a += 2;                                        900    a += 2;
882    if(a == 0)                                     901    if(a == 0) 
883    {                                              902    {
884      result = b*log(x1/x0);                       903      result = b*log(x1/x0);
885    }                                              904    }
886    else                                           905    else
887    {                                              906    {
888      result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;    907      result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
889    }                                              908    }
890    return result;                                 909    return result;
891                                                   910 
892 } //  end of SumOverInterval                      911 } //  end of SumOverInterval
893                                                   912 
894 //////////////////////////////////////////////    913 //////////////////////////////////////////////////////////////////////
895 //                                                914 //
896 // Calculation the PAI Cerenkov integral cross    915 // Calculation the PAI Cerenkov integral cross-section inside
897 // of interval of continuous values of photo-i    916 // of interval of continuous values of photo-ionisation Cerenkov
898 // cross-section. Parameter  'i' is the number    917 // cross-section. Parameter  'i' is the number of interval.
899                                                   918 
900 G4double G4PAIySection::SumOverInterCerenkov(     919 G4double G4PAIySection::SumOverInterCerenkov( G4int i )
901 {                                                 920 {         
902    G4double x0,x1,y0,yy1,a,c,result;              921    G4double x0,x1,y0,yy1,a,c,result;
903                                                   922 
904    x0  = fSplineEnergy[i];                        923    x0  = fSplineEnergy[i];
905    x1  = fSplineEnergy[i+1];                      924    x1  = fSplineEnergy[i+1];
906                                                   925 
907    if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)    926    if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
908                                                   927 
909    y0  = fdNdxCerenkov[i];                        928    y0  = fdNdxCerenkov[i];
910    yy1 = fdNdxCerenkov[i+1];                      929    yy1 = fdNdxCerenkov[i+1];
911    // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"    930    // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
912    //   <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4en    931    //   <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
913                                                   932 
914    c = x1/x0;                                     933    c = x1/x0;
915    a = log10(yy1/y0)/log10(c);                    934    a = log10(yy1/y0)/log10(c);
916    G4double b = 0.0;                              935    G4double b = 0.0;
917    if(a < 20.) b = y0/pow(x0,a);                  936    if(a < 20.) b = y0/pow(x0,a);
918                                                   937 
919    a += 1.0;                                      938    a += 1.0;
920    if(a == 0) result = b*log(c);                  939    if(a == 0) result = b*log(c);
921    else       result = y0*(x1*pow(c,a-1) - x0)    940    else       result = y0*(x1*pow(c,a-1) - x0)/a;   
922    a += 1.0;                                      941    a += 1.0;
923                                                   942 
924    if( a == 0 ) fIntegralCerenkov[0] += b*log(    943    if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
925    else         fIntegralCerenkov[0] += y0*(x1    944    else         fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
926    //  G4cout<<"a = "<<a<<"; b = "<<b<<"; resu    945    //  G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;   
927    return result;                                 946    return result;
928                                                   947 
929 } //  end of SumOverInterCerenkov                 948 } //  end of SumOverInterCerenkov
930                                                   949 
931 //////////////////////////////////////////////    950 //////////////////////////////////////////////////////////////////////
932 //                                                951 //
933 // Calculation the PAI Plasmon integral cross-    952 // Calculation the PAI Plasmon integral cross-section inside
934 // of interval of continuous values of photo-i    953 // of interval of continuous values of photo-ionisation Plasmon
935 // cross-section. Parameter  'i' is the number    954 // cross-section. Parameter  'i' is the number of interval.
936                                                   955 
937 G4double G4PAIySection::SumOverInterPlasmon( G    956 G4double G4PAIySection::SumOverInterPlasmon( G4int i )
938 {                                                 957 {         
939   G4double x0,x1,y0,yy1,a,c,result;               958   G4double x0,x1,y0,yy1,a,c,result;
940                                                   959 
941    x0  = fSplineEnergy[i];                        960    x0  = fSplineEnergy[i];
942    x1  = fSplineEnergy[i+1];                      961    x1  = fSplineEnergy[i+1];
943                                                   962 
944    if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6)    963    if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
945                                                   964 
946    y0  = fdNdxPlasmon[i];                         965    y0  = fdNdxPlasmon[i];
947    yy1 = fdNdxPlasmon[i+1];                       966    yy1 = fdNdxPlasmon[i+1];
948    c = x1/x0;                                     967    c = x1/x0;
949    a = log10(yy1/y0)/log10(c);                    968    a = log10(yy1/y0)/log10(c);
950                                                   969 
951    G4double b = 0.0;                              970    G4double b = 0.0;
952    if(a < 20.) b = y0/pow(x0,a);                  971    if(a < 20.) b = y0/pow(x0,a);
953                                                   972 
954    a += 1.0;                                      973    a += 1.0;
955    if(a == 0) result = b*log(x1/x0);              974    if(a == 0) result = b*log(x1/x0);
956    else       result = y0*(x1*pow(c,a-1) - x0)    975    else       result = y0*(x1*pow(c,a-1) - x0)/a;   
957    a += 1.0;                                      976    a += 1.0;
958                                                   977 
959    if( a == 0 ) fIntegralPlasmon[0] += b*log(x    978    if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
960    else         fIntegralPlasmon[0] += y0*(x1*    979    else         fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
961                                                   980    
962    return result;                                 981    return result;
963                                                   982 
964 } //  end of SumOverInterPlasmon                  983 } //  end of SumOverInterPlasmon
965                                                   984 
966 //////////////////////////////////////////////    985 ///////////////////////////////////////////////////////////////////////////////
967 //                                                986 //
968 // Integration of PAI cross-section for the ca    987 // Integration of PAI cross-section for the case of
969 // passing across border between intervals        988 // passing across border between intervals
970                                                   989 
971 G4double G4PAIySection::SumOverBorder( G4int      990 G4double G4PAIySection::SumOverBorder( G4int      i , 
972                                        G4doubl    991                                        G4double en0    )
973 {                                                 992 {               
974   G4double x0,x1,y0,yy1,a,d,e0,result;         << 993   G4double x0,x1,y0,yy1,a,/*c,*/d,e0,result;
975                                                   994 
976    e0 = en0;                                      995    e0 = en0;
977    x0 = fSplineEnergy[i];                         996    x0 = fSplineEnergy[i];
978    x1 = fSplineEnergy[i+1];                       997    x1 = fSplineEnergy[i+1];
979    y0 = fDifPAIySection[i];                       998    y0 = fDifPAIySection[i];
980    yy1 = fDifPAIySection[i+1];                    999    yy1 = fDifPAIySection[i+1];
981                                                   1000 
                                                   >> 1001    //c = x1/x0;
982    d = e0/x0;                                     1002    d = e0/x0;   
983    a = log10(yy1/y0)/log10(x1/x0);                1003    a = log10(yy1/y0)/log10(x1/x0);
984                                                   1004 
985    G4double b = 0.0;                              1005    G4double b = 0.0;
986    if(a < 20.) b = y0/pow(x0,a);                  1006    if(a < 20.) b = y0/pow(x0,a);
987                                                   1007    
988    a += 1;                                        1008    a += 1;
989    if(a == 0)                                     1009    if(a == 0)
990    {                                              1010    {
991       result = b*log(x0/e0);                      1011       result = b*log(x0/e0);
992    }                                              1012    }
993    else                                           1013    else
994    {                                              1014    {
995       result = y0*(x0 - e0*pow(d,a-1))/a;         1015       result = y0*(x0 - e0*pow(d,a-1))/a;
996    }                                              1016    }
997    a++;                                           1017    a++;
998    if(a == 0)                                     1018    if(a == 0)
999    {                                              1019    {
1000       fIntegralPAIySection[0] += b*log(x0/e0)    1020       fIntegralPAIySection[0] += b*log(x0/e0);
1001    }                                             1021    }
1002    else                                          1022    else 
1003    {                                             1023    {
1004       fIntegralPAIySection[0] += y0*(x0*x0 -     1024       fIntegralPAIySection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1005    }                                             1025    }
1006    x0 = fSplineEnergy[i - 1];                    1026    x0 = fSplineEnergy[i - 1];
1007    x1 = fSplineEnergy[i - 2];                    1027    x1 = fSplineEnergy[i - 2];
1008    y0 = fDifPAIySection[i - 1];                  1028    y0 = fDifPAIySection[i - 1];
1009    yy1 = fDifPAIySection[i - 2];                 1029    yy1 = fDifPAIySection[i - 2];
1010                                                  1030 
1011    //c = x1/x0;                                  1031    //c = x1/x0;
1012    d = e0/x0;                                    1032    d = e0/x0;   
1013    a = log10(yy1/y0)/log10(x1/x0);               1033    a = log10(yy1/y0)/log10(x1/x0);
1014                                               << 1034    //  b0 = log10(y0) - a*log10(x0);
1015    b = 0.0;                                   << 1035    b = y0/pow(x0,a);
1016    if(a < 20.) b = y0/pow(x0,a);              << 
1017                                               << 
1018    a += 1;                                       1036    a += 1;
1019    if(a == 0)                                    1037    if(a == 0)
1020    {                                             1038    {
1021       result += b*log(e0/x0);                    1039       result += b*log(e0/x0);
1022    }                                             1040    }
1023    else                                          1041    else
1024    {                                             1042    {
1025       result += y0*(e0*pow(d,a-1) - x0)/a;       1043       result += y0*(e0*pow(d,a-1) - x0)/a;
1026    }                                             1044    }
1027    a++;                                          1045    a++;
1028    if(a == 0)                                    1046    if(a == 0) 
1029    {                                             1047    {
1030       fIntegralPAIySection[0] += b*log(e0/x0)    1048       fIntegralPAIySection[0] += b*log(e0/x0);
1031    }                                             1049    }
1032    else                                          1050    else
1033    {                                             1051    {
1034       fIntegralPAIySection[0] += y0*(e0*e0*po    1052       fIntegralPAIySection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1035    }                                             1053    }
1036    return result;                                1054    return result;
1037                                                  1055 
1038 }                                                1056 } 
1039                                                  1057 
1040 /////////////////////////////////////////////    1058 ///////////////////////////////////////////////////////////////////////
1041                                                  1059 
1042 G4double G4PAIySection::SumOverBorderdEdx( G4    1060 G4double G4PAIySection::SumOverBorderdEdx( G4int      i , 
1043                                        G4doub    1061                                        G4double en0    )
1044 {                                                1062 {               
1045   G4double x0,x1,y0,yy1,a,/*c,*/d,e0,result;     1063   G4double x0,x1,y0,yy1,a,/*c,*/d,e0,result;
1046                                                  1064 
1047    e0 = en0;                                     1065    e0 = en0;
1048    x0 = fSplineEnergy[i];                        1066    x0 = fSplineEnergy[i];
1049    x1 = fSplineEnergy[i+1];                      1067    x1 = fSplineEnergy[i+1];
1050    y0 = fDifPAIySection[i];                      1068    y0 = fDifPAIySection[i];
1051    yy1 = fDifPAIySection[i+1];                   1069    yy1 = fDifPAIySection[i+1];
1052                                                  1070 
                                                   >> 1071    //c = x1/x0;
1053    d = e0/x0;                                    1072    d = e0/x0;   
1054    a = log10(yy1/y0)/log10(x1/x0);               1073    a = log10(yy1/y0)/log10(x1/x0);
1055                                                  1074    
1056    G4double b = 0.0;                             1075    G4double b = 0.0;
1057    if(a < 20.) b = y0/pow(x0,a);                 1076    if(a < 20.) b = y0/pow(x0,a);
1058                                                  1077    
1059    a += 2;                                       1078    a += 2;
1060    if(a == 0)                                    1079    if(a == 0)
1061    {                                             1080    {
1062       result = b*log(x0/e0);                     1081       result = b*log(x0/e0);
1063    }                                             1082    }
1064    else                                          1083    else 
1065    {                                             1084    {
1066       result = y0*(x0*x0 - e0*e0*pow(d,a-2))/    1085       result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1067    }                                             1086    }
1068    x0 = fSplineEnergy[i - 1];                    1087    x0 = fSplineEnergy[i - 1];
1069    x1 = fSplineEnergy[i - 2];                    1088    x1 = fSplineEnergy[i - 2];
1070    y0 = fDifPAIySection[i - 1];                  1089    y0 = fDifPAIySection[i - 1];
1071    yy1 = fDifPAIySection[i - 2];                 1090    yy1 = fDifPAIySection[i - 2];
1072                                                  1091 
                                                   >> 1092    //c = x1/x0;
1073    d = e0/x0;                                    1093    d = e0/x0;   
1074    a = log10(yy1/y0)/log10(x1/x0);               1094    a = log10(yy1/y0)/log10(x1/x0);
1075                                                  1095 
1076    b = 0.0;                                   << 
1077    if(a < 20.) b = y0/pow(x0,a);                 1096    if(a < 20.) b = y0/pow(x0,a);
1078                                                  1097 
1079    a += 2;                                       1098    a += 2;
1080    if(a == 0)                                    1099    if(a == 0) 
1081    {                                             1100    {
1082       result += b*log(e0/x0);                    1101       result += b*log(e0/x0);
1083    }                                             1102    }
1084    else                                          1103    else
1085    {                                             1104    {
1086       result += y0*(e0*e0*pow(d,a-2) - x0*x0)    1105       result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1087    }                                             1106    }
1088    return result;                                1107    return result;
                                                   >> 1108 
1089 }                                                1109 } 
1090                                                  1110 
1091 /////////////////////////////////////////////    1111 ///////////////////////////////////////////////////////////////////////////////
1092 //                                               1112 //
1093 // Integration of Cerenkov cross-section for     1113 // Integration of Cerenkov cross-section for the case of
1094 // passing across border between intervals       1114 // passing across border between intervals
1095                                                  1115 
1096 G4double G4PAIySection::SumOverBordCerenkov(     1116 G4double G4PAIySection::SumOverBordCerenkov( G4int      i , 
1097                                                  1117                                              G4double en0    )
1098 {                                                1118 {               
1099    G4double x0,x1,y0,yy1,a,e0,c,d,result;        1119    G4double x0,x1,y0,yy1,a,e0,c,d,result;
1100                                                  1120 
1101    e0 = en0;                                     1121    e0 = en0;
1102    x0 = fSplineEnergy[i];                        1122    x0 = fSplineEnergy[i];
1103    x1 = fSplineEnergy[i+1];                      1123    x1 = fSplineEnergy[i+1];
1104    y0 = fdNdxCerenkov[i];                        1124    y0 = fdNdxCerenkov[i];
1105    yy1 = fdNdxCerenkov[i+1];                     1125    yy1 = fdNdxCerenkov[i+1];
1106                                                  1126 
1107    //  G4cout<<G4endl;                           1127    //  G4cout<<G4endl;
1108    //G4cout<<"SumBordC, i = "<<i<<"; en0 = "<    1128    //G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1109    //     <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G    1129    //     <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1110    c = x1/x0;                                    1130    c = x1/x0;
1111    d = e0/x0;                                    1131    d = e0/x0;
1112    a = log10(yy1/y0)/log10(c);                   1132    a = log10(yy1/y0)/log10(c);
1113                                                  1133  
1114    G4double b = 0.0;                             1134    G4double b = 0.0;
1115    if(a < 20.) b = y0/pow(x0,a);                 1135    if(a < 20.) b = y0/pow(x0,a);
1116                                                  1136    
1117    a += 1.0;                                     1137    a += 1.0;
1118    if( a == 0 ) result = b*log(x0/e0);           1138    if( a == 0 ) result = b*log(x0/e0);
1119    else         result = y0*(x0 - e0*pow(d,a-    1139    else         result = y0*(x0 - e0*pow(d,a-1))/a;   
1120    a += 1.0;                                     1140    a += 1.0;
1121                                                  1141 
1122    if( a == 0 ) fIntegralCerenkov[0] += b*log    1142    if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
1123    else         fIntegralCerenkov[0] += y0*(x    1143    else         fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1124                                                  1144 
1125    //G4cout<<"a = "<<a<<"; b = "<<b<<"; resul    1145    //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1126                                                  1146    
1127    x0  = fSplineEnergy[i - 1];                   1147    x0  = fSplineEnergy[i - 1];
1128    x1  = fSplineEnergy[i - 2];                   1148    x1  = fSplineEnergy[i - 2];
1129    y0  = fdNdxCerenkov[i - 1];                   1149    y0  = fdNdxCerenkov[i - 1];
1130    yy1 = fdNdxCerenkov[i - 2];                   1150    yy1 = fdNdxCerenkov[i - 2];
1131                                                  1151 
1132    //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1           1152    //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1133    //    <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4    1153    //    <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1134                                                  1154 
1135    c = x1/x0;                                    1155    c = x1/x0;
1136    d = e0/x0;                                    1156    d = e0/x0;
1137    a  = log10(yy1/y0)/log10(x1/x0);              1157    a  = log10(yy1/y0)/log10(x1/x0);
1138                                                  1158   
1139    //   G4cout << "a= " << a << G4endl;          1159    //   G4cout << "a= " << a << G4endl;
1140    if(a > 20.0) b = 0.0;                         1160    if(a > 20.0) b = 0.0;
1141    else         b = y0/pow(x0,a);                1161    else         b = y0/pow(x0,a); 
1142                                                  1162 
1143    //G4cout << "b= " << b << G4endl;             1163    //G4cout << "b= " << b << G4endl;
1144                                                  1164 
1145    a += 1.0;                                     1165    a += 1.0;
1146    if( a == 0 ) result += b*log(e0/x0);          1166    if( a == 0 ) result += b*log(e0/x0);
1147    else         result += y0*(e0*pow(d,a-1) -    1167    else         result += y0*(e0*pow(d,a-1) - x0 )/a;
1148    a += 1.0;                                     1168    a += 1.0;
1149    //G4cout << "result= " << result << G4endl    1169    //G4cout << "result= " << result << G4endl;
1150                                                  1170 
1151    if( a == 0 )   fIntegralCerenkov[0] += b*l    1171    if( a == 0 )   fIntegralCerenkov[0] += b*log(e0/x0);
1152    else           fIntegralCerenkov[0] += y0*    1172    else           fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1153                                                  1173 
1154    //G4cout<<"a = "<<a<<"; b = "<<b<<"; resul    1174    //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;    
1155                                                  1175 
1156    return result;                                1176    return result;
                                                   >> 1177 
1157 }                                                1178 } 
1158                                                  1179 
1159 /////////////////////////////////////////////    1180 ///////////////////////////////////////////////////////////////////////////////
1160 //                                               1181 //
1161 // Integration of Plasmon cross-section for t    1182 // Integration of Plasmon cross-section for the case of
1162 // passing across border between intervals       1183 // passing across border between intervals
1163                                                  1184 
1164 G4double G4PAIySection::SumOverBordPlasmon( G    1185 G4double G4PAIySection::SumOverBordPlasmon( G4int      i , 
1165                                                  1186                                              G4double en0    )
1166 {                                                1187 {               
1167    G4double x0,x1,y0,yy1,a,c,d,e0,result;        1188    G4double x0,x1,y0,yy1,a,c,d,e0,result;
1168                                                  1189 
1169    e0 = en0;                                     1190    e0 = en0;
1170    x0 = fSplineEnergy[i];                        1191    x0 = fSplineEnergy[i];
1171    x1 = fSplineEnergy[i+1];                      1192    x1 = fSplineEnergy[i+1];
1172    y0 = fdNdxPlasmon[i];                         1193    y0 = fdNdxPlasmon[i];
1173    yy1 = fdNdxPlasmon[i+1];                      1194    yy1 = fdNdxPlasmon[i+1];
1174                                                  1195 
1175    c = x1/x0;                                    1196    c = x1/x0;
1176    d = e0/x0;                                    1197    d = e0/x0;   
1177    a = log10(yy1/y0)/log10(c);                   1198    a = log10(yy1/y0)/log10(c);
1178                                                  1199 
1179    G4double b = 0.0;                             1200    G4double b = 0.0;
1180    if(a < 20.) b = y0/pow(x0,a);                 1201    if(a < 20.) b = y0/pow(x0,a);
1181                                                  1202    
1182    a += 1.0;                                     1203    a += 1.0;
1183    if( a == 0 ) result = b*log(x0/e0);           1204    if( a == 0 ) result = b*log(x0/e0);
1184    else         result = y0*(x0 - e0*pow(d,a-    1205    else         result = y0*(x0 - e0*pow(d,a-1))/a;   
1185    a += 1.0;                                     1206    a += 1.0;
1186                                                  1207 
1187    if( a == 0 ) fIntegralPlasmon[0] += b*log(    1208    if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
1188    else         fIntegralPlasmon[0] += y0*(x0    1209    else         fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1189                                                  1210    
1190    x0 = fSplineEnergy[i - 1];                    1211    x0 = fSplineEnergy[i - 1];
1191    x1 = fSplineEnergy[i - 2];                    1212    x1 = fSplineEnergy[i - 2];
1192    y0 = fdNdxPlasmon[i - 1];                     1213    y0 = fdNdxPlasmon[i - 1];
1193    yy1 = fdNdxPlasmon[i - 2];                    1214    yy1 = fdNdxPlasmon[i - 2];
1194                                                  1215 
1195    c = x1/x0;                                    1216    c = x1/x0;
1196    d = e0/x0;                                    1217    d = e0/x0;
1197    a = log10(yy1/y0)/log10(c);                   1218    a = log10(yy1/y0)/log10(c);
1198                                                  1219  
1199    if(a < 20.) b = y0/pow(x0,a);                 1220    if(a < 20.) b = y0/pow(x0,a);
1200                                                  1221 
1201    a += 1.0;                                     1222    a += 1.0;
1202    if( a == 0 ) result += b*log(e0/x0);          1223    if( a == 0 ) result += b*log(e0/x0);
1203    else         result += y0*(e0*pow(d,a-1) -    1224    else         result += y0*(e0*pow(d,a-1) - x0)/a;
1204    a += 1.0;                                     1225    a += 1.0;
1205                                                  1226 
1206    if( a == 0 )   fIntegralPlasmon[0] += b*lo    1227    if( a == 0 )   fIntegralPlasmon[0] += b*log(e0/x0);
1207    else           fIntegralPlasmon[0] += y0*(    1228    else           fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1208                                                  1229    
1209    return result;                                1230    return result;
1210                                                  1231 
1211 }                                                1232 } 
1212                                                  1233 
1213 /////////////////////////////////////////////    1234 /////////////////////////////////////////////////////////////////////////
1214 //                                               1235 //
1215 //                                               1236 //
1216                                                  1237 
1217 G4double G4PAIySection::GetStepEnergyLoss( G4    1238 G4double G4PAIySection::GetStepEnergyLoss( G4double step )
1218 {                                                1239 {  
1219   G4int iTransfer ;                              1240   G4int iTransfer ;
1220   G4long numOfCollisions;                        1241   G4long numOfCollisions;
1221   G4double loss = 0.0;                           1242   G4double loss = 0.0;
1222   G4double meanNumber, position;                 1243   G4double meanNumber, position;
1223                                                  1244 
1224   // G4cout<<" G4PAIySection::GetStepEnergyLo    1245   // G4cout<<" G4PAIySection::GetStepEnergyLoss "<<G4endl;
1225                                                  1246 
1226                                                  1247 
1227                                                  1248 
1228   meanNumber = fIntegralPAIySection[1]*step;     1249   meanNumber = fIntegralPAIySection[1]*step;
1229   numOfCollisions = G4Poisson(meanNumber);       1250   numOfCollisions = G4Poisson(meanNumber);
1230                                                  1251 
1231   //   G4cout<<"numOfCollisions = "<<numOfCol    1252   //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1232                                                  1253 
1233   while(numOfCollisions)                         1254   while(numOfCollisions)
1234   {                                              1255   {
1235     position = fIntegralPAIySection[1]*G4Unif    1256     position = fIntegralPAIySection[1]*G4UniformRand();
1236                                                  1257 
1237     for( iTransfer=1; iTransfer<=fSplineNumbe    1258     for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1238     {                                            1259     {
1239         if( position >= fIntegralPAIySection[    1260         if( position >= fIntegralPAIySection[iTransfer] ) break;
1240     }                                            1261     }
1241     loss += fSplineEnergy[iTransfer] ;           1262     loss += fSplineEnergy[iTransfer] ;
1242     numOfCollisions--;                           1263     numOfCollisions--;
1243     // Loop checking, 03-Aug-2015, Vladimir I    1264     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1244   }                                              1265   }
1245   // G4cout<<"PAI energy loss = "<<loss/keV<<    1266   // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl; 
1246                                                  1267 
1247   return loss;                                   1268   return loss;
1248 }                                                1269 }
1249                                                  1270 
1250 /////////////////////////////////////////////    1271 /////////////////////////////////////////////////////////////////////////
1251 //                                               1272 //
1252 //                                               1273 //
1253                                                  1274 
1254 G4double G4PAIySection::GetStepCerenkovLoss(     1275 G4double G4PAIySection::GetStepCerenkovLoss( G4double step )
1255 {                                                1276 {  
1256   G4int iTransfer ;                              1277   G4int iTransfer ;
1257   G4long numOfCollisions;                        1278   G4long numOfCollisions;
1258   G4double loss = 0.0;                           1279   G4double loss = 0.0;
1259   G4double meanNumber, position;                 1280   G4double meanNumber, position;
1260                                                  1281 
1261   // G4cout<<" G4PAIySection::GetStepCreLosnk    1282   // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl;
1262                                                  1283 
1263                                                  1284 
1264                                                  1285 
1265   meanNumber = fIntegralCerenkov[1]*step;        1286   meanNumber = fIntegralCerenkov[1]*step;
1266   numOfCollisions = G4Poisson(meanNumber);       1287   numOfCollisions = G4Poisson(meanNumber);
1267                                                  1288 
1268   //   G4cout<<"numOfCollisions = "<<numOfCol    1289   //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1269                                                  1290 
1270   while(numOfCollisions)                         1291   while(numOfCollisions)
1271   {                                              1292   {
1272     position = fIntegralCerenkov[1]*G4Uniform    1293     position = fIntegralCerenkov[1]*G4UniformRand();
1273                                                  1294 
1274     for( iTransfer=1; iTransfer<=fSplineNumbe    1295     for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1275     {                                            1296     {
1276         if( position >= fIntegralCerenkov[iTr    1297         if( position >= fIntegralCerenkov[iTransfer] ) break;
1277     }                                            1298     }
1278     loss += fSplineEnergy[iTransfer] ;           1299     loss += fSplineEnergy[iTransfer] ;
1279     numOfCollisions--;                           1300     numOfCollisions--;
1280     // Loop checking, 03-Aug-2015, Vladimir I    1301     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1281   }                                              1302   }
1282   // G4cout<<"PAI Cerenkov loss = "<<loss/keV    1303   // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl; 
1283                                                  1304 
1284   return loss;                                   1305   return loss;
1285 }                                                1306 }
1286                                                  1307 
1287 /////////////////////////////////////////////    1308 /////////////////////////////////////////////////////////////////////////
1288 //                                               1309 //
1289 //                                               1310 //
1290                                                  1311 
1291 G4double G4PAIySection::GetStepPlasmonLoss( G    1312 G4double G4PAIySection::GetStepPlasmonLoss( G4double step )
1292 {                                                1313 {  
1293   G4int iTransfer ;                              1314   G4int iTransfer ;
1294   G4long numOfCollisions;                        1315   G4long numOfCollisions;
1295   G4double loss = 0.0;                           1316   G4double loss = 0.0;
1296   G4double meanNumber, position;                 1317   G4double meanNumber, position;
1297                                                  1318 
1298   // G4cout<<" G4PAIySection::GetStepCreLosnk    1319   // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl;
1299                                                  1320 
1300                                                  1321 
1301                                                  1322 
1302   meanNumber = fIntegralPlasmon[1]*step;         1323   meanNumber = fIntegralPlasmon[1]*step;
1303   numOfCollisions = G4Poisson(meanNumber);       1324   numOfCollisions = G4Poisson(meanNumber);
1304                                                  1325 
1305   //   G4cout<<"numOfCollisions = "<<numOfCol    1326   //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1306                                                  1327 
1307   while(numOfCollisions)                         1328   while(numOfCollisions)
1308   {                                              1329   {
1309     position = fIntegralPlasmon[1]*G4UniformR    1330     position = fIntegralPlasmon[1]*G4UniformRand();
1310                                                  1331 
1311     for( iTransfer=1; iTransfer<=fSplineNumbe    1332     for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1312     {                                            1333     {
1313         if( position >= fIntegralPlasmon[iTra    1334         if( position >= fIntegralPlasmon[iTransfer] ) break;
1314     }                                            1335     }
1315     loss += fSplineEnergy[iTransfer] ;           1336     loss += fSplineEnergy[iTransfer] ;
1316     numOfCollisions--;                           1337     numOfCollisions--;
1317     // Loop checking, 03-Aug-2015, Vladimir I    1338     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1318   }                                              1339   }
1319   // G4cout<<"PAI Plasmon loss = "<<loss/keV<    1340   // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl; 
1320                                                  1341 
1321   return loss;                                   1342   return loss;
1322 }                                                1343 }
1323                                                  1344 
1324 /////////////////////////////////////////////    1345 /////////////////////////////////////////////////////////////////////////////
1325 //                                               1346 //
1326                                                  1347 
1327 void G4PAIySection::CallError(G4int i, const     1348 void G4PAIySection::CallError(G4int i, const G4String& methodName) const
1328 {                                                1349 {
1329   G4String head = "G4PAIySection::" + methodN    1350   G4String head = "G4PAIySection::" + methodName + "()";
1330   G4ExceptionDescription ed;                     1351   G4ExceptionDescription ed;
1331   ed << "Wrong index " << i << " fSplineNumbe    1352   ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
1332   G4Exception(head,"pai001",FatalException,ed    1353   G4Exception(head,"pai001",FatalException,ed);
1333 }                                                1354 }
1334                                                  1355 
1335 /////////////////////////////////////////////    1356 /////////////////////////////////////////////////////////////////////////////
1336 //                                               1357 //
1337 // Init  array of Lorentz factors                1358 // Init  array of Lorentz factors
1338 //                                               1359 //
1339                                                  1360 
1340 G4int G4PAIySection::fNumberOfGammas = 111;      1361 G4int G4PAIySection::fNumberOfGammas = 111;
1341                                                  1362 
1342 const G4double G4PAIySection::fLorentzFactor[    1363 const G4double G4PAIySection::fLorentzFactor[112] =     // fNumberOfGammas+1
1343 {                                                1364 {
1344 0.0,                                             1365 0.0,
1345 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.1    1366 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
1346 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.2    1367 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
1347 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.4    1368 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
1348 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.9    1369 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
1349 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.7    1370 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
1350 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.2    1371 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
1351 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.2    1372 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
1352 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.2    1373 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
1353 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.3    1374 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
1354 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.2    1375 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
1355 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.9    1376 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
1356 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.4    1377 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
1357 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.7    1378 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
1358 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.2    1379 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
1359 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.8    1380 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
1360 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.8    1381 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
1361 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.4    1382 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
1362 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.5    1383 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
1363 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.2    1384 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
1364 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.3    1385 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
1365 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.3    1386 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
1366 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.2    1387 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
1367 1.065799e+05                                     1388 1.065799e+05
1368 };                                               1389 };
1369                                                  1390 
1370 /////////////////////////////////////////////    1391 ///////////////////////////////////////////////////////////////////////
1371 //                                               1392 //
1372 // The number of gamma for creation of  splin    1393 // The number of gamma for creation of  spline (near ion-min , G ~ 4 )
1373 //                                               1394 //
1374                                                  1395 
1375 const G4int G4PAIySection::fRefGammaNumber =     1396 const G4int G4PAIySection::fRefGammaNumber = 29; 
1376                                                  1397 
1377 //                                               1398 //   
1378 // end of G4PAIySection implementation file      1399 // end of G4PAIySection implementation file 
1379 //                                               1400 //
1380 /////////////////////////////////////////////    1401 ////////////////////////////////////////////////////////////////////////////
1381                                                  1402 
1382                                                  1403