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 9.3.p1)


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