Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 28 #include "G4DNADingfelderChargeIncreaseModel.hh" 29 #include "G4PhysicalConstants.hh" 30 #include "G4SystemOfUnits.hh" 31 #include "G4DNAMolecularMaterial.hh" 32 #include "G4Log.hh" 33 #include "G4Pow.hh" 34 #include "G4Alpha.hh" 35 36 static G4Pow * gpow = G4Pow::GetInstance(); 37 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 39 40 using namespace std; 41 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 43 44 G4DNADingfelderChargeIncreaseModel::G4DNADingfelderChargeIncreaseModel(const G4ParticleDefinition*, 45 const G4String& nam) : 46 G4VEmModel(nam) 47 { 48 if (verboseLevel > 0) 49 { 50 G4cout << "Dingfelder charge increase model is constructed " << G4endl; 51 } 52 } 53 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 55 56 void G4DNADingfelderChargeIncreaseModel::Initialise(const G4ParticleDefinition* particle, 57 const G4DataVector& /*cuts*/) 58 { 59 60 if (verboseLevel > 3) 61 { 62 G4cout << "Calling G4DNADingfelderChargeIncreaseModel::Initialise()" 63 << G4endl; 64 } 65 66 // Energy limits 67 68 G4DNAGenericIonsManager *instance; 69 instance = G4DNAGenericIonsManager::Instance(); 70 hydrogenDef = instance->GetIon("hydrogen"); 71 alphaPlusPlusDef = G4Alpha::Alpha(); 72 alphaPlusDef = instance->GetIon("alpha+"); 73 heliumDef = instance->GetIon("helium"); 74 75 G4String hydrogen; 76 G4String alphaPlus; 77 G4String helium; 78 79 // Limits 80 81 hydrogen = hydrogenDef->GetParticleName(); 82 lowEnergyLimit[hydrogen] = 100. * eV; 83 highEnergyLimit[hydrogen] = 100. * MeV; 84 85 alphaPlus = alphaPlusDef->GetParticleName(); 86 lowEnergyLimit[alphaPlus] = 1. * keV; 87 highEnergyLimit[alphaPlus] = 400. * MeV; 88 89 helium = heliumDef->GetParticleName(); 90 lowEnergyLimit[helium] = 1. * keV; 91 highEnergyLimit[helium] = 400. * MeV; 92 93 // 94 95 if (particle==hydrogenDef) 96 { 97 SetLowEnergyLimit(lowEnergyLimit[hydrogen]); 98 SetHighEnergyLimit(highEnergyLimit[hydrogen]); 99 } 100 101 if (particle==alphaPlusDef) 102 { 103 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]); 104 SetHighEnergyLimit(highEnergyLimit[alphaPlus]); 105 } 106 107 if (particle==heliumDef) 108 { 109 SetLowEnergyLimit(lowEnergyLimit[helium]); 110 SetHighEnergyLimit(highEnergyLimit[helium]); 111 } 112 113 // Final state 114 115 //ALPHA+ 116 117 f0[0][0]=1.; 118 a0[0][0]=2.25; 119 a1[0][0]=-0.75; 120 b0[0][0]=-32.10; 121 c0[0][0]=0.600; 122 d0[0][0]=2.40; 123 x0[0][0]=4.60; 124 125 x1[0][0]=-1.; 126 b1[0][0]=-1.; 127 128 numberOfPartialCrossSections[0]=1; 129 130 //HELIUM 131 132 f0[0][1]=1.; 133 a0[0][1]=2.25; 134 a1[0][1]=-0.75; 135 b0[0][1]=-30.93; 136 c0[0][1]=0.590; 137 d0[0][1]=2.35; 138 x0[0][1]=4.29; 139 140 f0[1][1]=1.; 141 a0[1][1]=2.25; 142 a1[1][1]=-0.75; 143 b0[1][1]=-32.61; 144 c0[1][1]=0.435; 145 d0[1][1]=2.70; 146 x0[1][1]=4.45; 147 148 x1[0][1]=-1.; 149 b1[0][1]=-1.; 150 151 x1[1][1]=-1.; 152 b1[1][1]=-1.; 153 154 numberOfPartialCrossSections[1]=2; 155 156 // 157 158 if( verboseLevel>0 ) 159 { 160 G4cout << "Dingfelder charge increase model is initialized " << G4endl 161 << "Energy range: " 162 << LowEnergyLimit() / keV << " keV - " 163 << HighEnergyLimit() / MeV << " MeV for " 164 << particle->GetParticleName() 165 << G4endl; 166 } 167 168 // Initialize water density pointer 169 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER")); 170 171 if (isInitialised) return; 172 173 fParticleChangeForGamma = GetParticleChangeForGamma(); 174 isInitialised = true; 175 } 176 177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 178 179 G4double G4DNADingfelderChargeIncreaseModel::CrossSectionPerVolume(const G4Material* material, 180 const G4ParticleDefinition* particleDefinition, 181 G4double k, 182 G4double, 183 G4double) 184 { 185 if (verboseLevel > 3) 186 { 187 G4cout 188 << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeIncreaseModel" 189 << G4endl; 190 } 191 192 // Calculate total cross section for model 193 194 if ( 195 particleDefinition != hydrogenDef 196 && 197 particleDefinition != alphaPlusDef 198 && 199 particleDefinition != heliumDef 200 ) 201 202 return 0; 203 204 G4double lowLim = 0; 205 G4double highLim = 0; 206 G4double totalCrossSection = 0.; 207 208 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()]; 209 210 const G4String& particleName = particleDefinition->GetParticleName(); 211 212 std::map< G4String,G4double,std::less<G4String> >::iterator pos1; 213 pos1 = lowEnergyLimit.find(particleName); 214 215 if (pos1 != lowEnergyLimit.end()) 216 { 217 lowLim = pos1->second; 218 } 219 220 std::map< G4String,G4double,std::less<G4String> >::iterator pos2; 221 pos2 = highEnergyLimit.find(particleName); 222 223 if (pos2 != highEnergyLimit.end()) 224 { 225 highLim = pos2->second; 226 } 227 228 if (k >= lowLim && k <= highLim) 229 { 230 //HYDROGEN 231 if (particleDefinition == hydrogenDef) 232 { 233 const G4double aa = 2.835; 234 const G4double bb = 0.310; 235 const G4double cc = 2.100; 236 const G4double dd = 0.760; 237 const G4double fac = 1.0e-18; 238 const G4double rr = 13.606 * eV; 239 240 G4double t = k / (proton_mass_c2/electron_mass_c2); 241 G4double x = t / rr; 242 G4double temp = 4.0 * pi * Bohr_radius/nm * Bohr_radius/nm * fac; 243 G4double sigmal = temp * cc * (gpow->powA(x,dd)); 244 G4double sigmah = temp * (aa * G4Log(1.0 + x) + bb) / x; 245 totalCrossSection = 1.0/(1.0/sigmal + 1.0/sigmah) *m*m; 246 } 247 else 248 { 249 totalCrossSection = Sum(k,particleDefinition); 250 } 251 } 252 253 if (verboseLevel > 2) 254 { 255 G4cout << "__________________________________" << G4endl; 256 G4cout << "G4DNADingfelderChargeIncreaseModel - XS INFO START" << G4endl; 257 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleName << G4endl; 258 G4cout << "Cross section per water molecule (cm^2)=" << totalCrossSection/cm/cm << G4endl; 259 G4cout << "Cross section per water molecule (cm^-1)=" << totalCrossSection*waterDensity/(1./cm) << G4endl; 260 // G4cout << " - Cross section per water molecule (cm^-1)=" 261 // << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl; 262 G4cout << "G4DNADingfelderChargeIncreaseModel - XS INFO END" << G4endl; 263 } 264 265 return totalCrossSection*waterDensity; 266 // return totalCrossSection*material->GetAtomicNumDensityVector()[1]; 267 268 } 269 270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 271 272 void G4DNADingfelderChargeIncreaseModel::SampleSecondaries(std::vector< 273 G4DynamicParticle*>* fvect, 274 const G4MaterialCutsCouple* /*couple*/, 275 const G4DynamicParticle* aDynamicParticle, 276 G4double, 277 G4double) 278 { 279 if (verboseLevel > 3) 280 { 281 G4cout 282 << "Calling SampleSecondaries() of G4DNADingfelderChargeIncreaseModel" 283 << G4endl; 284 } 285 286 if (!statCode) fParticleChangeForGamma->ProposeLocalEnergyDeposit(0.); 287 288 G4ParticleDefinition* definition = aDynamicParticle->GetDefinition(); 289 290 G4double particleMass = definition->GetPDGMass(); 291 292 G4double inK = aDynamicParticle->GetKineticEnergy(); 293 294 G4int finalStateIndex = RandomSelect(inK,definition); 295 296 G4int n = NumberOfFinalStates(definition,finalStateIndex); 297 298 G4double outK = 0.; 299 300 if (!statCode) outK = inK - IncomingParticleBindingEnergyConstant(definition,finalStateIndex); 301 302 else outK = inK; 303 304 if (statCode) 305 fParticleChangeForGamma-> 306 ProposeLocalEnergyDeposit(IncomingParticleBindingEnergyConstant(definition,finalStateIndex)); 307 308 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill); 309 310 G4double electronK; 311 if (definition == hydrogenDef) electronK = inK*electron_mass_c2/proton_mass_c2; 312 else electronK = inK*electron_mass_c2/(particleMass); 313 314 if (outK<0) 315 { 316 G4Exception("G4DNADingfelderChargeIncreaseModel::SampleSecondaries","em0004", 317 FatalException,"Final kinetic energy is negative."); 318 } 319 320 auto dp = new G4DynamicParticle(OutgoingParticleDefinition(definition,finalStateIndex), 321 aDynamicParticle->GetMomentumDirection(), 322 outK); 323 324 fvect->push_back(dp); 325 326 n = n - 1; 327 328 while (n>0) 329 { 330 n--; 331 fvect->push_back(new G4DynamicParticle 332 (G4Electron::Electron(), aDynamicParticle->GetMomentumDirection(), electronK) ); 333 } 334 } 335 336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 337 338 G4int G4DNADingfelderChargeIncreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition, 339 G4int finalStateIndex) 340 341 { 342 343 if (particleDefinition == hydrogenDef) 344 return 2; 345 346 if (particleDefinition == alphaPlusDef) 347 return 2; 348 349 if (particleDefinition == heliumDef) 350 { 351 if (finalStateIndex == 0) 352 return 2; 353 return 3; 354 } 355 356 return 0; 357 } 358 359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 360 361 G4ParticleDefinition* G4DNADingfelderChargeIncreaseModel::OutgoingParticleDefinition(G4ParticleDefinition* particleDefinition, 362 G4int finalStateIndex) 363 { 364 365 if (particleDefinition == hydrogenDef) 366 return G4Proton::Proton(); 367 368 if (particleDefinition == alphaPlusDef) 369 return alphaPlusPlusDef; 370 371 if (particleDefinition == heliumDef) 372 { 373 if (finalStateIndex == 0) 374 return alphaPlusDef; 375 return alphaPlusPlusDef; 376 } 377 378 return nullptr; 379 } 380 381 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 382 383 G4double G4DNADingfelderChargeIncreaseModel::IncomingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition, 384 G4int finalStateIndex) 385 { 386 387 if (particleDefinition == hydrogenDef) 388 return 13.6 * eV; 389 390 if (particleDefinition == alphaPlusDef) 391 { 392 // Binding energy for He+ -> He++ + e- 54.509 eV 393 // Binding energy for He -> He+ + e- 24.587 eV 394 return 54.509 * eV; 395 } 396 397 if (particleDefinition == heliumDef) 398 { 399 // Binding energy for He+ -> He++ + e- 54.509 eV 400 // Binding energy for He -> He+ + e- 24.587 eV 401 402 if (finalStateIndex == 0) 403 return 24.587 * eV; 404 return (54.509 + 24.587) * eV; 405 } 406 407 return 0; 408 } 409 410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 411 412 G4double G4DNADingfelderChargeIncreaseModel::PartialCrossSection(const G4double& k, 413 const G4int& index, 414 const G4ParticleDefinition* particleDefinition) 415 { 416 G4int particleTypeIndex = 0; 417 418 if (particleDefinition == alphaPlusDef) 419 particleTypeIndex = 0; 420 421 if (particleDefinition == heliumDef) 422 particleTypeIndex = 1; 423 424 // 425 // sigma(T) = f0 10 ^ y(log10(T/eV)) 426 // 427 // / a0 x + b0 x < x0 428 // | 429 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1 430 // | 431 // \ a1 x + b1 x >= x1 432 // 433 // 434 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++) 435 // 436 // f0 has been added to the code in order to manage partial (shell-dependent) cross sections 437 // (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1) 438 // 439 // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al. 440 // Inelastic-collision cross sections of liquid water for interactions of energetic proton 441 // 442 443 if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex]) 444 { 445 // 446 // if x1 < x0 means that x1 and b1 will be calculated with the following formula 447 // (this piece of code is run on all alphas and not on protons) 448 // 449 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1)) 450 // 451 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0 452 // 453 454 x1[index][particleTypeIndex] = x0[index][particleTypeIndex] 455 + gpow->powA((a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) 456 / (c0[index][particleTypeIndex] 457 * d0[index][particleTypeIndex]), 458 1. / (d0[index][particleTypeIndex] - 1.)); 459 b1[index][particleTypeIndex] = (a0[index][particleTypeIndex] 460 - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex] 461 + b0[index][particleTypeIndex] 462 - c0[index][particleTypeIndex] 463 * gpow->powA(x1[index][particleTypeIndex] 464 - x0[index][particleTypeIndex], 465 d0[index][particleTypeIndex]); 466 } 467 468 G4double x(G4Log(k / eV)/gpow->logZ(10)); 469 G4double y; 470 471 if (x < x0[index][particleTypeIndex]) 472 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]; 473 else if (x < x1[index][particleTypeIndex]) 474 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] 475 - c0[index][particleTypeIndex] 476 * gpow->powA(x - x0[index][particleTypeIndex], 477 d0[index][particleTypeIndex]); 478 else 479 y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex]; 480 481 return f0[index][particleTypeIndex] * gpow->powA(10., y) * m * m; 482 483 } 484 485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 486 487 G4int G4DNADingfelderChargeIncreaseModel::RandomSelect(const G4double& k, 488 const G4ParticleDefinition* particleDefinition) 489 { 490 G4int particleTypeIndex = 0; 491 492 if (particleDefinition == hydrogenDef) 493 return 0; 494 495 if (particleDefinition == alphaPlusDef) 496 particleTypeIndex = 0; 497 498 if (particleDefinition == heliumDef) 499 particleTypeIndex = 1; 500 501 const G4int n = numberOfPartialCrossSections[particleTypeIndex]; 502 auto values(new G4double[n]); 503 G4double value = 0; 504 G4int i = n; 505 506 while (i > 0) 507 { 508 i--; 509 values[i] = PartialCrossSection(k, i, particleDefinition); 510 value += values[i]; 511 } 512 513 value *= G4UniformRand(); 514 515 i = n; 516 while (i > 0) 517 { 518 i--; 519 520 if (values[i] > value) 521 break; 522 523 value -= values[i]; 524 } 525 526 delete[] values; 527 528 return i; 529 } 530 531 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 532 533 G4double G4DNADingfelderChargeIncreaseModel::Sum(const G4double& k, 534 const G4ParticleDefinition* particleDefinition) 535 { 536 G4int particleTypeIndex = 0; 537 538 if (particleDefinition == alphaPlusDef) 539 particleTypeIndex = 0; 540 541 if (particleDefinition == heliumDef) 542 particleTypeIndex = 1; 543 544 G4double totalCrossSection = 0.; 545 546 for (G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++) 547 { 548 totalCrossSection += PartialCrossSection(k, i, particleDefinition); 549 } 550 551 return totalCrossSection; 552 } 553