Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4eIonisationSpectrum.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/lowenergy/src/G4eIonisationSpectrum.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4eIonisationSpectrum.cc (Version 10.4.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 // $Id: G4eIonisationSpectrum.cc 95950 2016-03-03 10:42:48Z gcosmo $
 26 //                                                 27 //
 27 // -------------------------------------------     28 // -------------------------------------------------------------------
 28 //                                                 29 //
 29 // GEANT4 Class file                               30 // GEANT4 Class file
 30 //                                                 31 //
 31 //                                                 32 //
 32 // File name:     G4eIonisationSpectrum            33 // File name:     G4eIonisationSpectrum
 33 //                                                 34 //
 34 // Author:        V.Ivanchenko (Vladimir.Ivanc     35 // Author:        V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
 35 //                                                 36 //
 36 // Creation date: 29 September 2001                37 // Creation date: 29 September 2001
 37 //                                                 38 //
 38 // Modifications:                                  39 // Modifications:
 39 // 10.10.2001 MGP           Revision to improv     40 // 10.10.2001 MGP           Revision to improve code quality and
 40 //                          consistency with d     41 //                          consistency with design
 41 // 02.11.2001 VI            Optimize sampling      42 // 02.11.2001 VI            Optimize sampling of energy
 42 // 29.11.2001 VI            New parametrisatio     43 // 29.11.2001 VI            New parametrisation
 43 // 19.04.2002 VI            Add protection in      44 // 19.04.2002 VI            Add protection in case of energy below binding
 44 // 30.05.2002 VI            Update to 24-param     45 // 30.05.2002 VI            Update to 24-parameters data
 45 // 11.07.2002 VI            Fix in integration     46 // 11.07.2002 VI            Fix in integration over spectrum
 46 // 23.03.2009 LP            Added protection a     47 // 23.03.2009 LP            Added protection against division by zero (for
 47 //                          faulty database fi     48 //                          faulty database files), Bug Report 1042
 48 //                                                 49 //
 49 // -------------------------------------------     50 // -------------------------------------------------------------------
 50 //                                                 51 //
 51                                                    52 
 52 #include "G4eIonisationSpectrum.hh"                53 #include "G4eIonisationSpectrum.hh"
 53 #include "G4AtomicTransitionManager.hh"            54 #include "G4AtomicTransitionManager.hh"
 54 #include "G4AtomicShell.hh"                        55 #include "G4AtomicShell.hh"
 55 #include "G4DataVector.hh"                         56 #include "G4DataVector.hh"
 56 #include "Randomize.hh"                            57 #include "Randomize.hh"
 57 #include "G4PhysicalConstants.hh"                  58 #include "G4PhysicalConstants.hh"
 58 #include "G4SystemOfUnits.hh"                      59 #include "G4SystemOfUnits.hh"
 59 #include "G4Exp.hh"                                60 #include "G4Exp.hh"
 60                                                    61 
 61 G4eIonisationSpectrum::G4eIonisationSpectrum()     62 G4eIonisationSpectrum::G4eIonisationSpectrum():G4VEnergySpectrum(),
 62   lowestE(0.1*eV),                                 63   lowestE(0.1*eV),
 63   factor(1.3),                                     64   factor(1.3),
 64   iMax(24),                                        65   iMax(24),
 65   verbose(0)                                       66   verbose(0)
 66 {                                                  67 {
 67   theParam = new G4eIonisationParameters();        68   theParam = new G4eIonisationParameters();
 68 }                                                  69 }
 69                                                    70 
 70                                                    71 
 71 G4eIonisationSpectrum::~G4eIonisationSpectrum(     72 G4eIonisationSpectrum::~G4eIonisationSpectrum()
 72 {                                                  73 {
 73   delete theParam;                                 74   delete theParam;
 74 }                                                  75 }
 75                                                    76 
 76                                                    77 
 77 G4double G4eIonisationSpectrum::Probability(G4     78 G4double G4eIonisationSpectrum::Probability(G4int Z,
 78                     G4double tMin,                 79                     G4double tMin,
 79               G4double tMax,                       80               G4double tMax,
 80               G4double e,                          81               G4double e,
 81               G4int shell,                         82               G4int shell,
 82               const G4ParticleDefinition* ) co     83               const G4ParticleDefinition* ) const
 83 {                                                  84 {
 84   // Please comment what Probability does and      85   // Please comment what Probability does and what are the three
 85   // functions mentioned below                     86   // functions mentioned below
 86   // Describe the algorithms used                  87   // Describe the algorithms used
 87                                                    88 
 88   G4double eMax = MaxEnergyOfSecondaries(e);       89   G4double eMax = MaxEnergyOfSecondaries(e);
 89   G4double t0 = std::max(tMin, lowestE);           90   G4double t0 = std::max(tMin, lowestE);
 90   G4double tm = std::min(tMax, eMax);              91   G4double tm = std::min(tMax, eMax);
 91   if(t0 >= tm) return 0.0;                         92   if(t0 >= tm) return 0.0;
 92                                                    93 
 93   G4double bindingEnergy = (G4AtomicTransition     94   G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
 94                            Shell(Z, shell)->Bi     95                            Shell(Z, shell)->BindingEnergy();
 95                                                    96 
 96   if(e <= bindingEnergy) return 0.0;               97   if(e <= bindingEnergy) return 0.0;
 97                                                    98 
 98   G4double energy = e + bindingEnergy;             99   G4double energy = e + bindingEnergy;
 99                                                   100 
100   G4double x1 = std::min(0.5,(t0 + bindingEner    101   G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
101   G4double x2 = std::min(0.5,(tm + bindingEner    102   G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
102                                                   103 
103   if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.    104   if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
104     G4cout << "G4eIonisationSpectrum::Probabil    105     G4cout << "G4eIonisationSpectrum::Probability: Z= " << Z
105            << "; shell= " << shell                106            << "; shell= " << shell
106            << "; E(keV)= " << e/keV               107            << "; E(keV)= " << e/keV
107            << "; Eb(keV)= " << bindingEnergy/k    108            << "; Eb(keV)= " << bindingEnergy/keV
108            << "; x1= " << x1                      109            << "; x1= " << x1
109            << "; x2= " << x2                      110            << "; x2= " << x2
110            << G4endl;                             111            << G4endl;
111                                                   112 
112   }                                               113   }
113                                                   114 
114   G4DataVector p;                                 115   G4DataVector p;
115                                                   116 
116   // Access parameters                            117   // Access parameters
117   for (G4int i=0; i<iMax; i++)                    118   for (G4int i=0; i<iMax; i++)
118   {                                               119   {
119     G4double x = theParam->Parameter(Z, shell,    120     G4double x = theParam->Parameter(Z, shell, i, e);
120     if(i<4) x /= energy;                          121     if(i<4) x /= energy;
121     p.push_back(x);                               122     p.push_back(x);
122   }                                               123   }
123                                                   124 
124   if(p[3] > 0.5) p[3] = 0.5;                      125   if(p[3] > 0.5) p[3] = 0.5;
125                                                   126 
126   G4double gLocal = energy/electron_mass_c2 +     127   G4double gLocal = energy/electron_mass_c2 + 1.;
127   p.push_back((2.0*gLocal - 1.0)/(gLocal*gLoca    128   p.push_back((2.0*gLocal - 1.0)/(gLocal*gLocal));
128                                                   129 
129   //Add protection against division by zero: a    130   //Add protection against division by zero: actually in Function()
130   //parameter p[3] appears in the denominator.    131   //parameter p[3] appears in the denominator. Bug report 1042
131   if (p[3] > 0)                                   132   if (p[3] > 0)
132     p[iMax-1] = Function(p[3], p);                133     p[iMax-1] = Function(p[3], p);
133   else                                            134   else
134     {                                             135     {
135       G4cout << "WARNING: G4eIonisationSpectru    136       G4cout << "WARNING: G4eIonisationSpectrum::Probability "
136        << "parameter p[3] <= 0. G4LEDATA dabat    137        << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
137        << Z << ". Please check and/or update i    138        << Z << ". Please check and/or update it " << G4endl;
138     }                                             139     }
139                                                   140 
140   if(e >= 1. && e <= 0. && Z == 4) p.push_back    141   if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
141                                                   142 
142                                                   143 
143   G4double val = IntSpectrum(x1, x2, p);          144   G4double val = IntSpectrum(x1, x2, p);
144   G4double x0  = (lowestE + bindingEnergy)/ene    145   G4double x0  = (lowestE + bindingEnergy)/energy;
145   G4double nor = IntSpectrum(x0, 0.5, p);         146   G4double nor = IntSpectrum(x0, 0.5, p);
146                                                   147 
147   if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.    148   if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
148     G4cout << "tcut= " << tMin                    149     G4cout << "tcut= " << tMin
149            << "; tMax= " << tMax                  150            << "; tMax= " << tMax
150            << "; x0= " << x0                      151            << "; x0= " << x0
151            << "; x1= " << x1                      152            << "; x1= " << x1
152            << "; x2= " << x2                      153            << "; x2= " << x2
153            << "; val= " << val                    154            << "; val= " << val
154            << "; nor= " << nor                    155            << "; nor= " << nor
155            << "; sum= " << p[0]                   156            << "; sum= " << p[0]
156            << "; a= " << p[1]                     157            << "; a= " << p[1]
157            << "; b= " << p[2]                     158            << "; b= " << p[2]
158            << "; c= " << p[3]                     159            << "; c= " << p[3]
159            << G4endl;                             160            << G4endl;
160     if(shell == 1) G4cout << "============" <<    161     if(shell == 1) G4cout << "============" << G4endl;
161   }                                               162   }
162                                                   163 
163   p.clear();                                      164   p.clear();
164                                                   165 
165   if(nor > 0.0) val /= nor;                       166   if(nor > 0.0) val /= nor;
166   else          val  = 0.0;                       167   else          val  = 0.0;
167                                                   168 
168   return val;                                     169   return val;
169 }                                                 170 }
170                                                   171 
171                                                   172 
172 G4double G4eIonisationSpectrum::AverageEnergy(    173 G4double G4eIonisationSpectrum::AverageEnergy(G4int Z,
173                       G4double tMin,              174                       G4double tMin,
174                 G4double tMax,                    175                 G4double tMax,
175                 G4double e,                       176                 G4double e,
176                 G4int shell,                      177                 G4int shell,
177                 const G4ParticleDefinition* )     178                 const G4ParticleDefinition* ) const
178 {                                                 179 {
179   // Please comment what AverageEnergy does an    180   // Please comment what AverageEnergy does and what are the three
180   // functions mentioned below                    181   // functions mentioned below
181   // Describe the algorithms used                 182   // Describe the algorithms used
182                                                   183 
183   G4double eMax = MaxEnergyOfSecondaries(e);      184   G4double eMax = MaxEnergyOfSecondaries(e);
184   G4double t0 = std::max(tMin, lowestE);          185   G4double t0 = std::max(tMin, lowestE);
185   G4double tm = std::min(tMax, eMax);             186   G4double tm = std::min(tMax, eMax);
186   if(t0 >= tm) return 0.0;                        187   if(t0 >= tm) return 0.0;
187                                                   188 
188   G4double bindingEnergy = (G4AtomicTransition    189   G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
189                            Shell(Z, shell)->Bi    190                            Shell(Z, shell)->BindingEnergy();
190                                                   191 
191   if(e <= bindingEnergy) return 0.0;              192   if(e <= bindingEnergy) return 0.0;
192                                                   193 
193   G4double energy = e + bindingEnergy;            194   G4double energy = e + bindingEnergy;
194                                                   195 
195   G4double x1 = std::min(0.5,(t0 + bindingEner    196   G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
196   G4double x2 = std::min(0.5,(tm + bindingEner    197   G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
197                                                   198 
198   if(verbose > 1) {                               199   if(verbose > 1) {
199     G4cout << "G4eIonisationSpectrum::AverageE    200     G4cout << "G4eIonisationSpectrum::AverageEnergy: Z= " << Z
200            << "; shell= " << shell                201            << "; shell= " << shell
201            << "; E(keV)= " << e/keV               202            << "; E(keV)= " << e/keV
202            << "; bindingE(keV)= " << bindingEn    203            << "; bindingE(keV)= " << bindingEnergy/keV
203            << "; x1= " << x1                      204            << "; x1= " << x1
204            << "; x2= " << x2                      205            << "; x2= " << x2
205            << G4endl;                             206            << G4endl;
206   }                                               207   }
207                                                   208 
208   G4DataVector p;                                 209   G4DataVector p;
209                                                   210 
210   // Access parameters                            211   // Access parameters
211   for (G4int i=0; i<iMax; i++)                    212   for (G4int i=0; i<iMax; i++)
212   {                                               213   {
213     G4double x = theParam->Parameter(Z, shell,    214     G4double x = theParam->Parameter(Z, shell, i, e);
214     if(i<4) x /= energy;                          215     if(i<4) x /= energy;
215     p.push_back(x);                               216     p.push_back(x);
216   }                                               217   }
217                                                   218 
218   if(p[3] > 0.5) p[3] = 0.5;                      219   if(p[3] > 0.5) p[3] = 0.5;
219                                                   220 
220   G4double gLocal2 = energy/electron_mass_c2 +    221   G4double gLocal2 = energy/electron_mass_c2 + 1.;
221   p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLo    222   p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLocal2));
222                                                   223 
223                                                   224 
224   //Add protection against division by zero: a    225   //Add protection against division by zero: actually in Function()
225   //parameter p[3] appears in the denominator.    226   //parameter p[3] appears in the denominator. Bug report 1042
226   if (p[3] > 0)                                   227   if (p[3] > 0)
227     p[iMax-1] = Function(p[3], p);                228     p[iMax-1] = Function(p[3], p);
228   else                                            229   else
229     {                                             230     {
230       G4cout << "WARNING: G4eIonisationSpectru    231       G4cout << "WARNING: G4eIonisationSpectrum::AverageEnergy "
231        << "parameter p[3] <= 0. G4LEDATA dabat    232        << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
232        << Z << ". Please check and/or update i    233        << Z << ". Please check and/or update it " << G4endl;
233     }                                             234     }
234                                                   235 
235   G4double val = AverageValue(x1, x2, p);         236   G4double val = AverageValue(x1, x2, p);
236   G4double x0  = (lowestE + bindingEnergy)/ene    237   G4double x0  = (lowestE + bindingEnergy)/energy;
237   G4double nor = IntSpectrum(x0, 0.5, p);         238   G4double nor = IntSpectrum(x0, 0.5, p);
238   val *= energy;                                  239   val *= energy;
239                                                   240 
240   if(verbose > 1) {                               241   if(verbose > 1) {
241     G4cout << "tcut(MeV)= " << tMin/MeV           242     G4cout << "tcut(MeV)= " << tMin/MeV
242            << "; tMax(MeV)= " << tMax/MeV         243            << "; tMax(MeV)= " << tMax/MeV
243            << "; x0= " << x0                      244            << "; x0= " << x0
244            << "; x1= " << x1                      245            << "; x1= " << x1
245            << "; x2= " << x2                      246            << "; x2= " << x2
246            << "; val= " << val                    247            << "; val= " << val
247            << "; nor= " << nor                    248            << "; nor= " << nor
248            << "; sum= " << p[0]                   249            << "; sum= " << p[0]
249            << "; a= " << p[1]                     250            << "; a= " << p[1]
250            << "; b= " << p[2]                     251            << "; b= " << p[2]
251            << "; c= " << p[3]                     252            << "; c= " << p[3]
252            << G4endl;                             253            << G4endl;
253   }                                               254   }
254                                                   255 
255   p.clear();                                      256   p.clear();
256                                                   257 
257   if(nor > 0.0) val /= nor;                       258   if(nor > 0.0) val /= nor;
258   else          val  = 0.0;                       259   else          val  = 0.0;
259                                                   260 
260   return val;                                     261   return val;
261 }                                                 262 }
262                                                   263 
263                                                   264 
264 G4double G4eIonisationSpectrum::SampleEnergy(G    265 G4double G4eIonisationSpectrum::SampleEnergy(G4int Z,
265                G4double tMin,                     266                G4double tMin,
266                G4double tMax,                     267                G4double tMax,
267                      G4double e,                  268                      G4double e,
268                      G4int shell,                 269                      G4int shell,
269                const G4ParticleDefinition* ) c    270                const G4ParticleDefinition* ) const
270 {                                                 271 {
271   // Please comment what SampleEnergy does        272   // Please comment what SampleEnergy does
272   G4double tDelta = 0.0;                          273   G4double tDelta = 0.0;
273   G4double t0 = std::max(tMin, lowestE);          274   G4double t0 = std::max(tMin, lowestE);
274   G4double tm = std::min(tMax, MaxEnergyOfSeco    275   G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
275   if(t0 > tm) return tDelta;                      276   if(t0 > tm) return tDelta;
276                                                   277 
277   G4double bindingEnergy = (G4AtomicTransition    278   G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
278                            Shell(Z, shell)->Bi    279                            Shell(Z, shell)->BindingEnergy();
279                                                   280 
280   if(e <= bindingEnergy) return 0.0;              281   if(e <= bindingEnergy) return 0.0;
281                                                   282 
282   G4double energy = e + bindingEnergy;            283   G4double energy = e + bindingEnergy;
283                                                   284 
284   G4double x1 = std::min(0.5,(t0 + bindingEner    285   G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
285   G4double x2 = std::min(0.5,(tm + bindingEner    286   G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
286   if(x1 >= x2) return tDelta;                     287   if(x1 >= x2) return tDelta;
287                                                   288 
288   if(verbose > 1) {                               289   if(verbose > 1) {
289     G4cout << "G4eIonisationSpectrum::SampleEn    290     G4cout << "G4eIonisationSpectrum::SampleEnergy: Z= " << Z
290            << "; shell= " << shell                291            << "; shell= " << shell
291            << "; E(keV)= " << e/keV               292            << "; E(keV)= " << e/keV
292            << G4endl;                             293            << G4endl;
293   }                                               294   }
294                                                   295 
295   // Access parameters                            296   // Access parameters
296   G4DataVector p;                                 297   G4DataVector p;
297                                                   298 
298   // Access parameters                            299   // Access parameters
299   for (G4int i=0; i<iMax; i++)                    300   for (G4int i=0; i<iMax; i++)
300   {                                               301   {
301     G4double x = theParam->Parameter(Z, shell,    302     G4double x = theParam->Parameter(Z, shell, i, e);
302     if(i<4) x /= energy;                          303     if(i<4) x /= energy;
303     p.push_back(x);                               304     p.push_back(x);
304   }                                               305   }
305                                                   306 
306   if(p[3] > 0.5) p[3] = 0.5;                      307   if(p[3] > 0.5) p[3] = 0.5;
307                                                   308 
308   G4double gLocal3 = energy/electron_mass_c2 +    309   G4double gLocal3 = energy/electron_mass_c2 + 1.;
309   p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLo    310   p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLocal3));
310                                                   311 
311                                                   312 
312   //Add protection against division by zero: a    313   //Add protection against division by zero: actually in Function()
313   //parameter p[3] appears in the denominator.    314   //parameter p[3] appears in the denominator. Bug report 1042
314   if (p[3] > 0)                                   315   if (p[3] > 0)
315     p[iMax-1] = Function(p[3], p);                316     p[iMax-1] = Function(p[3], p);
316   else                                            317   else
317     {                                             318     {
318       G4cout << "WARNING: G4eIonisationSpectru    319       G4cout << "WARNING: G4eIonisationSpectrum::SampleSpectrum "
319        << "parameter p[3] <= 0. G4LEDATA dabat    320        << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
320        << Z << ". Please check and/or update i    321        << Z << ". Please check and/or update it " << G4endl;
321     }                                             322     }
322                                                   323 
323   G4double aria1 = 0.0;                           324   G4double aria1 = 0.0;
324   G4double a1 = std::max(x1,p[1]);                325   G4double a1 = std::max(x1,p[1]);
325   G4double a2 = std::min(x2,p[3]);                326   G4double a2 = std::min(x2,p[3]);
326   if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);     327   if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
327   G4double aria2 = 0.0;                           328   G4double aria2 = 0.0;
328   G4double a3 = std::max(x1,p[3]);                329   G4double a3 = std::max(x1,p[3]);
329   G4double a4 = x2;                               330   G4double a4 = x2;
330   if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);     331   if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
331                                                   332 
332   G4double aria = (aria1 + aria2)*G4UniformRan    333   G4double aria = (aria1 + aria2)*G4UniformRand();
333   G4double amaj, fun, q, x, z1, z2, dx, dx1;      334   G4double amaj, fun, q, x, z1, z2, dx, dx1;
334                                                   335 
335   //======= First aria to sample =====            336   //======= First aria to sample =====
336                                                   337 
337   if(aria <= aria1) {                             338   if(aria <= aria1) {
338                                                   339 
339     amaj = p[4];                                  340     amaj = p[4];
340     for (G4int j=5; j<iMax; j++) {                341     for (G4int j=5; j<iMax; j++) {
341       if(p[j] > amaj) amaj = p[j];                342       if(p[j] > amaj) amaj = p[j];
342     }                                             343     }
343                                                   344 
344     a1 = 1./a1;                                   345     a1 = 1./a1;
345     a2 = 1./a2;                                   346     a2 = 1./a2;
346                                                   347 
347     G4int i;                                      348     G4int i;
348     do {                                          349     do {
349                                                   350 
350       x = 1./(a2 + G4UniformRand()*(a1 - a2));    351       x = 1./(a2 + G4UniformRand()*(a1 - a2));
351       z1 = p[1];                                  352       z1 = p[1];
352       z2 = p[3];                                  353       z2 = p[3];
353       dx = (p[2] - p[1]) / 3.0;                   354       dx = (p[2] - p[1]) / 3.0;
354       dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);     355       dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
355       for (i=4; i<iMax-1; i++) {                  356       for (i=4; i<iMax-1; i++) {
356                                                   357 
357         if (i < 7) {                              358         if (i < 7) {
358           z2 = z1 + dx;                           359           z2 = z1 + dx;
359         } else if(iMax-2 == i) {                  360         } else if(iMax-2 == i) {
360           z2 = p[3];                              361           z2 = p[3];
361           break;                                  362           break;
362         } else {                                  363         } else {
363           z2 = z1*dx1;                            364           z2 = z1*dx1;
364         }                                         365         }
365         if(x >= z1 && x <= z2) break;             366         if(x >= z1 && x <= z2) break;
366         z1 = z2;                                  367         z1 = z2;
367       }                                           368       }
368       fun = p[i] + (x - z1) * (p[i+1] - p[i])/    369       fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
369                                                   370 
370       if(fun > amaj) {                            371       if(fun > amaj) {
371           G4cout << "WARNING in G4eIonisationS    372           G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
372                  << " Majoranta " << amaj         373                  << " Majoranta " << amaj
373                  << " < " << fun                  374                  << " < " << fun
374                  << " in the first aria at x=     375                  << " in the first aria at x= " << x
375                  << G4endl;                       376                  << G4endl;
376       }                                           377       }
377                                                   378 
378       q = amaj*G4UniformRand();                   379       q = amaj*G4UniformRand();
379                                                   380 
380     } while (q >= fun);                           381     } while (q >= fun);
381                                                   382 
382   //======= Second aria to sample =====           383   //======= Second aria to sample =====
383                                                   384 
384   } else {                                        385   } else {
385                                                   386 
386     amaj = std::max(p[iMax-1], Function(0.5, p    387     amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
387     a1 = 1./a3;                                   388     a1 = 1./a3;
388     a2 = 1./a4;                                   389     a2 = 1./a4;
389                                                   390 
390     do {                                          391     do {
391                                                   392 
392       x = 1./(a2 + G4UniformRand()*(a1 - a2));    393       x = 1./(a2 + G4UniformRand()*(a1 - a2));
393       fun = Function(x, p);                       394       fun = Function(x, p);
394                                                   395 
395       if(fun > amaj) {                            396       if(fun > amaj) {
396           G4cout << "WARNING in G4eIonisationS    397           G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
397                  << " Majoranta " << amaj         398                  << " Majoranta " << amaj
398                  << " < " << fun                  399                  << " < " << fun
399                  << " in the second aria at x=    400                  << " in the second aria at x= " << x
400                  << G4endl;                       401                  << G4endl;
401       }                                           402       }
402                                                   403 
403       q = amaj*G4UniformRand();                   404       q = amaj*G4UniformRand();
404                                                   405 
405     } while (q >= fun);                           406     } while (q >= fun);
406                                                   407 
407   }                                               408   }
408                                                   409 
409   p.clear();                                      410   p.clear();
410                                                   411 
411   tDelta = x*energy - bindingEnergy;              412   tDelta = x*energy - bindingEnergy;
412                                                   413 
413   if(verbose > 1) {                               414   if(verbose > 1) {
414     G4cout << "tcut(MeV)= " << tMin/MeV           415     G4cout << "tcut(MeV)= " << tMin/MeV
415            << "; tMax(MeV)= " << tMax/MeV         416            << "; tMax(MeV)= " << tMax/MeV
416            << "; x1= " << x1                      417            << "; x1= " << x1
417            << "; x2= " << x2                      418            << "; x2= " << x2
418            << "; a1= " << a1                      419            << "; a1= " << a1
419            << "; a2= " << a2                      420            << "; a2= " << a2
420            << "; x= " << x                        421            << "; x= " << x
421            << "; be= " << bindingEnergy           422            << "; be= " << bindingEnergy
422            << "; e= " << e                        423            << "; e= " << e
423            << "; tDelta= " << tDelta              424            << "; tDelta= " << tDelta
424            << G4endl;                             425            << G4endl;
425   }                                               426   }
426                                                   427 
427                                                   428 
428   return tDelta;                                  429   return tDelta;
429 }                                                 430 }
430                                                   431 
431                                                   432 
432 G4double G4eIonisationSpectrum::IntSpectrum(G4    433 G4double G4eIonisationSpectrum::IntSpectrum(G4double xMin,
433               G4double xMax,                      434               G4double xMax,
434               const G4DataVector& p) const        435               const G4DataVector& p) const
435 {                                                 436 {
436   // Please comment what IntSpectrum does         437   // Please comment what IntSpectrum does
437   G4double sum = 0.0;                             438   G4double sum = 0.0;
438   if(xMin >= xMax) return sum;                    439   if(xMin >= xMax) return sum;
439                                                   440 
440   G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2,    441   G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;
441                                                   442 
442   // Integral over interpolation aria             443   // Integral over interpolation aria
443   if(xMin < p[3]) {                               444   if(xMin < p[3]) {
444                                                   445 
445     x1 = p[1];                                    446     x1 = p[1];
446     y1 = p[4];                                    447     y1 = p[4];
447                                                   448 
448     G4double dx = (p[2] - p[1]) / 3.0;            449     G4double dx = (p[2] - p[1]) / 3.0;
449     G4double dx1= G4Exp(std::log(p[3]/p[2]) /     450     G4double dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
450                                                   451 
451     for (size_t i=0; i<19; i++) {                 452     for (size_t i=0; i<19; i++) {
452                                                   453 
453       q = 0.0;                                    454       q = 0.0;
454       if (i < 3) {                                455       if (i < 3) {
455         x2 = x1 + dx;                             456         x2 = x1 + dx;
456       } else if(18 == i) {                        457       } else if(18 == i) {
457         x2 = p[3];                                458         x2 = p[3];
458       } else {                                    459       } else {
459         x2 = x1*dx1;                              460         x2 = x1*dx1;
460       }                                           461       }
461                                                   462 
462       y2 = p[5 + i];                              463       y2 = p[5 + i];
463                                                   464 
464       if (xMax <= x1) {                           465       if (xMax <= x1) {
465         break;                                    466         break;
466       } else if (xMin < x2) {                     467       } else if (xMin < x2) {
467                                                   468 
468         xs1 = x1;                                 469         xs1 = x1;
469         xs2 = x2;                                 470         xs2 = x2;
470         ys1 = y1;                                 471         ys1 = y1;
471         ys2 = y2;                                 472         ys2 = y2;
472                                                   473 
473         if (x2 > x1) {                            474         if (x2 > x1) {
474           if (xMin > x1) {                        475           if (xMin > x1) {
475             xs1 = xMin;                           476             xs1 = xMin;
476             ys1 += (xs1 - x1)*(y2 - y1)/(x2 -     477             ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
477     }                                             478     }
478           if (xMax < x2) {                        479           if (xMax < x2) {
479             xs2 = xMax;                           480             xs2 = xMax;
480             ys2 += (xs2 - x2)*(y1 - y2)/(x1 -     481             ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
481     }                                             482     }
482           if (xs2 > xs1) {                        483           if (xs2 > xs1) {
483             q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)     484             q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
484               +  std::log(xs2/xs1)*(ys2 - ys1)    485               +  std::log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
485             sum += q;                             486             sum += q;
486             if(p.size() == 26) G4cout << "i= "    487             if(p.size() == 26) G4cout << "i= " << i << "  q= " << q << " sum= " << sum << G4endl;
487     }                                             488     }
488   }                                               489   }
489       }                                           490       }
490       x1 = x2;                                    491       x1 = x2;
491       y1 = y2;                                    492       y1 = y2;
492     }                                             493     }
493   }                                               494   }
494                                                   495 
495   // Integral over aria with parametrised form    496   // Integral over aria with parametrised formula
496                                                   497 
497   x1 = std::max(xMin, p[3]);                      498   x1 = std::max(xMin, p[3]);
498   if(x1 >= xMax) return sum;                      499   if(x1 >= xMax) return sum;
499   x2 = xMax;                                      500   x2 = xMax;
500                                                   501 
501   xs1 = 1./x1;                                    502   xs1 = 1./x1;
502   xs2 = 1./x2;                                    503   xs2 = 1./x2;
503   q = (xs1 - xs2)*(1.0 - p[0])                    504   q = (xs1 - xs2)*(1.0 - p[0])
504        - p[iMax]*std::log(x2/x1)                  505        - p[iMax]*std::log(x2/x1)
505        + (1. - p[iMax])*(x2 - x1)                 506        + (1. - p[iMax])*(x2 - x1)
506        + 1./(1. - x2) - 1./(1. - x1)              507        + 1./(1. - x2) - 1./(1. - x1)
507        + p[iMax]*std::log((1. - x2)/(1. - x1))    508        + p[iMax]*std::log((1. - x2)/(1. - x1))
508        + 0.25*p[0]*(xs1*xs1 - xs2*xs2);           509        + 0.25*p[0]*(xs1*xs1 - xs2*xs2);
509   sum += q;                                       510   sum += q;
510   if(p.size() == 26) G4cout << "param...  q= "    511   if(p.size() == 26) G4cout << "param...  q= " << q <<  " sum= " << sum << G4endl;
511                                                   512 
512   return sum;                                     513   return sum;
513 }                                                 514 }
514                                                   515 
515                                                   516 
516 G4double G4eIonisationSpectrum::AverageValue(G    517 G4double G4eIonisationSpectrum::AverageValue(G4double xMin,
517                      G4double xMax,               518                      G4double xMax,
518                const G4DataVector& p) const       519                const G4DataVector& p) const
519 {                                                 520 {
520   G4double sum = 0.0;                             521   G4double sum = 0.0;
521   if(xMin >= xMax) return sum;                    522   if(xMin >= xMax) return sum;
522                                                   523 
523   G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;    524   G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;
524                                                   525 
525   // Integral over interpolation aria             526   // Integral over interpolation aria
526   if(xMin < p[3]) {                               527   if(xMin < p[3]) {
527                                                   528 
528     x1 = p[1];                                    529     x1 = p[1];
529     y1 = p[4];                                    530     y1 = p[4];
530                                                   531 
531     G4double dx = (p[2] - p[1]) / 3.0;            532     G4double dx = (p[2] - p[1]) / 3.0;
532     G4double dx1= G4Exp(std::log(p[3]/p[2]) /     533     G4double dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
533                                                   534 
534     for (size_t i=0; i<19; i++) {                 535     for (size_t i=0; i<19; i++) {
535                                                   536 
536       if (i < 3) {                                537       if (i < 3) {
537         x2 = x1 + dx;                             538         x2 = x1 + dx;
538       } else if(18 == i) {                        539       } else if(18 == i) {
539         x2 = p[3];                                540         x2 = p[3];
540       } else {                                    541       } else {
541         x2 = x1*dx1;                              542         x2 = x1*dx1;
542       }                                           543       }
543                                                   544 
544       y2 = p[5 + i];                              545       y2 = p[5 + i];
545                                                   546 
546       if (xMax <= x1) {                           547       if (xMax <= x1) {
547         break;                                    548         break;
548       } else if (xMin < x2) {                     549       } else if (xMin < x2) {
549                                                   550 
550         xs1 = x1;                                 551         xs1 = x1;
551         xs2 = x2;                                 552         xs2 = x2;
552         ys1 = y1;                                 553         ys1 = y1;
553         ys2 = y2;                                 554         ys2 = y2;
554                                                   555 
555         if (x2 > x1) {                            556         if (x2 > x1) {
556           if (xMin > x1) {                        557           if (xMin > x1) {
557             xs1 = xMin;                           558             xs1 = xMin;
558             ys1 += (xs1 - x1)*(y2 - y1)/(x2 -     559             ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
559     }                                             560     }
560           if (xMax < x2) {                        561           if (xMax < x2) {
561             xs2 = xMax;                           562             xs2 = xMax;
562             ys2 += (xs2 - x2)*(y1 - y2)/(x1 -     563             ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
563     }                                             564     }
564           if (xs2 > xs1) {                        565           if (xs2 > xs1) {
565             sum += std::log(xs2/xs1)*(ys1*xs2     566             sum += std::log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
566                 +  ys2 - ys1;                     567                 +  ys2 - ys1;
567     }                                             568     }
568   }                                               569   }
569       }                                           570       }
570       x1 = x2;                                    571       x1 = x2;
571       y1 = y2;                                    572       y1 = y2;
572                                                   573 
573     }                                             574     }
574   }                                               575   }
575                                                   576 
576   // Integral over aria with parametrised form    577   // Integral over aria with parametrised formula
577                                                   578 
578   x1 = std::max(xMin, p[3]);                      579   x1 = std::max(xMin, p[3]);
579   if(x1 >= xMax) return sum;                      580   if(x1 >= xMax) return sum;
580   x2 = xMax;                                      581   x2 = xMax;
581                                                   582 
582   xs1 = 1./x1;                                    583   xs1 = 1./x1;
583   xs2 = 1./x2;                                    584   xs2 = 1./x2;
584                                                   585 
585   sum  += std::log(x2/x1)*(1.0 - p[0])            586   sum  += std::log(x2/x1)*(1.0 - p[0])
586         + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)      587         + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)
587         + 1./(1. - x2) - 1./(1. - x1)             588         + 1./(1. - x2) - 1./(1. - x1)
588         + (1. + p[iMax])*std::log((1. - x2)/(1    589         + (1. + p[iMax])*std::log((1. - x2)/(1. - x1))
589         + 0.5*p[0]*(xs1 - xs2);                   590         + 0.5*p[0]*(xs1 - xs2);
590                                                   591 
591   return sum;                                     592   return sum;
592 }                                                 593 }
593                                                   594 
594                                                   595 
595 void G4eIonisationSpectrum::PrintData() const     596 void G4eIonisationSpectrum::PrintData() const
596 {                                                 597 {
597   theParam->PrintData();                          598   theParam->PrintData();
598 }                                                 599 }
599                                                   600 
600 G4double G4eIonisationSpectrum::MaxEnergyOfSec    601 G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
601                    G4int, // Z = 0,               602                    G4int, // Z = 0,
602                    const G4ParticleDefinition*    603                    const G4ParticleDefinition* ) const
603 {                                                 604 {
604   return 0.5 * kineticEnergy;                     605   return 0.5 * kineticEnergy;
605 }                                                 606 }
606                                                   607