Geant4 Cross Reference

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


  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 // $Id: G4InitXscPAI.cc,v 1.9 2006/06/29 19:53:00 gunter Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-02 $
                                                   >>  29 //
 27 //                                                 30 // 
 28 // G4InitXscPAI.cc -- class implementation fil     31 // G4InitXscPAI.cc -- class implementation file
 29 //                                                 32 //
 30 // GEANT 4 class implementation file               33 // GEANT 4 class implementation file
 31 //                                                 34 //
 32 // For information related to this code, pleas     35 // For information related to this code, please, contact
 33 // the Geant4 Collaboration.                       36 // the Geant4 Collaboration.
 34 //                                                 37 //
 35 // R&D: Vladimir.Grichine@cern.ch                  38 // R&D: Vladimir.Grichine@cern.ch
 36 //                                                 39 //
 37 // History:                                        40 // History:
 38 //                                                 41 //
 39                                                    42 
 40                                                    43 
 41                                                    44 
 42 #include "G4InitXscPAI.hh"                         45 #include "G4InitXscPAI.hh"
 43                                                    46 
 44 #include "globals.hh"                              47 #include "globals.hh"
 45 #include "G4PhysicalConstants.hh"              << 
 46 #include "G4SystemOfUnits.hh"                  << 
 47 #include "G4ios.hh"                                48 #include "G4ios.hh"
 48 #include "G4Poisson.hh"                            49 #include "G4Poisson.hh"
 49 #include "G4Integrator.hh"                         50 #include "G4Integrator.hh"
 50 #include "G4Material.hh"                           51 #include "G4Material.hh"
 51 #include "G4MaterialCutsCouple.hh"                 52 #include "G4MaterialCutsCouple.hh"
 52 #include "G4SandiaTable.hh"                        53 #include "G4SandiaTable.hh"
 53                                                    54 
 54                                                    55 
 55                                                    56 
 56 // Local class constants                           57 // Local class constants
 57                                                    58 
 58 const G4double G4InitXscPAI::fDelta        = 0     59 const G4double G4InitXscPAI::fDelta        = 0.005 ; // energy shift from interval border
 59 const G4int G4InitXscPAI::fPAIbin          = 1     60 const G4int G4InitXscPAI::fPAIbin          = 100 ;     // size of energy transfer vectors
 60 const G4double G4InitXscPAI::fSolidDensity = 0     61 const G4double G4InitXscPAI::fSolidDensity = 0.05*g/cm3 ; // ~gas-solid border
 61                                                    62 
 62 //////////////////////////////////////////////     63 //////////////////////////////////////////////////////////////////
 63 //                                                 64 //
 64 // Constructor                                     65 // Constructor
 65 //                                                 66 //
 66                                                    67 
 67 using namespace std;                               68 using namespace std;
 68                                                    69 
 69 G4InitXscPAI::G4InitXscPAI( const G4MaterialCu     70 G4InitXscPAI::G4InitXscPAI( const G4MaterialCutsCouple* matCC)
 70   : fPAIxscVector(nullptr),                    <<  71   : fPAIxscVector(NULL),
 71     fPAIdEdxVector(nullptr),                   <<  72     fPAIdEdxVector(NULL),
 72     fPAIphotonVector(nullptr),                 <<  73     fPAIphotonVector(NULL),
 73     fPAIelectronVector(nullptr),               <<  74     fPAIelectronVector(NULL),
 74     fChCosSqVector(nullptr),                   <<  75     fChCosSqVector(NULL),
 75     fChWidthVector(nullptr)                    <<  76     fChWidthVector(NULL)
 76 {                                                  77 {
 77   G4int i, j, matIndex;                            78   G4int i, j, matIndex;
 78                                                    79  
 79   fDensity         = matCC->GetMaterial()->Get     80   fDensity         = matCC->GetMaterial()->GetDensity();
 80   fElectronDensity = matCC->GetMaterial()->Get     81   fElectronDensity = matCC->GetMaterial()->GetElectronDensity();
 81   matIndex         = (G4int)matCC->GetMaterial <<  82   matIndex         = matCC->GetMaterial()->GetIndex();
 82                                                    83 
 83   fSandia          = new G4SandiaTable(matInde     84   fSandia          = new G4SandiaTable(matIndex);
 84   fIntervalNumber  = fSandia->GetMaxInterval()     85   fIntervalNumber  = fSandia->GetMaxInterval()-1;
 85                                                    86 
 86   fMatSandiaMatrix = new G4OrderedTable();         87   fMatSandiaMatrix = new G4OrderedTable();
 87                                                    88  
 88   for (i = 0; i < fIntervalNumber; ++i)        <<  89   for (i = 0; i < fIntervalNumber; i++)
 89   {                                                90   {
 90     fMatSandiaMatrix->push_back(new G4DataVect     91     fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
 91   }                                                92   }           
 92   for (i = 0; i < fIntervalNumber; ++i)        <<  93   for (G4int i = 0; i < fIntervalNumber; i++)
 93   {                                                94   {
 94     (*(*fMatSandiaMatrix)[i])[0] = fSandia->Ge     95     (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
 95                                                    96 
 96     for(j = 1; j < 5 ; ++j)                    <<  97     for(j = 1; j < 5 ; j++)
 97     {                                              98     {
 98       (*(*fMatSandiaMatrix)[i])[j] = fSandia->     99       (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
 99     }                                             100     }     
100   }                                               101   }
101   KillCloseIntervals();                           102   KillCloseIntervals();
102   Normalisation();                                103   Normalisation();
103   fBetaGammaSq = fTmax = 0.0;                  << 104 
104   fIntervalTmax = fCurrentInterval = 0;        << 
105 }                                                 105 }
106                                                   106 
107                                                   107 
108                                                   108 
109                                                   109 
110 //////////////////////////////////////////////    110 ////////////////////////////////////////////////////////////////////////////
111 //                                                111 //
112 // Destructor                                     112 // Destructor
113                                                   113 
114 G4InitXscPAI::~G4InitXscPAI()                     114 G4InitXscPAI::~G4InitXscPAI()
115 {                                                 115 {
116   delete fPAIxscVector;                        << 116   if(fPAIxscVector)      delete fPAIxscVector;  
117   delete fPAIdEdxVector;                       << 117   if(fPAIdEdxVector)     delete fPAIdEdxVector;  
118   delete fPAIphotonVector;                     << 118   if(fPAIphotonVector)   delete fPAIphotonVector;  
119   delete fPAIelectronVector;                   << 119   if(fPAIelectronVector) delete fPAIelectronVector;  
120   delete fChCosSqVector;                       << 120   if(fChCosSqVector)     delete fChCosSqVector;  
121   delete fChWidthVector;                       << 121   if(fChWidthVector)     delete fChWidthVector;  
122   delete fSandia;                              << 
123   delete fMatSandiaMatrix;                     << 
124 }                                                 122 }
125                                                   123 
126 //////////////////////////////////////////////    124 ////////////////////////////////////////////////////////////////////////
127 //                                                125 //
128 // Kill close intervals, recalculate fInterval    126 // Kill close intervals, recalculate fIntervalNumber
129                                                   127 
130 void G4InitXscPAI::KillCloseIntervals()           128 void G4InitXscPAI::KillCloseIntervals()
131 {                                                 129 {
132   G4int i, j, k;                                  130   G4int i, j, k;
133   G4double energy1, energy2;                      131   G4double energy1, energy2; 
134                                                   132             
135   for( i = 0 ; i < fIntervalNumber - 1 ; i++ )    133   for( i = 0 ; i < fIntervalNumber - 1 ; i++ )
136   {                                               134   {
137     energy1 = (*(*fMatSandiaMatrix)[i])[0];       135     energy1 = (*(*fMatSandiaMatrix)[i])[0];
138     energy2 = (*(*fMatSandiaMatrix)[i+1])[0];     136     energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
139                                                   137 
140     if( energy2 - energy1 > 1.5*fDelta*(energy    138     if( energy2 - energy1 > 1.5*fDelta*(energy1 + energy2) )  continue ;
141     else                                          139     else
142     {                                             140     {
143       for(j = i; j < fIntervalNumber-1; j++)      141       for(j = i; j < fIntervalNumber-1; j++)
144       {                                           142       {
145         for( k = 0; k < 5; k++ )                  143         for( k = 0; k < 5; k++ )
146         {                                         144         {
147           (*(*fMatSandiaMatrix)[j])[k] = (*(*f    145           (*(*fMatSandiaMatrix)[j])[k] = (*(*fMatSandiaMatrix)[j+1])[k];
148   }                                               146   }
149       }                                           147       }
150       fIntervalNumber-- ;                         148       fIntervalNumber-- ;
151       i-- ;                                       149       i-- ;
152     }                                             150     }
153   }                                               151   }
154                                                   152   
155 }                                                 153 }
156                                                   154 
157 //////////////////////////////////////////////    155 ////////////////////////////////////////////////////////////////////////
158 //                                                156 //
159 // Kill close intervals, recalculate fInterval    157 // Kill close intervals, recalculate fIntervalNumber
160                                                   158 
161 void G4InitXscPAI::Normalisation()                159 void G4InitXscPAI::Normalisation()
162 {                                                 160 {
163   G4int i, j;                                     161   G4int i, j;
164   G4double energy1, energy2, /*delta,*/ cof; / << 162   G4double energy1, energy2, delta, cof; // , shift;
165                                                   163 
166   energy1 = (*(*fMatSandiaMatrix)[fIntervalNum    164   energy1 = (*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
167   energy2 = 2.*(*(*fMatSandiaMatrix)[fInterval    165   energy2 = 2.*(*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
168                                                   166 
169                                                   167  
170   cof = RutherfordIntegral(fIntervalNumber-1,e    168   cof = RutherfordIntegral(fIntervalNumber-1,energy1,energy2);
171                                                   169             
172   for( i = fIntervalNumber-2; i >= 0; i-- )       170   for( i = fIntervalNumber-2; i >= 0; i-- )
173   {                                               171   {
174     energy1 = (*(*fMatSandiaMatrix)[i])[0];       172     energy1 = (*(*fMatSandiaMatrix)[i])[0];
175     energy2 = (*(*fMatSandiaMatrix)[i+1])[0];     173     energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
176                                                   174 
177     cof += RutherfordIntegral(i,energy1,energy    175     cof += RutherfordIntegral(i,energy1,energy2);
178     // G4cout<<"norm. cof = "<<cof<<G4endl;       176     // G4cout<<"norm. cof = "<<cof<<G4endl;
179   }                                               177   }
180   fNormalizationCof  = 2*pi*pi*hbarc*hbarc*fin    178   fNormalizationCof  = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
181   fNormalizationCof *= fElectronDensity;          179   fNormalizationCof *= fElectronDensity;
182   //delta = fNormalizationCof - cof;           << 180   delta = fNormalizationCof - cof;
183   fNormalizationCof /= cof;                       181   fNormalizationCof /= cof;
184   //  G4cout<<"G4InitXscPAI::fNormalizationCof    182   //  G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof
185   //    <<";  at delta ="<<delta<<G4endl ;        183   //    <<";  at delta ="<<delta<<G4endl ;
186                                                   184 
187   for (i = 0; i < fIntervalNumber; i++) // ren << 185   for (G4int i = 0; i < fIntervalNumber; i++) // renormalisation on QM sum rule
188   {                                               186   {
189     for(j = 1; j < 5 ; j++)                       187     for(j = 1; j < 5 ; j++)
190     {                                             188     {
191       (*(*fMatSandiaMatrix)[i])[j] *= fNormali    189       (*(*fMatSandiaMatrix)[i])[j] *= fNormalizationCof;
192     }                                             190     }     
193   }                                               191   }
194   /*                                              192   /* 
195   if(delta > 0) // shift the first energy inte    193   if(delta > 0) // shift the first energy interval
196   {                                               194   {
197     for(i=1;i<100;i++)                            195     for(i=1;i<100;i++)
198     {                                             196     {
199       energy1 = (1.-i/100.)*(*(*fMatSandiaMatr    197       energy1 = (1.-i/100.)*(*(*fMatSandiaMatrix)[0])[0];
200       energy2 = (*(*fMatSandiaMatrix)[0])[0];     198       energy2 = (*(*fMatSandiaMatrix)[0])[0];
201       shift   = RutherfordIntegral(0,energy1,e    199       shift   = RutherfordIntegral(0,energy1,energy2);
202       G4cout<<shift<<"\t";                        200       G4cout<<shift<<"\t";
203       if(shift >= delta) break;                   201       if(shift >= delta) break;
204     }                                             202     }
205     (*(*fMatSandiaMatrix)[0])[0] = energy1;       203     (*(*fMatSandiaMatrix)[0])[0] = energy1;
206     cof += shift;                                 204     cof += shift;
207   }                                               205   }
208   else if(delta < 0)                              206   else if(delta < 0)
209   {                                               207   {
210     for(i=1;i<100;i++)                            208     for(i=1;i<100;i++)
211     {                                             209     {
212       energy1 = (*(*fMatSandiaMatrix)[0])[0];     210       energy1 = (*(*fMatSandiaMatrix)[0])[0];
213       energy2 = (*(*fMatSandiaMatrix)[0])[0] +    211       energy2 = (*(*fMatSandiaMatrix)[0])[0] + 
214               ( (*(*fMatSandiaMatrix)[0])[0] -    212               ( (*(*fMatSandiaMatrix)[0])[0] - (*(*fMatSandiaMatrix)[0])[0] )*i/100.;
215       shift   = RutherfordIntegral(0,energy1,e    213       shift   = RutherfordIntegral(0,energy1,energy2);
216       if( shift >= std::abs(delta) ) break;       214       if( shift >= std::abs(delta) ) break;
217     }                                             215     }
218     (*(*fMatSandiaMatrix)[0])[0] = energy2;       216     (*(*fMatSandiaMatrix)[0])[0] = energy2;
219     cof -= shift;                                 217     cof -= shift;
220   }                                               218   }
221   G4cout<<G4cout<<"G4InitXscPAI::fNormalizatio    219   G4cout<<G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof/cof
222         <<";  at delta ="<<delta<<"  and i = "    220         <<";  at delta ="<<delta<<"  and i = "<<i<<G4endl ;
223  */                                               221  */ 
224 }                                                 222 }
225                                                   223 
                                                   >> 224 
                                                   >> 225 
                                                   >> 226 
                                                   >> 227 
226 //////////////////////////////////////////////    228 ////////////////////////////////////////////////////////////////////
227 //                                                229 //
228 // Integration over electrons that could be co    230 // Integration over electrons that could be considered
229 // quasi-free at energy transfer of interest      231 // quasi-free at energy transfer of interest
230                                                   232 
231 G4double G4InitXscPAI::RutherfordIntegral( G4i    233 G4double G4InitXscPAI::RutherfordIntegral( G4int k,
232                     G4double x1,                  234                     G4double x1,
233                       G4double x2   )             235                       G4double x2   )
234 {                                                 236 {
235    G4double  c1, c2, c3, a1, a2, a3, a4 ;         237    G4double  c1, c2, c3, a1, a2, a3, a4 ;
236                                                   238 
237    a1 = (*(*fMatSandiaMatrix)[k])[1];             239    a1 = (*(*fMatSandiaMatrix)[k])[1]; 
238    a2 = (*(*fMatSandiaMatrix)[k])[2];             240    a2 = (*(*fMatSandiaMatrix)[k])[2]; 
239    a3 = (*(*fMatSandiaMatrix)[k])[3];             241    a3 = (*(*fMatSandiaMatrix)[k])[3]; 
240    a4 = (*(*fMatSandiaMatrix)[k])[4];             242    a4 = (*(*fMatSandiaMatrix)[k])[4]; 
241    // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<    243    // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;   
242    c1 = (x2 - x1)/x1/x2 ;                         244    c1 = (x2 - x1)/x1/x2 ;
243    c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;         245    c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
244    c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x    246    c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
245    // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<    247    // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;   
246                                                   248    
247    return  a1*log(x2/x1) + a2*c1 + a3*c2/2 + a    249    return  a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;
248                                                   250 
249 }   // end of RutherfordIntegral                  251 }   // end of RutherfordIntegral 
250                                                   252 
251 //////////////////////////////////////////////    253 ///////////////////////////////////////////////////////////////
252 //                                                254 //
253 //  Integrate photo-absorption cross-section f    255 //  Integrate photo-absorption cross-section from I1 up to omega
254                                                   256 
255 G4double G4InitXscPAI::IntegralTerm(G4double o    257 G4double G4InitXscPAI::IntegralTerm(G4double omega)
256 {                                                 258 {
257   G4int i;                                        259   G4int i;
258   G4double energy1, energy2, result = 0.;         260   G4double energy1, energy2, result = 0.; 
259                                                   261             
260   for( i = 0; i <= fIntervalTmax; i++ )           262   for( i = 0; i <= fIntervalTmax; i++ )
261   {                                               263   {
262     if(i == fIntervalTmax)                        264     if(i == fIntervalTmax) 
263     {                                             265     {
264       energy1 = (*(*fMatSandiaMatrix)[i])[0];     266       energy1 = (*(*fMatSandiaMatrix)[i])[0];
265       result += RutherfordIntegral(i,energy1,o    267       result += RutherfordIntegral(i,energy1,omega);
266     }                                             268     }
267     else                                          269     else 
268     {                                             270     {
269       if( omega <= (*(*fMatSandiaMatrix)[i+1])    271       if( omega <= (*(*fMatSandiaMatrix)[i+1])[0])
270       {                                           272       {
271         energy1 = (*(*fMatSandiaMatrix)[i])[0]    273         energy1 = (*(*fMatSandiaMatrix)[i])[0];
272         result += RutherfordIntegral(i,energy1    274         result += RutherfordIntegral(i,energy1,omega);
273         break;                                    275         break;
274       }                                           276       }
275       else                                        277       else
276       {                                           278       {
277         energy1 = (*(*fMatSandiaMatrix)[i])[0]    279         energy1 = (*(*fMatSandiaMatrix)[i])[0];
278         energy2 = (*(*fMatSandiaMatrix)[i+1])[    280         energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
279         result += RutherfordIntegral(i,energy1    281         result += RutherfordIntegral(i,energy1,energy2);
280       }                                           282       }
281     }                                             283     }
282     // G4cout<<"IntegralTerm<<"("<<omega<<")"<    284     // G4cout<<"IntegralTerm<<"("<<omega<<")"<<" = "<<result<<G4endl;
283   }                                               285   }
284   return result;                                  286   return result;
285 }                                                 287 }
286                                                   288 
287                                                   289 
288 //////////////////////////////////////////////    290 ////////////////////////////////////////////////////////////////
289 //                                                291 //
290 // Imaginary part of dielectric constant          292 // Imaginary part of dielectric constant
291 // (G4int k - interval number, G4double en1 -     293 // (G4int k - interval number, G4double en1 - energy point)
292                                                   294 
293 G4double G4InitXscPAI::ImPartDielectricConst(     295 G4double G4InitXscPAI::ImPartDielectricConst( G4int    k ,
294                              G4double energy1     296                              G4double energy1 )
295 {                                                 297 {
296    G4double energy2,energy3,energy4,a1,a2,a3,a    298    G4double energy2,energy3,energy4,a1,a2,a3,a4,result;
297                                                   299 
298    a1 = (*(*fMatSandiaMatrix)[k])[1];             300    a1 = (*(*fMatSandiaMatrix)[k])[1]; 
299    a2 = (*(*fMatSandiaMatrix)[k])[2];             301    a2 = (*(*fMatSandiaMatrix)[k])[2]; 
300    a3 = (*(*fMatSandiaMatrix)[k])[3];             302    a3 = (*(*fMatSandiaMatrix)[k])[3]; 
301    a4 = (*(*fMatSandiaMatrix)[k])[4];             303    a4 = (*(*fMatSandiaMatrix)[k])[4]; 
302                                                   304 
303    energy2 = energy1*energy1;                     305    energy2 = energy1*energy1;
304    energy3 = energy2*energy1;                     306    energy3 = energy2*energy1;
305    energy4 = energy3*energy1;                     307    energy4 = energy3*energy1;
306                                                   308    
307    result  = a1/energy1+a2/energy2+a3/energy3+    309    result  = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ;  
308    result *= hbarc/energy1 ;                      310    result *= hbarc/energy1 ;
309                                                   311    
310    return result ;                                312    return result ;
311                                                   313 
312 }  // end of ImPartDielectricConst                314 }  // end of ImPartDielectricConst 
313                                                   315 
314 //////////////////////////////////////////////    316 ////////////////////////////////////////////////////////////////
315 //                                                317 //
316 // Modulus squared of dielectric constant         318 // Modulus squared of dielectric constant
317 // (G4int k - interval number, G4double omega     319 // (G4int k - interval number, G4double omega - energy point)
318                                                   320 
319 G4double G4InitXscPAI::ModuleSqDielectricConst    321 G4double G4InitXscPAI::ModuleSqDielectricConst( G4int    k ,
320                              G4double omega )     322                              G4double omega )
321 {                                                 323 {
322    G4double eIm2, eRe2, result;                   324    G4double eIm2, eRe2, result;
323                                                   325 
324    result = ImPartDielectricConst(k,omega);       326    result = ImPartDielectricConst(k,omega);
325    eIm2   = result*result;                        327    eIm2   = result*result;
326                                                   328 
327    result = RePartDielectricConst(omega);         329    result = RePartDielectricConst(omega);
328    eRe2   = result*result;                        330    eRe2   = result*result;
329                                                   331 
330    result = eIm2 + eRe2;                          332    result = eIm2 + eRe2;
331                                                   333   
332    return result ;                                334    return result ;
333 }                                                 335 }  
334                                                   336 
335                                                   337 
336 //////////////////////////////////////////////    338 //////////////////////////////////////////////////////////////////////////////
337 //                                                339 //
338 // Real part of dielectric constant minus unit    340 // Real part of dielectric constant minus unit: epsilon_1 - 1
339 // (G4double enb - energy point)                  341 // (G4double enb - energy point)
340 //                                                342 //
341                                                   343 
342 G4double G4InitXscPAI::RePartDielectricConst(G    344 G4double G4InitXscPAI::RePartDielectricConst(G4double enb)
343 {                                                 345 {
344   G4int i;                                        346   G4int i;       
345    G4double x0, x02, x03, x04, x05, x1, x2, a1    347    G4double x0, x02, x03, x04, x05, x1, x2, a1,a2,a3,a4,xx1 ,xx2 , xx12,
346             c1, c2, c3, cof1, cof2, xln1, xln2    348             c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
347                                                   349 
348    x0 = enb ;                                     350    x0 = enb ;
349    result = 0 ;                                   351    result = 0 ;
350                                                   352    
351    for( i = 0; i < fIntervalNumber-1; i++)        353    for( i = 0; i < fIntervalNumber-1; i++)
352    {                                              354    {
353       x1 = (*(*fMatSandiaMatrix)[i])[0];          355       x1 = (*(*fMatSandiaMatrix)[i])[0];
354       x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;       356       x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;
355                                                   357 
356       a1 = (*(*fMatSandiaMatrix)[i])[1];          358       a1 = (*(*fMatSandiaMatrix)[i])[1]; 
357       a2 = (*(*fMatSandiaMatrix)[i])[2];          359       a2 = (*(*fMatSandiaMatrix)[i])[2]; 
358       a3 = (*(*fMatSandiaMatrix)[i])[3];          360       a3 = (*(*fMatSandiaMatrix)[i])[3]; 
359       a4 = (*(*fMatSandiaMatrix)[i])[4];          361       a4 = (*(*fMatSandiaMatrix)[i])[4];
360                                                   362  
361       if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta    363       if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta ) 
362       {                                           364       {
363         if(x0 >= x1) x0 = x1*(1+fDelta);          365         if(x0 >= x1) x0 = x1*(1+fDelta);
364         else         x0 = x1*(1-fDelta);          366         else         x0 = x1*(1-fDelta);
365       }                                           367       } 
366       if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta    368       if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta ) 
367       {                                           369       {
368         if(x0 >= x2) x0 = x2*(1+fDelta);          370         if(x0 >= x2) x0 = x2*(1+fDelta);
369         else         x0 = x2*(1-fDelta);          371         else         x0 = x2*(1-fDelta);
370       }                                           372       }
371       xx1 = x1 - x0 ;                             373       xx1 = x1 - x0 ;
372       xx2 = x2 - x0 ;                             374       xx2 = x2 - x0 ;
373       xx12 = xx2/xx1 ;                            375       xx12 = xx2/xx1 ;
374                                                   376       
375       if( xx12 < 0 ) xx12 = -xx12;                377       if( xx12 < 0 ) xx12 = -xx12;
376                                                   378       
377       xln1 = log(x2/x1) ;                         379       xln1 = log(x2/x1) ;
378       xln2 = log(xx12) ;                          380       xln2 = log(xx12) ;
379       xln3 = log((x2 + x0)/(x1 + x0)) ;           381       xln3 = log((x2 + x0)/(x1 + x0)) ;
380                                                   382 
381       x02 = x0*x0 ;                               383       x02 = x0*x0 ;
382       x03 = x02*x0 ;                              384       x03 = x02*x0 ;
383       x04 = x03*x0 ;                              385       x04 = x03*x0 ;
384       x05 = x04*x0;                               386       x05 = x04*x0;
385                                                   387 
386       c1  = (x2 - x1)/x1/x2 ;                     388       c1  = (x2 - x1)/x1/x2 ;
387       c2  = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;      389       c2  = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
388       c3  = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x    390       c3  = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
389                                                   391 
390       result -= (a1/x02 + a3/x04)*xln1 ;          392       result -= (a1/x02 + a3/x04)*xln1 ;
391       result -= (a2/x02 + a4/x04)*c1 ;            393       result -= (a2/x02 + a4/x04)*c1 ;
392       result -= a3*c2/2/x02 ;                     394       result -= a3*c2/2/x02 ;
393       result -= a4*c3/3/x02 ;                     395       result -= a4*c3/3/x02 ;
394                                                   396 
395       cof1 = a1/x02 + a3/x04 ;                    397       cof1 = a1/x02 + a3/x04 ;
396       cof2 = a2/x03 + a4/x05 ;                    398       cof2 = a2/x03 + a4/x05 ;
397                                                   399 
398       result += 0.5*(cof1 +cof2)*xln2 ;           400       result += 0.5*(cof1 +cof2)*xln2 ;
399       result += 0.5*(cof1 - cof2)*xln3 ;          401       result += 0.5*(cof1 - cof2)*xln3 ;
400    }                                              402    } 
401    result *= 2*hbarc/pi ;                         403    result *= 2*hbarc/pi ;
402                                                   404    
403    return result ;                                405    return result ;
404                                                   406 
405 }   // end of RePartDielectricConst               407 }   // end of RePartDielectricConst 
406                                                   408 
407 //////////////////////////////////////////////    409 //////////////////////////////////////////////////////////////////////
408 //                                                410 //
409 // PAI differential cross-section in terms of     411 // PAI differential cross-section in terms of
410 // simplified Allison's equation                  412 // simplified Allison's equation
411 //                                                413 //
412                                                   414 
413 G4double G4InitXscPAI::DifPAIxSection( G4doubl    415 G4double G4InitXscPAI::DifPAIxSection( G4double omega )
414 {                                                 416 {
415   G4int i = fCurrentInterval;                     417   G4int i = fCurrentInterval;
416   G4double  betaGammaSq = fBetaGammaSq;           418   G4double  betaGammaSq = fBetaGammaSq;       
417   G4double integralTerm = IntegralTerm(omega);    419   G4double integralTerm = IntegralTerm(omega);
418   G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,res    420   G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
419   G4double epsilonRe = RePartDielectricConst(o    421   G4double epsilonRe = RePartDielectricConst(omega);
420   G4double epsilonIm = ImPartDielectricConst(i    422   G4double epsilonIm = ImPartDielectricConst(i,omega);
421   G4double be4 ;                                  423   G4double be4 ;
422   static const G4double betaBohr2 = fine_struc << 424   G4double betaBohr2 = fine_structure_const*fine_structure_const ;
423   static const G4double betaBohr4 = betaBohr2* << 425   G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
424   be2 = betaGammaSq/(1 + betaGammaSq) ;           426   be2 = betaGammaSq/(1 + betaGammaSq) ;
425   be4 = be2*be2 ;                                 427   be4 = be2*be2 ;
426                                                   428  
427    cof = 1 ;                                      429    cof = 1 ;
428    x1 = log(2*electron_mass_c2/omega) ;           430    x1 = log(2*electron_mass_c2/omega) ;
429                                                   431 
430    if( betaGammaSq < 0.01 ) x2 = log(be2) ;       432    if( betaGammaSq < 0.01 ) x2 = log(be2) ;
431    else                                           433    else
432    {                                              434    {
433      x2 = -log( (1/betaGammaSq - epsilonRe)*      435      x2 = -log( (1/betaGammaSq - epsilonRe)*
434           (1/betaGammaSq - epsilonRe) +           436           (1/betaGammaSq - epsilonRe) + 
435           epsilonIm*epsilonIm )/2 ;               437           epsilonIm*epsilonIm )/2 ;
436    }                                              438    }
437    if( epsilonIm == 0.0 || betaGammaSq < 0.01     439    if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
438    {                                              440    {
439      x6=0 ;                                       441      x6=0 ;
440    }                                              442    }
441    else                                           443    else
442    {                                              444    {
443      x3 = -epsilonRe + 1/betaGammaSq ;            445      x3 = -epsilonRe + 1/betaGammaSq ;
444      x5 = -1 - epsilonRe + be2*((1 +epsilonRe)    446      x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) +
445      epsilonIm*epsilonIm) ;                       447      epsilonIm*epsilonIm) ;
446                                                   448 
447      x7 = atan2(epsilonIm,x3) ;                   449      x7 = atan2(epsilonIm,x3) ;
448      x6 = x5 * x7 ;                               450      x6 = x5 * x7 ;
449    }                                              451    }
450     // if(fImPartDielectricConst[i] == 0) x6 =    452     // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
451                                                   453    
452    x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ;        454    x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ;
453    //   if( x4 < 0.0 ) x4 = 0.0 ;                 455    //   if( x4 < 0.0 ) x4 = 0.0 ;
454    x8 = (1 + epsilonRe)*(1 + epsilonRe) +         456    x8 = (1 + epsilonRe)*(1 + epsilonRe) + 
455         epsilonIm*epsilonIm;                      457         epsilonIm*epsilonIm;
456                                                   458 
457    result = (x4 + cof*integralTerm/omega/omega    459    result = (x4 + cof*integralTerm/omega/omega) ;
458    if(result < 1.0e-8) result = 1.0e-8 ;          460    if(result < 1.0e-8) result = 1.0e-8 ;
459    result *= fine_structure_const/be2/pi ;        461    result *= fine_structure_const/be2/pi ;
460    //   result *= (1-exp(-beta/betaBohr))*(1-e    462    //   result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
461    //  result *= (1-exp(-be2/betaBohr2)) ;        463    //  result *= (1-exp(-be2/betaBohr2)) ;
462    result *= (1-exp(-be4/betaBohr4)) ;            464    result *= (1-exp(-be4/betaBohr4)) ;
463    if(fDensity >= fSolidDensity)                  465    if(fDensity >= fSolidDensity)
464    {                                              466    { 
465       result /= x8 ;                              467       result /= x8 ;
466    }                                              468    }
467    return result ;                                469    return result ;
468                                                   470 
469 } // end of DifPAIxSection                        471 } // end of DifPAIxSection 
470                                                   472 
471 //////////////////////////////////////////////    473 //////////////////////////////////////////////////////////////////////
472 //                                                474 //
473 // Differential PAI dEdx(omega)=omega*dNdx(ome    475 // Differential PAI dEdx(omega)=omega*dNdx(omega)
474 //                                                476 //
475                                                   477 
476 G4double G4InitXscPAI::DifPAIdEdx( G4double om    478 G4double G4InitXscPAI::DifPAIdEdx( G4double omega )
477 {                                                 479 {
478   G4double dEdx = omega*DifPAIxSection(omega);    480   G4double dEdx = omega*DifPAIxSection(omega);
479   return dEdx;                                    481   return dEdx;
480 }                                                 482 }
481                                                   483 
482 //////////////////////////////////////////////    484 //////////////////////////////////////////////////////////////////////////
483 //                                                485 //
484 // Calculation od dN/dx of collisions with cre    486 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
485                                                   487 
486 G4double G4InitXscPAI::PAIdNdxCherenkov( G4dou    488 G4double G4InitXscPAI::PAIdNdxCherenkov( G4double omega  )
487 {                                                 489 {        
488   G4int i = fCurrentInterval;                     490   G4int i = fCurrentInterval;
489   G4double  betaGammaSq = fBetaGammaSq;           491   G4double  betaGammaSq = fBetaGammaSq;       
490   G4double epsilonRe = RePartDielectricConst(o    492   G4double epsilonRe = RePartDielectricConst(omega);
491   G4double epsilonIm = ImPartDielectricConst(i    493   G4double epsilonIm = ImPartDielectricConst(i,omega);
492                                                   494 
493   G4double /*cof,*/ logarithm, x3, x5, argumen << 495   G4double cof, logarithm, x3, x5, argument, modul2, dNdxC ; 
494   G4double be2, be4;                           << 496   G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr ;
495                                                   497 
496   //cof         = 1.0 ;                        << 498    cof         = 1.0 ;
497   static const G4double cofBetaBohr = 4.0 ;    << 499    cofBetaBohr = 4.0 ;
498   static const G4double betaBohr2   = fine_str << 500    betaBohr2   = fine_structure_const*fine_structure_const ;
499   static const G4double betaBohr4   = betaBohr << 501    betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr ;
500                                                   502 
501    be2 = betaGammaSq/(1 + betaGammaSq) ;          503    be2 = betaGammaSq/(1 + betaGammaSq) ;
502    be4 = be2*be2 ;                                504    be4 = be2*be2 ;
503                                                   505 
504    if( betaGammaSq < 0.01 ) logarithm = log(1.    506    if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
505    else                                           507    else
506    {                                              508    {
507      logarithm  = -log( (1/betaGammaSq - epsil    509      logarithm  = -log( (1/betaGammaSq - epsilonRe)*
508                   (1/betaGammaSq - epsilonRe)     510                   (1/betaGammaSq - epsilonRe) + 
509                   epsilonIm*epsilonIm )*0.5 ;     511                   epsilonIm*epsilonIm )*0.5 ;
510      logarithm += log(1+1.0/betaGammaSq) ;        512      logarithm += log(1+1.0/betaGammaSq) ;
511    }                                              513    }
512                                                   514 
513    if( epsilonIm == 0.0 || betaGammaSq < 0.01     515    if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
514    {                                              516    {
515      argument = 0.0 ;                             517      argument = 0.0 ;
516    }                                              518    }
517    else                                           519    else
518    {                                              520    {
519      x3 = -epsilonRe + 1.0/betaGammaSq ;          521      x3 = -epsilonRe + 1.0/betaGammaSq ;
520      x5 = -1.0 - epsilonRe +                      522      x5 = -1.0 - epsilonRe +
521           be2*((1.0 +epsilonRe)*(1.0 + epsilon    523           be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) +
522     epsilonIm*epsilonIm) ;                        524     epsilonIm*epsilonIm) ;
523      if( x3 == 0.0 ) argument = 0.5*pi;           525      if( x3 == 0.0 ) argument = 0.5*pi;
524      else            argument = atan2(epsilonI    526      else            argument = atan2(epsilonIm,x3) ;
525      argument *= x5  ;                            527      argument *= x5  ;
526    }                                              528    }   
527    dNdxC = ( logarithm*epsilonIm + argument )/    529    dNdxC = ( logarithm*epsilonIm + argument )/hbarc ;
528                                                   530   
529    if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;            531    if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
530                                                   532 
531    dNdxC *= fine_structure_const/be2/pi ;         533    dNdxC *= fine_structure_const/be2/pi ;
532                                                   534 
533    dNdxC *= (1-exp(-be4/betaBohr4)) ;             535    dNdxC *= (1-exp(-be4/betaBohr4)) ;
534                                                   536 
535    if(fDensity >= fSolidDensity)                  537    if(fDensity >= fSolidDensity)
536    {                                              538    { 
537       modul2 = (1.0 + epsilonRe)*(1.0 + epsilo    539       modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) + 
538                     epsilonIm*epsilonIm;          540                     epsilonIm*epsilonIm;
539       dNdxC /= modul2 ;                           541       dNdxC /= modul2 ;
540    }                                              542    }
541    return dNdxC ;                                 543    return dNdxC ;
542                                                   544 
543 } // end of PAIdNdxCerenkov                       545 } // end of PAIdNdxCerenkov 
544                                                   546 
545 //////////////////////////////////////////////    547 //////////////////////////////////////////////////////////////////////////
546 //                                                548 //
547 // Calculation od dN/dx of collisions with cre    549 // Calculation od dN/dx of collisions with creation of longitudinal EM
548 // excitations (plasmons, delta-electrons)        550 // excitations (plasmons, delta-electrons)
549                                                   551 
550 G4double G4InitXscPAI::PAIdNdxPlasmon( G4doubl    552 G4double G4InitXscPAI::PAIdNdxPlasmon( G4double omega )
551 {                                                 553 {        
552   G4int i = fCurrentInterval;                     554   G4int i = fCurrentInterval;
553   G4double  betaGammaSq = fBetaGammaSq;           555   G4double  betaGammaSq = fBetaGammaSq;       
554   G4double integralTerm = IntegralTerm(omega);    556   G4double integralTerm = IntegralTerm(omega);
555   G4double epsilonRe = RePartDielectricConst(o    557   G4double epsilonRe = RePartDielectricConst(omega);
556   G4double epsilonIm = ImPartDielectricConst(i    558   G4double epsilonIm = ImPartDielectricConst(i,omega);
557                                                   559 
558    G4double cof, resonance, modul2, dNdxP ;       560    G4double cof, resonance, modul2, dNdxP ;
559    G4double be2, be4;                          << 561    G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr ;
560                                                   562 
561    cof = 1 ;                                      563    cof = 1 ;
562    static const G4double cofBetaBohr = 4.0 ;   << 564    cofBetaBohr = 4.0 ;
563    static const G4double betaBohr2   = fine_st << 565    betaBohr2   = fine_structure_const*fine_structure_const ;
564    static const G4double betaBohr4   = betaBoh << 566    betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr ;
565                                                   567 
566    be2 = betaGammaSq/(1 + betaGammaSq) ;          568    be2 = betaGammaSq/(1 + betaGammaSq) ;
567    be4 = be2*be2 ;                                569    be4 = be2*be2 ;
568                                                   570  
569    resonance  = log(2*electron_mass_c2*be2/ome    571    resonance  = log(2*electron_mass_c2*be2/omega) ;  
570    resonance *= epsilonIm/hbarc ;                 572    resonance *= epsilonIm/hbarc ;
571                                                   573 
572                                                   574 
573    dNdxP = ( resonance + cof*integralTerm/omeg    575    dNdxP = ( resonance + cof*integralTerm/omega/omega ) ;
574                                                   576 
575    if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;          577    if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
576                                                   578 
577    dNdxP *= fine_structure_const/be2/pi ;         579    dNdxP *= fine_structure_const/be2/pi ;
578    dNdxP *= (1-exp(-be4/betaBohr4)) ;             580    dNdxP *= (1-exp(-be4/betaBohr4)) ;
579                                                   581 
580    if( fDensity >= fSolidDensity )                582    if( fDensity >= fSolidDensity )
581    {                                              583    { 
582      modul2 = (1 + epsilonRe)*(1 + epsilonRe)     584      modul2 = (1 + epsilonRe)*(1 + epsilonRe) + 
583                epsilonIm*epsilonIm;               585                epsilonIm*epsilonIm;
584      dNdxP /= modul2 ;                            586      dNdxP /= modul2 ;
585    }                                              587    }
586    return dNdxP ;                                 588    return dNdxP ;
587                                                   589 
588 } // end of PAIdNdxPlasmon                        590 } // end of PAIdNdxPlasmon 
589                                                   591 
590 //////////////////////////////////////////////    592 ////////////////////////////////////////////////////////////////////////
591 //                                                593 //
592 // Calculation of the PAI integral cross-secti    594 // Calculation of the PAI integral cross-section
593 // = specific primary ionisation, 1/cm            595 // = specific primary ionisation, 1/cm
594 //                                                596 // 
595                                                   597 
596 void G4InitXscPAI::IntegralPAIxSection(G4doubl    598 void G4InitXscPAI::IntegralPAIxSection(G4double bg2, G4double Tmax)
597 {                                                 599 {
598   G4int i,k,i1,i2;                                600   G4int i,k,i1,i2;
599   G4double energy1, energy2, result = 0.;         601   G4double energy1, energy2, result = 0.;
600                                                   602  
601   fBetaGammaSq = bg2;                             603   fBetaGammaSq = bg2;
602   fTmax        = Tmax;                            604   fTmax        = Tmax;
603                                                   605 
604   delete fPAIxscVector;                        << 606   if(fPAIxscVector) delete fPAIxscVector;  
605                                                   607   
606   fPAIxscVector = new G4PhysicsLogVector( (*(*    608   fPAIxscVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
607   fPAIxscVector->PutValue(fPAIbin-1,result);      609   fPAIxscVector->PutValue(fPAIbin-1,result);
608                                                   610             
609   for( i = fIntervalNumber - 1; i >= 0; i-- )     611   for( i = fIntervalNumber - 1; i >= 0; i-- )
610   {                                               612   {
611     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] )    613     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
612   }                                               614   }
613   if (i < 0) i = 0; // Tmax should be more tha    615   if (i < 0) i = 0; // Tmax should be more than 
614                     // first ionisation potent    616                     // first ionisation potential
615   fIntervalTmax = i;                              617   fIntervalTmax = i;
616                                                   618 
617   G4Integrator<G4InitXscPAI,G4double(G4InitXsc    619   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
618                                                   620 
619   for( k = fPAIbin - 2; k >= 0; k-- )             621   for( k = fPAIbin - 2; k >= 0; k-- )
620   {                                               622   {
621     energy1 = fPAIxscVector->GetLowEdgeEnergy(    623     energy1 = fPAIxscVector->GetLowEdgeEnergy(k);
622     energy2 = fPAIxscVector->GetLowEdgeEnergy(    624     energy2 = fPAIxscVector->GetLowEdgeEnergy(k+1);
623                                                   625 
624     for( i = fIntervalTmax; i >= 0; i-- )         626     for( i = fIntervalTmax; i >= 0; i-- ) 
625     {                                             627     {
626       if( energy2 > (*(*fMatSandiaMatrix)[i])[    628       if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
627     }                                             629     }
628     if(i < 0) i = 0;                              630     if(i < 0) i = 0;
629     i2 = i;                                       631     i2 = i;
630                                                   632 
631     for( i = fIntervalTmax; i >= 0; i-- )         633     for( i = fIntervalTmax; i >= 0; i-- ) 
632     {                                             634     {
633       if( energy1 > (*(*fMatSandiaMatrix)[i])[    635       if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
634     }                                             636     }
635     if(i < 0) i = 0;                              637     if(i < 0) i = 0;
636     i1 = i;                                       638     i1 = i;
637                                                   639 
638     if( i1 == i2 )                                640     if( i1 == i2 )
639     {                                             641     {
640       fCurrentInterval = i1;                      642       fCurrentInterval = i1;
641       result += integral.Legendre10(this,&G4In    643       result += integral.Legendre10(this,&G4InitXscPAI::DifPAIxSection,
642                                     energy1,en    644                                     energy1,energy2);
643       fPAIxscVector->PutValue(k,result);          645       fPAIxscVector->PutValue(k,result);
644     }                                             646     }
645     else                                          647     else
646     {                                             648     {
647       for( i = i2; i >= i1; i-- )                 649       for( i = i2; i >= i1; i-- ) 
648       {                                           650       {
649         fCurrentInterval = i;                     651         fCurrentInterval = i;
650                                                   652 
651         if( i==i2 )        result += integral.    653         if( i==i2 )        result += integral.Legendre10(this,
652                            &G4InitXscPAI::DifP    654                            &G4InitXscPAI::DifPAIxSection,
653                            (*(*fMatSandiaMatri    655                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
654                                                   656 
655   else if( i == i1 ) result += integral.Legend    657   else if( i == i1 ) result += integral.Legendre10(this,
656                            &G4InitXscPAI::DifP    658                            &G4InitXscPAI::DifPAIxSection,energy1,
657                            (*(*fMatSandiaMatri    659                            (*(*fMatSandiaMatrix)[i+1])[0]);
658                                                   660 
659         else               result += integral.    661         else               result += integral.Legendre10(this,
660                            &G4InitXscPAI::DifP    662                            &G4InitXscPAI::DifPAIxSection,
661                        (*(*fMatSandiaMatrix)[i    663                        (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
662       }                                           664       }
663       fPAIxscVector->PutValue(k,result);          665       fPAIxscVector->PutValue(k,result);
664     }                                             666     }
665     // G4cout<<k<<"\t"<<result<<G4endl;           667     // G4cout<<k<<"\t"<<result<<G4endl; 
666   }                                               668   }
667   return ;                                        669   return ;
668 }                                                 670 }
669                                                   671 
670                                                   672 
671 //////////////////////////////////////////////    673 ////////////////////////////////////////////////////////////////////////
672 //                                                674 //
673 // Calculation of the PAI integral dEdx           675 // Calculation of the PAI integral dEdx
674 // = mean energy loss per unit length, keV/cm     676 // = mean energy loss per unit length, keV/cm
675 //                                                677 // 
676                                                   678 
677 void G4InitXscPAI::IntegralPAIdEdx(G4double bg    679 void G4InitXscPAI::IntegralPAIdEdx(G4double bg2, G4double Tmax)
678 {                                                 680 {
679   G4int i,k,i1,i2;                                681   G4int i,k,i1,i2;
680   G4double energy1, energy2, result = 0.;         682   G4double energy1, energy2, result = 0.;
681                                                   683  
682   fBetaGammaSq = bg2;                             684   fBetaGammaSq = bg2;
683   fTmax        = Tmax;                            685   fTmax        = Tmax;
684                                                   686 
685   delete fPAIdEdxVector;                       << 687   if(fPAIdEdxVector) delete fPAIdEdxVector;  
686                                                   688   
687   fPAIdEdxVector = new G4PhysicsLogVector( (*(    689   fPAIdEdxVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
688   fPAIdEdxVector->PutValue(fPAIbin-1,result);     690   fPAIdEdxVector->PutValue(fPAIbin-1,result);
689                                                   691             
690   for( i = fIntervalNumber - 1; i >= 0; i-- )     692   for( i = fIntervalNumber - 1; i >= 0; i-- )
691   {                                               693   {
692     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] )    694     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
693   }                                               695   }
694   if (i < 0) i = 0; // Tmax should be more tha    696   if (i < 0) i = 0; // Tmax should be more than 
695                     // first ionisation potent    697                     // first ionisation potential
696   fIntervalTmax = i;                              698   fIntervalTmax = i;
697                                                   699 
698   G4Integrator<G4InitXscPAI,G4double(G4InitXsc    700   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
699                                                   701 
700   for( k = fPAIbin - 2; k >= 0; k-- )             702   for( k = fPAIbin - 2; k >= 0; k-- )
701   {                                               703   {
702     energy1 = fPAIdEdxVector->GetLowEdgeEnergy    704     energy1 = fPAIdEdxVector->GetLowEdgeEnergy(k);
703     energy2 = fPAIdEdxVector->GetLowEdgeEnergy    705     energy2 = fPAIdEdxVector->GetLowEdgeEnergy(k+1);
704                                                   706 
705     for( i = fIntervalTmax; i >= 0; i-- )         707     for( i = fIntervalTmax; i >= 0; i-- ) 
706     {                                             708     {
707       if( energy2 > (*(*fMatSandiaMatrix)[i])[    709       if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
708     }                                             710     }
709     if(i < 0) i = 0;                              711     if(i < 0) i = 0;
710     i2 = i;                                       712     i2 = i;
711                                                   713 
712     for( i = fIntervalTmax; i >= 0; i-- )         714     for( i = fIntervalTmax; i >= 0; i-- ) 
713     {                                             715     {
714       if( energy1 > (*(*fMatSandiaMatrix)[i])[    716       if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
715     }                                             717     }
716     if(i < 0) i = 0;                              718     if(i < 0) i = 0;
717     i1 = i;                                       719     i1 = i;
718                                                   720 
719     if( i1 == i2 )                                721     if( i1 == i2 )
720     {                                             722     {
721       fCurrentInterval = i1;                      723       fCurrentInterval = i1;
722       result += integral.Legendre10(this,&G4In    724       result += integral.Legendre10(this,&G4InitXscPAI::DifPAIdEdx,
723                                     energy1,en    725                                     energy1,energy2);
724       fPAIdEdxVector->PutValue(k,result);         726       fPAIdEdxVector->PutValue(k,result);
725     }                                             727     }
726     else                                          728     else
727     {                                             729     {
728       for( i = i2; i >= i1; i-- )                 730       for( i = i2; i >= i1; i-- ) 
729       {                                           731       {
730         fCurrentInterval = i;                     732         fCurrentInterval = i;
731                                                   733 
732         if( i==i2 )        result += integral.    734         if( i==i2 )        result += integral.Legendre10(this,
733                            &G4InitXscPAI::DifP    735                            &G4InitXscPAI::DifPAIdEdx,
734                            (*(*fMatSandiaMatri    736                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
735                                                   737 
736   else if( i == i1 ) result += integral.Legend    738   else if( i == i1 ) result += integral.Legendre10(this,
737                            &G4InitXscPAI::DifP    739                            &G4InitXscPAI::DifPAIdEdx,energy1,
738                            (*(*fMatSandiaMatri    740                            (*(*fMatSandiaMatrix)[i+1])[0]);
739                                                   741 
740         else               result += integral.    742         else               result += integral.Legendre10(this,
741                            &G4InitXscPAI::DifP    743                            &G4InitXscPAI::DifPAIdEdx,
742                        (*(*fMatSandiaMatrix)[i    744                        (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
743       }                                           745       }
744       fPAIdEdxVector->PutValue(k,result);         746       fPAIdEdxVector->PutValue(k,result);
745     }                                             747     }
746     // G4cout<<k<<"\t"<<result<<G4endl;           748     // G4cout<<k<<"\t"<<result<<G4endl; 
747   }                                               749   }
748   return ;                                        750   return ;
749 }                                                 751 }
750                                                   752 
751 //////////////////////////////////////////////    753 ////////////////////////////////////////////////////////////////////////
752 //                                                754 //
753 // Calculation of the PAI Cerenkov integral cr    755 // Calculation of the PAI Cerenkov integral cross-section
754 // fIntegralCrenkov[1] = specific Crenkov ioni    756 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
755 // and fIntegralCerenkov[0] = mean Cerenkov lo    757 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm  in keV/cm
756                                                   758 
757 void G4InitXscPAI::IntegralCherenkov(G4double     759 void G4InitXscPAI::IntegralCherenkov(G4double bg2, G4double Tmax)
758 {                                                 760 {
759   G4int i,k,i1,i2;                                761   G4int i,k,i1,i2;
760   G4double energy1, energy2, beta2, module2, c    762   G4double energy1, energy2, beta2, module2, cos2, width, result = 0.;
761                                                   763  
762   fBetaGammaSq = bg2;                             764   fBetaGammaSq = bg2;
763   fTmax        = Tmax;                            765   fTmax        = Tmax;
764   beta2        = bg2/(1+bg2);                     766   beta2        = bg2/(1+bg2);
765                                                   767 
766   delete fPAIphotonVector;                     << 768   if(fPAIphotonVector) delete fPAIphotonVector;  
767   delete fChCosSqVector;                       << 769   if(fChCosSqVector)   delete fChCosSqVector;  
768   delete fChWidthVector;                       << 770   if(fChWidthVector)   delete fChWidthVector;  
769                                                   771   
770   fPAIphotonVector = new G4PhysicsLogVector( (    772   fPAIphotonVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
771   fChCosSqVector = new G4PhysicsLogVector( (*(    773   fChCosSqVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
772   fChWidthVector = new G4PhysicsLogVector( (*(    774   fChWidthVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
773                                                   775 
774   fPAIphotonVector->PutValue(fPAIbin-1,result)    776   fPAIphotonVector->PutValue(fPAIbin-1,result);
775   fChCosSqVector->PutValue(fPAIbin-1,1.);         777   fChCosSqVector->PutValue(fPAIbin-1,1.);
776   fChWidthVector->PutValue(fPAIbin-1,1e-7);       778   fChWidthVector->PutValue(fPAIbin-1,1e-7);
777                                                   779             
778   for( i = fIntervalNumber - 1; i >= 0; i-- )     780   for( i = fIntervalNumber - 1; i >= 0; i-- )
779   {                                               781   {
780     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] )    782     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
781   }                                               783   }
782   if (i < 0) i = 0; // Tmax should be more tha    784   if (i < 0) i = 0; // Tmax should be more than 
783                     // first ionisation potent    785                     // first ionisation potential
784   fIntervalTmax = i;                              786   fIntervalTmax = i;
785                                                   787 
786   G4Integrator<G4InitXscPAI,G4double(G4InitXsc    788   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
787                                                   789 
788   for( k = fPAIbin - 2; k >= 0; k-- )             790   for( k = fPAIbin - 2; k >= 0; k-- )
789   {                                               791   {
790     energy1 = fPAIphotonVector->GetLowEdgeEner    792     energy1 = fPAIphotonVector->GetLowEdgeEnergy(k);
791     energy2 = fPAIphotonVector->GetLowEdgeEner    793     energy2 = fPAIphotonVector->GetLowEdgeEnergy(k+1);
792                                                   794 
793     for( i = fIntervalTmax; i >= 0; i-- )         795     for( i = fIntervalTmax; i >= 0; i-- ) 
794     {                                             796     {
795       if( energy2 > (*(*fMatSandiaMatrix)[i])[    797       if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
796     }                                             798     }
797     if(i < 0) i = 0;                              799     if(i < 0) i = 0;
798     i2 = i;                                       800     i2 = i;
799                                                   801 
800     for( i = fIntervalTmax; i >= 0; i-- )         802     for( i = fIntervalTmax; i >= 0; i-- ) 
801     {                                             803     {
802       if( energy1 > (*(*fMatSandiaMatrix)[i])[    804       if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
803     }                                             805     }
804     if(i < 0) i = 0;                              806     if(i < 0) i = 0;
805     i1 = i;                                       807     i1 = i;
806                                                   808 
807     module2 = ModuleSqDielectricConst(i1,energ    809     module2 = ModuleSqDielectricConst(i1,energy1);
808     cos2    = RePartDielectricConst(energy1)/m    810     cos2    = RePartDielectricConst(energy1)/module2/beta2;      
809     width   = ImPartDielectricConst(i1,energy1    811     width   = ImPartDielectricConst(i1,energy1)/module2/beta2;
810                                                   812       
811     fChCosSqVector->PutValue(k,cos2);             813     fChCosSqVector->PutValue(k,cos2);
812     fChWidthVector->PutValue(k,width);            814     fChWidthVector->PutValue(k,width);
813                                                   815 
814     if( i1 == i2 )                                816     if( i1 == i2 )
815     {                                             817     {
816       fCurrentInterval = i1;                      818       fCurrentInterval = i1;
817       result += integral.Legendre10(this,&G4In    819       result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxCherenkov,
818                                     energy1,en    820                                     energy1,energy2);
819       fPAIphotonVector->PutValue(k,result);       821       fPAIphotonVector->PutValue(k,result);
820                                                   822 
821     }                                             823     }
822     else                                          824     else
823     {                                             825     {
824       for( i = i2; i >= i1; i-- )                 826       for( i = i2; i >= i1; i-- ) 
825       {                                           827       {
826         fCurrentInterval = i;                     828         fCurrentInterval = i;
827                                                   829 
828         if( i==i2 )        result += integral.    830         if( i==i2 )        result += integral.Legendre10(this,
829                            &G4InitXscPAI::PAId    831                            &G4InitXscPAI::PAIdNdxCherenkov,
830                            (*(*fMatSandiaMatri    832                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
831                                                   833 
832   else if( i == i1 ) result += integral.Legend    834   else if( i == i1 ) result += integral.Legendre10(this,
833                            &G4InitXscPAI::PAId    835                            &G4InitXscPAI::PAIdNdxCherenkov,energy1,
834                            (*(*fMatSandiaMatri    836                            (*(*fMatSandiaMatrix)[i+1])[0]);
835                                                   837 
836         else               result += integral.    838         else               result += integral.Legendre10(this,
837                            &G4InitXscPAI::PAId    839                            &G4InitXscPAI::PAIdNdxCherenkov,
838                        (*(*fMatSandiaMatrix)[i    840                        (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
839       }                                           841       }
840       fPAIphotonVector->PutValue(k,result);       842       fPAIphotonVector->PutValue(k,result);
841     }                                             843     }
842     // G4cout<<k<<"\t"<<result<<G4endl;           844     // G4cout<<k<<"\t"<<result<<G4endl; 
843   }                                               845   }
844   return;                                         846   return;
845 }   // end of IntegralCerenkov                    847 }   // end of IntegralCerenkov 
846                                                   848 
847 //////////////////////////////////////////////    849 ////////////////////////////////////////////////////////////////////////
848 //                                                850 //
849 // Calculation of the PAI Plasmon integral cro    851 // Calculation of the PAI Plasmon integral cross-section
850 // fIntegralPlasmon[1] = splasmon primary ioni    852 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
851 // and fIntegralPlasmon[0] = mean plasmon loss    853 // and fIntegralPlasmon[0] = mean plasmon loss per cm  in keV/cm
852                                                   854 
853 void G4InitXscPAI::IntegralPlasmon(G4double bg    855 void G4InitXscPAI::IntegralPlasmon(G4double bg2, G4double Tmax)
854 {                                                 856 {
855   G4int i,k,i1,i2;                                857   G4int i,k,i1,i2;
856   G4double energy1, energy2, result = 0.;         858   G4double energy1, energy2, result = 0.;
857                                                   859  
858   fBetaGammaSq = bg2;                             860   fBetaGammaSq = bg2;
859   fTmax        = Tmax;                            861   fTmax        = Tmax;
860                                                   862 
861   delete fPAIelectronVector;                   << 863   if(fPAIelectronVector) delete fPAIelectronVector;  
862                                                   864   
863   fPAIelectronVector = new G4PhysicsLogVector(    865   fPAIelectronVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
864   fPAIelectronVector->PutValue(fPAIbin-1,resul    866   fPAIelectronVector->PutValue(fPAIbin-1,result);
865                                                   867             
866   for( i = fIntervalNumber - 1; i >= 0; i-- )     868   for( i = fIntervalNumber - 1; i >= 0; i-- )
867   {                                               869   {
868     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] )    870     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
869   }                                               871   }
870   if (i < 0) i = 0; // Tmax should be more tha    872   if (i < 0) i = 0; // Tmax should be more than 
871                     // first ionisation potent    873                     // first ionisation potential
872   fIntervalTmax = i;                              874   fIntervalTmax = i;
873                                                   875 
874   G4Integrator<G4InitXscPAI,G4double(G4InitXsc    876   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
875                                                   877 
876   for( k = fPAIbin - 2; k >= 0; k-- )             878   for( k = fPAIbin - 2; k >= 0; k-- )
877   {                                               879   {
878     energy1 = fPAIelectronVector->GetLowEdgeEn    880     energy1 = fPAIelectronVector->GetLowEdgeEnergy(k);
879     energy2 = fPAIelectronVector->GetLowEdgeEn    881     energy2 = fPAIelectronVector->GetLowEdgeEnergy(k+1);
880                                                   882 
881     for( i = fIntervalTmax; i >= 0; i-- )         883     for( i = fIntervalTmax; i >= 0; i-- ) 
882     {                                             884     {
883       if( energy2 > (*(*fMatSandiaMatrix)[i])[    885       if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
884     }                                             886     }
885     if(i < 0) i = 0;                              887     if(i < 0) i = 0;
886     i2 = i;                                       888     i2 = i;
887                                                   889 
888     for( i = fIntervalTmax; i >= 0; i-- )         890     for( i = fIntervalTmax; i >= 0; i-- ) 
889     {                                             891     {
890       if( energy1 > (*(*fMatSandiaMatrix)[i])[    892       if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
891     }                                             893     }
892     if(i < 0) i = 0;                              894     if(i < 0) i = 0;
893     i1 = i;                                       895     i1 = i;
894                                                   896 
895     if( i1 == i2 )                                897     if( i1 == i2 )
896     {                                             898     {
897       fCurrentInterval = i1;                      899       fCurrentInterval = i1;
898       result += integral.Legendre10(this,&G4In    900       result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxPlasmon,
899                                     energy1,en    901                                     energy1,energy2);
900       fPAIelectronVector->PutValue(k,result);     902       fPAIelectronVector->PutValue(k,result);
901     }                                             903     }
902     else                                          904     else
903     {                                             905     {
904       for( i = i2; i >= i1; i-- )                 906       for( i = i2; i >= i1; i-- ) 
905       {                                           907       {
906         fCurrentInterval = i;                     908         fCurrentInterval = i;
907                                                   909 
908         if( i==i2 )        result += integral.    910         if( i==i2 )        result += integral.Legendre10(this,
909                            &G4InitXscPAI::PAId    911                            &G4InitXscPAI::PAIdNdxPlasmon,
910                            (*(*fMatSandiaMatri    912                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
911                                                   913 
912   else if( i == i1 ) result += integral.Legend    914   else if( i == i1 ) result += integral.Legendre10(this,
913                            &G4InitXscPAI::PAId    915                            &G4InitXscPAI::PAIdNdxPlasmon,energy1,
914                            (*(*fMatSandiaMatri    916                            (*(*fMatSandiaMatrix)[i+1])[0]);
915                                                   917 
916         else               result += integral.    918         else               result += integral.Legendre10(this,
917                            &G4InitXscPAI::PAId    919                            &G4InitXscPAI::PAIdNdxPlasmon,
918                        (*(*fMatSandiaMatrix)[i    920                        (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
919       }                                           921       }
920       fPAIelectronVector->PutValue(k,result);     922       fPAIelectronVector->PutValue(k,result);
921     }                                             923     }
922     // G4cout<<k<<"\t"<<result<<G4endl;           924     // G4cout<<k<<"\t"<<result<<G4endl; 
923   }                                               925   }
924   return;                                         926   return;
925 }   // end of IntegralPlasmon                     927 }   // end of IntegralPlasmon
926                                                   928 
927                                                   929 
928 //////////////////////////////////////////////    930 /////////////////////////////////////////////////////////////////////////
929 //                                                931 //
930 //                                                932 //
931                                                   933 
932 G4double G4InitXscPAI::GetPhotonLambda( G4doub    934 G4double G4InitXscPAI::GetPhotonLambda( G4double omega )
933 {                                                 935 {  
934   G4int i ;                                       936   G4int i ;
935   G4double omega2, omega3, omega4, a1, a2, a3,    937   G4double omega2, omega3, omega4, a1, a2, a3, a4, lambda ;
936                                                   938 
937   omega2 = omega*omega ;                          939   omega2 = omega*omega ;
938   omega3 = omega2*omega ;                         940   omega3 = omega2*omega ;
939   omega4 = omega2*omega2 ;                        941   omega4 = omega2*omega2 ;
940                                                   942 
941   for(i = 0; i < fIntervalNumber;i++)             943   for(i = 0; i < fIntervalNumber;i++)
942   {                                               944   {
943     if( omega < (*(*fMatSandiaMatrix)[i])[0] )    945     if( omega < (*(*fMatSandiaMatrix)[i])[0] ) break ;
944   }                                               946   }
945   if( i == 0 )                                    947   if( i == 0 )
946   {                                               948   {
947     G4cout<<"Warning: energy in G4InitXscPAI::    949     G4cout<<"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<G4endl;
948   }                                               950   }
949   else i-- ;                                      951   else i-- ;
950                                                   952 
951   a1 = (*(*fMatSandiaMatrix)[i])[1];              953   a1 = (*(*fMatSandiaMatrix)[i])[1]; 
952   a2 = (*(*fMatSandiaMatrix)[i])[2];              954   a2 = (*(*fMatSandiaMatrix)[i])[2]; 
953   a3 = (*(*fMatSandiaMatrix)[i])[3];              955   a3 = (*(*fMatSandiaMatrix)[i])[3]; 
954   a4 = (*(*fMatSandiaMatrix)[i])[4];              956   a4 = (*(*fMatSandiaMatrix)[i])[4]; 
955                                                   957 
956   lambda = 1./(a1/omega + a2/omega2 + a3/omega    958   lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4);
957                                                   959 
958   return lambda ;                                 960   return lambda ;
959 }                                                 961 }
960                                                   962 
961 //////////////////////////////////////////////    963 /////////////////////////////////////////////////////////////////////////
962 //                                                964 //
963 //                                                965 //
964                                                   966 
965 //////////////////////////////////////////////    967 /////////////////////////////////////////////////////////////////////////
966 //                                                968 //
967 //                                                969 //
968                                                   970 
969 G4double G4InitXscPAI::GetStepEnergyLoss( G4do    971 G4double G4InitXscPAI::GetStepEnergyLoss( G4double step )
970 {                                                 972 {  
971   G4double loss = 0.0 ;                           973   G4double loss = 0.0 ;
972   loss *= step;                                   974   loss *= step;
973                                                   975 
974   return loss ;                                   976   return loss ;
975 }                                                 977 }
976                                                   978 
977 //////////////////////////////////////////////    979 /////////////////////////////////////////////////////////////////////////
978 //                                                980 //
979 //                                                981 //
980                                                   982 
981 G4double G4InitXscPAI::GetStepCerenkovLoss( G4    983 G4double G4InitXscPAI::GetStepCerenkovLoss( G4double step )
982 {                                                 984 {  
983   G4double loss = 0.0 ;                           985   G4double loss = 0.0 ;
984   loss *= step;                                   986   loss *= step;
985                                                   987 
986   return loss ;                                   988   return loss ;
987 }                                                 989 }
988                                                   990 
989 //////////////////////////////////////////////    991 /////////////////////////////////////////////////////////////////////////
990 //                                                992 //
991 //                                                993 //
992                                                   994 
993 G4double G4InitXscPAI::GetStepPlasmonLoss( G4d    995 G4double G4InitXscPAI::GetStepPlasmonLoss( G4double step )
994 {                                                 996 {  
995                                                   997 
996                                                   998 
997   G4double loss = 0.0 ;                           999   G4double loss = 0.0 ;
998   loss *= step;                                   1000   loss *= step;
999   return loss ;                                   1001   return loss ;
1000 }                                                1002 }
1001                                                  1003 
1002                                                  1004    
1003 //                                               1005 //   
1004 // end of G4InitXscPAI implementation file       1006 // end of G4InitXscPAI implementation file 
1005 //                                               1007 //
1006 /////////////////////////////////////////////    1008 ////////////////////////////////////////////////////////////////////////////
1007                                                  1009 
1008                                                  1010