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