Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4PAIxSection.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/G4PAIxSection.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4PAIxSection.cc (Version 11.2.2)


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