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 #include "G4DNAUeharaScreenedRutherfordElastic 26 #include "G4DNAUeharaScreenedRutherfordElasticModel.hh" 27 #include "G4PhysicalConstants.hh" 27 #include "G4PhysicalConstants.hh" 28 #include "G4SystemOfUnits.hh" 28 #include "G4SystemOfUnits.hh" 29 #include "G4DNAMolecularMaterial.hh" 29 #include "G4DNAMolecularMaterial.hh" 30 #include "G4Exp.hh" 30 #include "G4Exp.hh" 31 31 32 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 33 33 34 using namespace std; 34 using namespace std; 35 35 36 // #define UEHARA_VERBOSE 36 // #define UEHARA_VERBOSE 37 37 38 //....oooOO0OOooo........oooOO0OOooo........oo 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 39 39 40 G4DNAUeharaScreenedRutherfordElasticModel:: 40 G4DNAUeharaScreenedRutherfordElasticModel:: 41 G4DNAUeharaScreenedRutherfordElasticModel(cons 41 G4DNAUeharaScreenedRutherfordElasticModel(const G4ParticleDefinition*, 42 cons << 42 const G4String& nam) : >> 43 G4VEmModel(nam), isInitialised(false) 43 { 44 { 44 // Energy limits of the models << 45 fpWaterDensity = 0; 45 intermediateEnergyLimit = 200. * eV; << 46 intermediateEnergyLimit = 200. * eV; // Switch between two final state models 46 iLowEnergyLimit = 9.*eV; << 47 47 iHighEnergyLimit = 1.*MeV; << 48 SetLowEnergyLimit(9.*eV); >> 49 // SetHighEnergyLimit(10.*keV); >> 50 SetHighEnergyLimit(1.*MeV); 48 51 49 verboseLevel = 0; 52 verboseLevel = 0; 50 // Verbosity scale: 53 // Verbosity scale: 51 // 0 = nothing 54 // 0 = nothing 52 // 1 = warning for energy non-conservation 55 // 1 = warning for energy non-conservation 53 // 2 = details of energy budget 56 // 2 = details of energy budget 54 // 3 = calculation of cross sections, file o 57 // 3 = calculation of cross sections, file openings, sampling of atoms 55 // 4 = entering in methods 58 // 4 = entering in methods >> 59 >> 60 #ifdef UEHARA_VERBOSE >> 61 if (verboseLevel) >> 62 { >> 63 G4cout << "Screened Rutherford Elastic model is constructed " << G4endl >> 64 << "Energy range: " >> 65 << LowEnergyLimit() / eV << " eV - " >> 66 << HighEnergyLimit() / MeV << " MeV" >> 67 << G4endl; >> 68 } >> 69 #endif 56 70 >> 71 fParticleChangeForGamma = 0; >> 72 >> 73 // Selection of computation method >> 74 // We do not recommend "true" usage with the current cumul. proba. settings 57 fasterCode = false; 75 fasterCode = false; 58 } 76 } 59 77 60 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 61 79 >> 80 G4DNAUeharaScreenedRutherfordElasticModel:: >> 81 ~G4DNAUeharaScreenedRutherfordElasticModel() >> 82 { >> 83 } >> 84 >> 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 86 62 void 87 void 63 G4DNAUeharaScreenedRutherfordElasticModel:: 88 G4DNAUeharaScreenedRutherfordElasticModel:: 64 Initialise(const G4ParticleDefinition* particl 89 Initialise(const G4ParticleDefinition* particle, 65 const G4DataVector& /*cuts*/) 90 const G4DataVector& /*cuts*/) 66 { 91 { 67 if(isInitialised) { return; } << 92 #ifdef UEHARA_VERBOSE 68 if (verboseLevel > 3) 93 if (verboseLevel > 3) 69 { 94 { >> 95 G4cout << "Calling G4DNAUeharaScreenedRutherfordElasticModel::Initialise()" >> 96 << G4endl; 70 } 97 } >> 98 #endif 71 99 72 if(particle->GetParticleName() != "e-") 100 if(particle->GetParticleName() != "e-") 73 { 101 { 74 G4Exception("*** WARNING: the G4DNAUeharaS 102 G4Exception("*** WARNING: the G4DNAUeharaScreenedRutherfordElasticModel is " 75 "not intented to be used with 103 "not intented to be used with another particle than the electron", 76 "",FatalException,"") ; 104 "",FatalException,"") ; 77 } 105 } 78 106 >> 107 // Energy limits >> 108 if(LowEnergyLimit() < 9.*CLHEP::eV) >> 109 { >> 110 G4Exception("*** WARNING : the G4DNAUeharaScreenedRutherfordElasticModel " >> 111 "class is not validated below 9 eV", >> 112 "",JustWarning,"") ; >> 113 } >> 114 >> 115 if (HighEnergyLimit() > 10.*CLHEP::keV) >> 116 { >> 117 G4Exception("*** WARNING: the G4DNAUeharaScreenedRutherfordElasticModel " >> 118 "class is used above 10 keV", >> 119 "",JustWarning,"") ; >> 120 } 79 121 80 if( verboseLevel>1 ) << 122 #ifdef UEHARA_VERBOSE >> 123 if( verboseLevel>0 ) 81 { 124 { 82 G4cout << "G4DNAUeharaScreenedRutherfordEl << 125 G4cout << "Screened Rutherford elastic model is initialized " << G4endl 83 << G4endl; << 126 << "Energy range: " 84 G4cout << "Energy range: " << 127 << LowEnergyLimit() / eV << " eV - " 85 << LowEnergyLimit() / eV << " eV - " << 128 << HighEnergyLimit() / MeV << " MeV" 86 << HighEnergyLimit() / MeV << " MeV" << 129 << G4endl; 87 << G4endl; << 88 } 130 } >> 131 #endif >> 132 >> 133 if (isInitialised){ return; } >> 134 89 // Constants for final state by Brenner & Za 135 // Constants for final state by Brenner & Zaider 90 // Note: the instantiation must be placed af 136 // Note: the instantiation must be placed after if (isInitialised) 91 137 92 betaCoeff= 138 betaCoeff= 93 { 139 { 94 7.51525, 140 7.51525, 95 -0.41912, 141 -0.41912, 96 7.2017E-3, 142 7.2017E-3, 97 -4.646E-5, 143 -4.646E-5, 98 1.02897E-7}; 144 1.02897E-7}; 99 145 100 deltaCoeff= 146 deltaCoeff= 101 { 147 { 102 2.9612, 148 2.9612, 103 -0.26376, 149 -0.26376, 104 4.307E-3, 150 4.307E-3, 105 -2.6895E-5, 151 -2.6895E-5, 106 5.83505E-8}; 152 5.83505E-8}; 107 153 108 gamma035_10Coeff= 154 gamma035_10Coeff= 109 { 155 { 110 -1.7013, 156 -1.7013, 111 -1.48284, 157 -1.48284, 112 0.6331, 158 0.6331, 113 -0.10911, 159 -0.10911, 114 8.358E-3, 160 8.358E-3, 115 -2.388E-4}; 161 -2.388E-4}; 116 162 117 gamma10_100Coeff = 163 gamma10_100Coeff = 118 { 164 { 119 -3.32517, 165 -3.32517, 120 0.10996, 166 0.10996, 121 -4.5255E-3, 167 -4.5255E-3, 122 5.8372E-5, 168 5.8372E-5, 123 -2.4659E-7}; 169 -2.4659E-7}; 124 170 125 gamma100_200Coeff= 171 gamma100_200Coeff= 126 { 172 { 127 2.4775E-2, 173 2.4775E-2, 128 -2.96264E-5, 174 -2.96264E-5, 129 -1.20655E-7}; 175 -1.20655E-7}; 130 176 131 // Initialize water density pointer 177 // Initialize water density pointer 132 fpWaterDensity = 178 fpWaterDensity = 133 G4DNAMolecularMaterial::Instance()-> 179 G4DNAMolecularMaterial::Instance()-> 134 GetNumMolPerVolTableFor(G4Material::GetMat 180 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER")); 135 181 136 fParticleChangeForGamma = GetParticleChangeF 182 fParticleChangeForGamma = GetParticleChangeForGamma(); 137 isInitialised = true; 183 isInitialised = true; 138 } 184 } 139 185 140 //....oooOO0OOooo........oooOO0OOooo........oo 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 141 187 142 G4double 188 G4double 143 G4DNAUeharaScreenedRutherfordElasticModel:: 189 G4DNAUeharaScreenedRutherfordElasticModel:: 144 CrossSectionPerVolume(const G4Material* materi 190 CrossSectionPerVolume(const G4Material* material, 145 const G4ParticleDefiniti 191 const G4ParticleDefinition* /*particleDefinition*/, 146 G4double ekin, 192 G4double ekin, 147 G4double, 193 G4double, 148 G4double) 194 G4double) 149 { 195 { 150 #ifdef UEHARA_VERBOSE 196 #ifdef UEHARA_VERBOSE 151 if (verboseLevel > 3) 197 if (verboseLevel > 3) 152 { 198 { 153 G4cout 199 G4cout 154 << "Calling CrossSectionPerVolume() of G4D 200 << "Calling CrossSectionPerVolume() of G4DNAUeharaScreenedRutherfordElasticModel" 155 << G4endl; 201 << G4endl; 156 } 202 } 157 #endif 203 #endif 158 204 159 // Calculate total cross section for model 205 // Calculate total cross section for model 160 206 161 G4double sigma = 0.; 207 G4double sigma = 0.; 162 if(ekin < iLowEnergyLimit || ekin > iHighEne << 163 << 164 G4double waterDensity = (*fpWaterDensity)[ma 208 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()]; 165 209 166 G4double z = 7.42; // FROM PMB 37 (1992) 184 << 210 // Use activation boundaries instead 167 G4double n = ScreeningFactor(ekin,z); << 211 //if(ekin > HighEnergyLimit() || ekin < LowEnergyLimit()) return 0.; 168 G4double crossSection = RutherfordCrossSecti << 212 169 sigma = pi * crossSection / (n * (n + 1.)); << 213 if(waterDensity!= 0.0) >> 214 { >> 215 G4double z = 7.42; // FROM PMB 37 (1992) 1841-1858 p1842 >> 216 G4double n = ScreeningFactor(ekin,z); >> 217 G4double crossSection = RutherfordCrossSection(ekin, z); >> 218 sigma = pi * crossSection / (n * (n + 1.)); 170 219 171 #ifdef UEHARA_VERBOSE 220 #ifdef UEHARA_VERBOSE 172 if (verboseLevel > 2) << 221 if (verboseLevel > 2) 173 { << 222 { 174 G4cout << "_______________________________ << 223 G4cout << "__________________________________" << G4endl; 175 G4cout << "=== G4DNAUeharaScreenedRutherfo << 224 G4cout << "=== G4DNAUeharaScreenedRutherfordElasticModel - XS INFO START" 176 << G4endl; << 225 << G4endl; 177 G4cout << "=== Kinetic energy(eV)=" << eki << 226 G4cout << "=== Kinetic energy(eV)=" << ekin/eV 178 << " particle : " << particleDefini << 227 << " particle : " << particleDefinition->GetParticleName() << G4endl; 179 G4cout << "=== Cross section per water mol << 228 G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm 180 << G4endl; << 229 << G4endl; 181 G4cout << "=== Cross section per water mol << 230 G4cout << "=== Cross section per water molecule (cm^-1)=" 182 << sigma*waterDensity/(1./cm) << G4 << 231 << sigma*waterDensity/(1./cm) << G4endl; 183 G4cout << "=== G4DNAUeharaScreenedRutherfo << 232 G4cout << "=== G4DNAUeharaScreenedRutherfordElasticModel - XS INFO END" 184 << G4endl; << 233 << G4endl; 185 } << 234 } 186 #endif 235 #endif >> 236 } 187 237 188 return sigma*waterDensity; 238 return sigma*waterDensity; 189 } 239 } 190 240 191 //....oooOO0OOooo........oooOO0OOooo........oo 241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 192 242 193 G4double 243 G4double 194 G4DNAUeharaScreenedRutherfordElasticModel::Rut 244 G4DNAUeharaScreenedRutherfordElasticModel::RutherfordCrossSection(G4double k, 195 245 G4double z) 196 { 246 { 197 // 247 // 198 // e^4 248 // e^4 / K + m_e c^2 \^2 199 // sigma_Ruth(K) = Z (Z+1) ----------------- 249 // sigma_Ruth(K) = Z (Z+1) -------------------- | --------------------- | 200 // (4 pi epsilon_0) 250 // (4 pi epsilon_0)^2 \ K * (K + 2 m_e c^2) / 201 // 251 // 202 // Where K is the electron non-relativistic 252 // Where K is the electron non-relativistic kinetic energy 203 // 253 // 204 // NIM 155, pp. 145-156, 1978 254 // NIM 155, pp. 145-156, 1978 205 255 206 G4double length = (e_squared * (k + electron 256 G4double length = (e_squared * (k + electron_mass_c2)) 207 / (4 * pi * epsilon0 * k * (k + 2 * elec 257 / (4 * pi * epsilon0 * k * (k + 2 * electron_mass_c2)); 208 G4double cross = z * (z + 1) * length * leng 258 G4double cross = z * (z + 1) * length * length; 209 259 210 return cross; 260 return cross; 211 } 261 } 212 262 213 //....oooOO0OOooo........oooOO0OOooo........oo 263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 214 264 215 G4double G4DNAUeharaScreenedRutherfordElasticM 265 G4double G4DNAUeharaScreenedRutherfordElasticModel::ScreeningFactor(G4double k, 216 266 G4double z) 217 { 267 { 218 // From Phys Med Biol 37 (1992) 1841-1858 268 // From Phys Med Biol 37 (1992) 1841-1858 219 // Z=7.42 for water 269 // Z=7.42 for water 220 270 221 const G4double constK(1.7E-5); 271 const G4double constK(1.7E-5); 222 272 223 G4double beta2; 273 G4double beta2; 224 beta2 = 1. - 1. / ((1. + k / electron_mass_c 274 beta2 = 1. - 1. / ((1. + k / electron_mass_c2) * (1. + k / electron_mass_c2)); 225 275 226 G4double etaC; 276 G4double etaC; 227 if (k < 50 * keV) 277 if (k < 50 * keV) 228 etaC = 1.198; 278 etaC = 1.198; 229 else 279 else 230 etaC = 1.13 + 3.76 * (z * z / (137 * 137 * 280 etaC = 1.13 + 3.76 * (z * z / (137 * 137 * beta2)); 231 281 232 G4double numerator = etaC * constK * std::po 282 G4double numerator = etaC * constK * std::pow(z, 2. / 3.); 233 283 234 k /= electron_mass_c2; 284 k /= electron_mass_c2; 235 285 236 G4double denominator = k * (2 + k); 286 G4double denominator = k * (2 + k); 237 287 238 G4double value = 0.; 288 G4double value = 0.; 239 if (denominator > 0.) 289 if (denominator > 0.) 240 value = numerator / denominator; // form 3 290 value = numerator / denominator; // form 3 241 291 242 return value; 292 return value; 243 } 293 } 244 294 245 //....oooOO0OOooo........oooOO0OOooo........oo 295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 246 296 247 void 297 void 248 G4DNAUeharaScreenedRutherfordElasticModel:: 298 G4DNAUeharaScreenedRutherfordElasticModel:: 249 SampleSecondaries(std::vector<G4DynamicParticl 299 SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/, 250 const G4MaterialCutsCouple* 300 const G4MaterialCutsCouple* /*couple*/, 251 const G4DynamicParticle* aDy 301 const G4DynamicParticle* aDynamicElectron, 252 G4double, 302 G4double, 253 G4double) 303 G4double) 254 { 304 { 255 #ifdef UEHARA_VERBOSE 305 #ifdef UEHARA_VERBOSE 256 if (verboseLevel > 3) 306 if (verboseLevel > 3) 257 { 307 { 258 G4cout 308 G4cout 259 << "Calling SampleSecondaries() of G4D 309 << "Calling SampleSecondaries() of G4DNAUeharaScreenedRutherfordElasticModel" 260 << G4endl; 310 << G4endl; 261 } 311 } 262 #endif 312 #endif 263 313 264 G4double electronEnergy0 = aDynamicElectron- 314 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy(); 265 if(electronEnergy0 < iLowEnergyLimit || elec << 266 return; << 267 315 268 G4double cosTheta = 0.; 316 G4double cosTheta = 0.; 269 317 270 if (electronEnergy0<intermediateEnergyLimit) 318 if (electronEnergy0<intermediateEnergyLimit) 271 { 319 { 272 #ifdef UEHARA_VERBOSE 320 #ifdef UEHARA_VERBOSE 273 if (verboseLevel > 3) 321 if (verboseLevel > 3) 274 G4cout << "---> Using Brenner & Zaider m 322 G4cout << "---> Using Brenner & Zaider model" << G4endl; 275 #endif 323 #endif 276 cosTheta = BrennerZaiderRandomizeCosTheta( 324 cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0); 277 } 325 } 278 else //if (electronEnergy0>=intermediateEner 326 else //if (electronEnergy0>=intermediateEnergyLimit) 279 { 327 { 280 #ifdef UEHARA_VERBOSE 328 #ifdef UEHARA_VERBOSE 281 if (verboseLevel > 3) 329 if (verboseLevel > 3) 282 G4cout << "---> Using Screened Rutherfor 330 G4cout << "---> Using Screened Rutherford model" << G4endl; 283 #endif 331 #endif 284 G4double z = 7.42; // FROM PMB 37 (1992) 332 G4double z = 7.42; // FROM PMB 37 (1992) 1841-1858 p1842 285 cosTheta = ScreenedRutherfordRandomizeCosT 333 cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z); 286 } 334 } 287 335 288 G4double phi = 2. * pi * G4UniformRand(); 336 G4double phi = 2. * pi * G4UniformRand(); 289 337 290 G4ThreeVector zVers = aDynamicElectron->GetM 338 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection(); 291 G4ThreeVector xVers = zVers.orthogonal(); 339 G4ThreeVector xVers = zVers.orthogonal(); 292 G4ThreeVector yVers = zVers.cross(xVers); 340 G4ThreeVector yVers = zVers.cross(xVers); 293 341 294 G4double xDir = std::sqrt(1. - cosTheta*cosT 342 G4double xDir = std::sqrt(1. - cosTheta*cosTheta); 295 G4double yDir = xDir; 343 G4double yDir = xDir; 296 xDir *= std::cos(phi); 344 xDir *= std::cos(phi); 297 yDir *= std::sin(phi); 345 yDir *= std::sin(phi); 298 346 299 G4ThreeVector zPrimeVers((xDir*xVers + yDir* 347 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers)); 300 348 301 fParticleChangeForGamma->ProposeMomentumDire 349 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()); 302 350 303 fParticleChangeForGamma->SetProposedKineticE 351 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0); 304 } 352 } 305 353 306 //....oooOO0OOooo........oooOO0OOooo........oo 354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 307 355 308 G4double 356 G4double 309 G4DNAUeharaScreenedRutherfordElasticModel:: 357 G4DNAUeharaScreenedRutherfordElasticModel:: 310 BrennerZaiderRandomizeCosTheta(G4double k) 358 BrennerZaiderRandomizeCosTheta(G4double k) 311 { 359 { 312 // d sigma_el 1 360 // d sigma_el 1 beta(K) 313 // ------------ (K) ~ ---------------------- 361 // ------------ (K) ~ --------------------------------- + --------------------------------- 314 // d Omega (1 + 2 gamma(K) - cos 362 // d Omega (1 + 2 gamma(K) - cos(theta))^2 (1 + 2 delta(K) + cos(theta))^2 315 // 363 // 316 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/( 364 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2) 317 // 365 // 318 // Phys. Med. Biol. 29 N.4 (1983) 443-447 366 // Phys. Med. Biol. 29 N.4 (1983) 443-447 319 367 320 // gamma(K), beta(K) and delta(K) are polyno 368 // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV 321 369 322 k /= eV; 370 k /= eV; 323 371 324 G4double beta = G4Exp(CalculatePolynomial(k, 372 G4double beta = G4Exp(CalculatePolynomial(k, betaCoeff)); 325 G4double delta = G4Exp(CalculatePolynomial(k 373 G4double delta = G4Exp(CalculatePolynomial(k, deltaCoeff)); 326 G4double gamma; 374 G4double gamma; 327 375 328 if (k > 100.) 376 if (k > 100.) 329 { 377 { 330 gamma = CalculatePolynomial(k, gamma100_20 378 gamma = CalculatePolynomial(k, gamma100_200Coeff); 331 // Only in this case it is not the exponen 379 // Only in this case it is not the exponent of the polynomial 332 } else 380 } else 333 { 381 { 334 if (k > 10) 382 if (k > 10) 335 { 383 { 336 gamma = G4Exp(CalculatePolynomial(k, gam 384 gamma = G4Exp(CalculatePolynomial(k, gamma10_100Coeff)); 337 } 385 } 338 else 386 else 339 { 387 { 340 gamma = G4Exp(CalculatePolynomial(k, gam 388 gamma = G4Exp(CalculatePolynomial(k, gamma035_10Coeff)); 341 } 389 } 342 } 390 } 343 391 344 // ***** Original method 392 // ***** Original method 345 393 346 if (!fasterCode) << 394 if (fasterCode == false) 347 { 395 { 348 G4double oneOverMax = 1. 396 G4double oneOverMax = 1. 349 / (1. / (4. * gamma * gamma) 397 / (1. / (4. * gamma * gamma) 350 + beta / ((2. + 2. * delta) * (2. + 2. 398 + beta / ((2. + 2. * delta) * (2. + 2. * delta))); 351 399 352 G4double cosTheta = 0.; 400 G4double cosTheta = 0.; 353 G4double leftDenominator = 0.; 401 G4double leftDenominator = 0.; 354 G4double rightDenominator = 0.; 402 G4double rightDenominator = 0.; 355 G4double fCosTheta = 0.; 403 G4double fCosTheta = 0.; 356 404 357 do 405 do 358 { 406 { 359 cosTheta = 2. * G4UniformRand()- 1.; 407 cosTheta = 2. * G4UniformRand()- 1.; 360 408 361 leftDenominator = (1. + 2.*gamma - cosTh 409 leftDenominator = (1. + 2.*gamma - cosTheta); 362 rightDenominator = (1. + 2.*delta + cosT 410 rightDenominator = (1. + 2.*delta + cosTheta); 363 if ( (leftDenominator * rightDenominator 411 if ( (leftDenominator * rightDenominator) != 0. ) 364 { 412 { 365 fCosTheta = oneOverMax * (1./(leftDeno 413 fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator) 366 + beta/(righ 414 + beta/(rightDenominator*rightDenominator)); 367 } 415 } 368 } 416 } 369 while (fCosTheta < G4UniformRand()); 417 while (fCosTheta < G4UniformRand()); 370 418 371 return cosTheta; 419 return cosTheta; 372 } 420 } 373 421 374 // ***** Alternative method using cumulative 422 // ***** Alternative method using cumulative probability 375 423 376 // if (fasterCode) << 424 else // if (fasterCode) 377 << 425 { >> 426 378 // 427 // 379 // modified by Shogo OKADA @ KEK, JP, 2016. 428 // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.) 380 // 429 // 381 // An integral of differential cross-sectio 430 // An integral of differential cross-section formula shown above this member function 382 // (integral variable: cos(theta), integral 431 // (integral variable: cos(theta), integral interval: [-1, x]) is as follows: 383 // 432 // 384 // 1.0 + x beta * ( 433 // 1.0 + x beta * (1 + x) 385 // I = --------------------- + ------------ 434 // I = --------------------- + ---------------------- (1) 386 // (a - x) * (a + 1.0) (b + x) * 435 // (a - x) * (a + 1.0) (b + x) * (b - 1.0) 387 // 436 // 388 // where a = 1.0 + 2.0 * gamma(K), b = 1.0 437 // where a = 1.0 + 2.0 * gamma(K), b = 1.0 + 2.0 * delta(K) 389 // 438 // 390 // Then, a cumulative probability (cp) is a 439 // Then, a cumulative probability (cp) is as follows: 391 // 440 // 392 // cp 1.0 + x beta * 441 // cp 1.0 + x beta * (1 + x) 393 // ---- = --------------------- + --------- 442 // ---- = --------------------- + ---------------------- (2) 394 // S (a - x) * (a + 1.0) (b + x) 443 // S (a - x) * (a + 1.0) (b + x) * (b - 1.0) 395 // 444 // 396 // where 1/S is the integral of differnetic 445 // where 1/S is the integral of differnetical cross-section (1) on interval [-1, 1] 397 // 446 // 398 // 1 2.0 2. 447 // 1 2.0 2.0 * beta 399 // --- = ----------------------- + ------- 448 // --- = ----------------------- + ----------------------- (3) 400 // S (a - 1.0) * (a + 1.0) (b + 1 449 // S (a - 1.0) * (a + 1.0) (b + 1.0) * (b - 1.0) 401 // 450 // 402 // x is calculated from the quadratic equat 451 // x is calculated from the quadratic equation derived from (2) and (3): 403 // 452 // 404 // A * x^2 + B * x + C = 0 453 // A * x^2 + B * x + C = 0 405 // 454 // 406 // where A, B, anc C are coefficients of th 455 // where A, B, anc C are coefficients of the equation: 407 // A = S * {(b - 1.0) - beta * (a + 1.0)} 456 // A = S * {(b - 1.0) - beta * (a + 1.0)} + cp * (a + 1.0) * (b - 1.0), 408 // B = S * {(b - 1.0) * (b + 1.0) + beta * 457 // B = S * {(b - 1.0) * (b + 1.0) + beta * (a - 1.0) * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * (a - b) 409 // C = S * {b * (b - 1.0) + beta * a * (a 458 // C = S * {b * (b - 1.0) + beta * a * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * ab 410 // 459 // 411 460 412 // sampling cumulative probability 461 // sampling cumulative probability 413 G4double cp = G4UniformRand(); 462 G4double cp = G4UniformRand(); 414 463 415 G4double a = 1.0 + 2.0 * gamma; 464 G4double a = 1.0 + 2.0 * gamma; 416 G4double b = 1.0 + 2.0 * delta; 465 G4double b = 1.0 + 2.0 * delta; 417 G4double a1 = a - 1.0; 466 G4double a1 = a - 1.0; 418 G4double a2 = a + 1.0; 467 G4double a2 = a + 1.0; 419 G4double b1 = b - 1.0; 468 G4double b1 = b - 1.0; 420 G4double b2 = b + 1.0; 469 G4double b2 = b + 1.0; 421 G4double c1 = a - b; 470 G4double c1 = a - b; 422 G4double c2 = a * b; 471 G4double c2 = a * b; 423 472 424 G4double S = 2.0 / (a1 * a2) + 2.0 * beta / 473 G4double S = 2.0 / (a1 * a2) + 2.0 * beta / (b1 * b2); S = 1.0 / S; 425 474 426 // coefficients for the quadratic equation 475 // coefficients for the quadratic equation 427 G4double A = S * (b1 - beta * a2) + cp * a2 476 G4double A = S * (b1 - beta * a2) + cp * a2 * b1; 428 G4double B = S * (b1 * b2 + beta * a1 * a2) 477 G4double B = S * (b1 * b2 + beta * a1 * a2) - cp * a2 * b1 * c1; 429 G4double C = S * (b * b1 + beta * a * a2) - 478 G4double C = S * (b * b1 + beta * a * a2) - cp * a2 * b1 * c2; 430 479 431 // calculate cos(theta) 480 // calculate cos(theta) 432 return (-1.0 * B + std::sqrt(B * B - 4.0 * 481 return (-1.0 * B + std::sqrt(B * B - 4.0 * A * C)) / (2.0 * A); 433 482 434 /* 483 /* 435 G4double cosTheta = -1; 484 G4double cosTheta = -1; 436 G4double cumul = 0; 485 G4double cumul = 0; 437 G4double value = 0; 486 G4double value = 0; 438 G4double leftDenominator = 0.; 487 G4double leftDenominator = 0.; 439 G4double rightDenominator = 0.; 488 G4double rightDenominator = 0.; 440 489 441 // Number of integration steps in the -1,1 490 // Number of integration steps in the -1,1 range 442 G4int iMax=200; 491 G4int iMax=200; 443 492 444 G4double random = G4UniformRand(); 493 G4double random = G4UniformRand(); 445 494 446 // Cumulate differential cross section 495 // Cumulate differential cross section 447 for (G4int i=0; i<iMax; i++) 496 for (G4int i=0; i<iMax; i++) 448 { 497 { 449 cosTheta = -1 + i*2./(iMax-1); 498 cosTheta = -1 + i*2./(iMax-1); 450 leftDenominator = (1. + 2.*gamma - cosTheta 499 leftDenominator = (1. + 2.*gamma - cosTheta); 451 rightDenominator = (1. + 2.*delta + cosThet 500 rightDenominator = (1. + 2.*delta + cosTheta); 452 if ( (leftDenominator * rightDenominator) ! 501 if ( (leftDenominator * rightDenominator) != 0. ) 453 { 502 { 454 cumul = cumul + (1./(leftDenominator*leftDe 503 cumul = cumul + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator)); 455 } 504 } 456 } 505 } 457 506 458 // Select cosTheta 507 // Select cosTheta 459 for (G4int i=0; i<iMax; i++) 508 for (G4int i=0; i<iMax; i++) 460 { 509 { 461 cosTheta = -1 + i*2./(iMax-1); 510 cosTheta = -1 + i*2./(iMax-1); 462 leftDenominator = (1. + 2.*gamma - cosTheta 511 leftDenominator = (1. + 2.*gamma - cosTheta); 463 rightDenominator = (1. + 2.*delta + cosThet 512 rightDenominator = (1. + 2.*delta + cosTheta); 464 if (cumul !=0 && (leftDenominator * rightDe 513 if (cumul !=0 && (leftDenominator * rightDenominator) != 0.) 465 value = value + (1./(leftDenominator*leftDe 514 value = value + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator)) / cumul; 466 if (random < value) break; 515 if (random < value) break; 467 } 516 } 468 517 469 return cosTheta; 518 return cosTheta; 470 */ 519 */ >> 520 } 471 521 472 << 522 return 0.; 473 //return 0.; << 474 523 475 } 524 } 476 525 477 //....oooOO0OOooo........oooOO0OOooo........oo 526 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 478 527 479 G4double 528 G4double 480 G4DNAUeharaScreenedRutherfordElasticModel:: 529 G4DNAUeharaScreenedRutherfordElasticModel:: 481 CalculatePolynomial(G4double k, 530 CalculatePolynomial(G4double k, 482 std::vector<G4double>& vec 531 std::vector<G4double>& vec) 483 { 532 { 484 // Sum_{i=0}^{size-1} vector_i k^i 533 // Sum_{i=0}^{size-1} vector_i k^i 485 // 534 // 486 // Phys. Med. Biol. 29 N.4 (1983) 443-447 535 // Phys. Med. Biol. 29 N.4 (1983) 443-447 487 536 488 G4double result = 0.; 537 G4double result = 0.; 489 size_t size = vec.size(); 538 size_t size = vec.size(); 490 539 491 while (size > 0) 540 while (size > 0) 492 { 541 { 493 size--; 542 size--; 494 543 495 result *= k; 544 result *= k; 496 result += vec[size]; 545 result += vec[size]; 497 } 546 } 498 547 499 return result; 548 return result; 500 } 549 } 501 550 502 //....oooOO0OOooo........oooOO0OOooo........oo 551 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 503 552 504 G4double 553 G4double 505 G4DNAUeharaScreenedRutherfordElasticModel:: 554 G4DNAUeharaScreenedRutherfordElasticModel:: 506 ScreenedRutherfordRandomizeCosTheta(G4double k 555 ScreenedRutherfordRandomizeCosTheta(G4double k, 507 G4double z 556 G4double z) 508 { 557 { 509 558 510 // d sigma_el sigma_Ruth(K) 559 // d sigma_el sigma_Ruth(K) 511 // ------------ (K) ~ ---------------------- 560 // ------------ (K) ~ ----------------------------- 512 // d Omega (1 + 2 n(K) - cos(the 561 // d Omega (1 + 2 n(K) - cos(theta))^2 513 // 562 // 514 // We extract cos(theta) distributed as (1 + 563 // We extract cos(theta) distributed as (1 + 2 n(K) - cos(theta))^-2 515 // 564 // 516 // Maximum is for theta=0: 1/(4 n(K)^2) (Whe 565 // Maximum is for theta=0: 1/(4 n(K)^2) (When n(K) is positive, that is always 517 // satisfied within the validity of the proc 566 // satisfied within the validity of the process) 518 // 567 // 519 // Phys. Med. Biol. 45 (2000) 3171-3194 568 // Phys. Med. Biol. 45 (2000) 3171-3194 520 569 521 // ***** Original method 570 // ***** Original method 522 571 523 if (!fasterCode) << 572 if (fasterCode == false) 524 { 573 { 525 G4double n = ScreeningFactor(k, z); 574 G4double n = ScreeningFactor(k, z); 526 575 527 G4double oneOverMax = (4. * n * n); 576 G4double oneOverMax = (4. * n * n); 528 577 529 G4double cosTheta = 0.; 578 G4double cosTheta = 0.; 530 G4double fCosTheta; 579 G4double fCosTheta; 531 580 532 do 581 do 533 { 582 { 534 cosTheta = 2. * G4UniformRand()- 1.; 583 cosTheta = 2. * G4UniformRand()- 1.; 535 fCosTheta = (1 + 2.*n - cosTheta); 584 fCosTheta = (1 + 2.*n - cosTheta); 536 if (fCosTheta !=0.) fCosTheta = oneOverM 585 if (fCosTheta !=0.) fCosTheta = oneOverMax / (fCosTheta*fCosTheta); 537 } 586 } 538 while (fCosTheta < G4UniformRand()); 587 while (fCosTheta < G4UniformRand()); 539 588 540 return cosTheta; 589 return cosTheta; 541 } 590 } 542 591 543 // ***** Alternative method using cumulative 592 // ***** Alternative method using cumulative probability 544 // if (fasterCode) << 593 else // if (fasterCode) 545 << 594 { 546 // << 595 547 // modified by Shogo OKADA @ KEK, JP, 2016.2 << 596 // 548 // << 597 // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.) 549 // The cumulative probability (cp) is calcul << 598 // 550 // the differential cross-section fomula wit << 599 // The cumulative probability (cp) is calculated by integrating 551 // << 600 // the differential cross-section fomula with cos(theta): 552 // n(K) * (1.0 + cos(theta)) << 601 // 553 // cp = --------------------------------- << 602 // n(K) * (1.0 + cos(theta)) 554 // 1.0 + 2.0 * n(K) - cos(theta) << 603 // cp = --------------------------------- 555 // << 604 // 1.0 + 2.0 * n(K) - cos(theta) 556 // Then, cos(theta) is as follows: << 605 // 557 // << 606 // Then, cos(theta) is as follows: 558 // cp * (1.0 + 2.0 * n(K)) - n << 607 // 559 // cos(theta) = ---------------------------- << 608 // cp * (1.0 + 2.0 * n(K)) - n(K) 560 // n(k) + cp << 609 // cos(theta) = -------------------------------- 561 // << 610 // n(k) + cp 562 // where, K is kinetic energy, n(K) is scree << 611 // 563 // << 612 // where, K is kinetic energy, n(K) is screeing factor, and cp is cumulative probability 564 << 613 // 565 G4double n = ScreeningFactor(k, z); << 614 566 G4double cp = G4UniformRand(); << 615 G4double n = ScreeningFactor(k, z); 567 G4double numerator = cp * (1.0 + 2.0 * n) - << 616 G4double cp = G4UniformRand(); 568 G4double denominator = n + cp; << 617 G4double numerator = cp * (1.0 + 2.0 * n) - n; 569 return numerator / denominator; << 618 G4double denominator = n + cp; 570 << 619 return numerator / denominator; 571 /* << 620 572 G4double cosTheta = -1; << 621 /* 573 G4double cumul = 0; << 622 G4double cosTheta = -1; 574 G4double value = 0; << 623 G4double cumul = 0; 575 G4double n = ScreeningFactor(k, z); << 624 G4double value = 0; 576 G4double fCosTheta; << 625 G4double n = ScreeningFactor(k, z); 577 << 626 G4double fCosTheta; 578 // Number of integration steps in the -1,1 << 627 579 G4int iMax=200; << 628 // Number of integration steps in the -1,1 range 580 << 629 G4int iMax=200; 581 G4double random = G4UniformRand(); << 630 582 << 631 G4double random = G4UniformRand(); 583 // Cumulate differential cross section << 632 584 for (G4int i=0; i<iMax; i++) << 633 // Cumulate differential cross section 585 { << 634 for (G4int i=0; i<iMax; i++) 586 cosTheta = -1 + i*2./(iMax-1); << 635 { 587 fCosTheta = (1 + 2.*n - cosTheta); << 636 cosTheta = -1 + i*2./(iMax-1); 588 if (fCosTheta !=0.) cumul = cumul + 1./(fCo << 637 fCosTheta = (1 + 2.*n - cosTheta); 589 } << 638 if (fCosTheta !=0.) cumul = cumul + 1./(fCosTheta*fCosTheta); 590 << 639 } 591 // Select cosTheta << 640 592 for (G4int i=0; i<iMax; i++) << 641 // Select cosTheta 593 { << 642 for (G4int i=0; i<iMax; i++) 594 cosTheta = -1 + i*2./(iMax-1); << 643 { 595 fCosTheta = (1 + 2.*n - cosTheta); << 644 cosTheta = -1 + i*2./(iMax-1); 596 if (cumul !=0.) value = value + (1./(fCosTh << 645 fCosTheta = (1 + 2.*n - cosTheta); 597 if (random < value) break; << 646 if (cumul !=0.) value = value + (1./(fCosTheta*fCosTheta)) / cumul; 598 } << 647 if (random < value) break; 599 return cosTheta; << 648 } 600 */ << 649 return cosTheta; 601 << 650 */ >> 651 } 602 652 603 //return 0.; << 653 return 0.; 604 } 654 } 605 655