Geant4 Cross Reference |
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$ 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