Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 28 #include "G4DNADingfelderChargeDecreaseModel.h 29 #include "G4PhysicalConstants.hh" 30 #include "G4SystemOfUnits.hh" 31 #include "G4DNAMolecularMaterial.hh" 32 #include "G4DNAChemistryManager.hh" 33 #include "G4Log.hh" 34 #include "G4Pow.hh" 35 #include "G4Alpha.hh" 36 37 38 static G4Pow * gpow = G4Pow::GetInstance(); 39 40 //....oooOO0OOooo........oooOO0OOooo........oo 41 42 using namespace std; 43 44 //....oooOO0OOooo........oooOO0OOooo........oo 45 46 G4DNADingfelderChargeDecreaseModel::G4DNADingf 47 48 G4VEmModel(nam) 49 { 50 numberOfPartialCrossSections[0] = 0; 51 numberOfPartialCrossSections[1] = 0; 52 numberOfPartialCrossSections[2] = 0; 53 54 verboseLevel = 0; 55 // Verbosity scale: 56 // 0 = nothing 57 // 1 = warning for energy non-conservation 58 // 2 = details of energy budget 59 // 3 = calculation of cross sections, file o 60 // 4 = entering in methods 61 62 if (verboseLevel > 0) 63 { 64 G4cout << "Dingfelder charge decrease mode 65 } 66 // Selection of stationary mode 67 68 statCode = false; 69 } 70 71 //....oooOO0OOooo........oooOO0OOooo........oo 72 73 void G4DNADingfelderChargeDecreaseModel::Initi 74 75 { 76 77 if (verboseLevel > 3) 78 { 79 G4cout << "Calling G4DNADingfelderChargeDe 80 << G4endl; 81 } 82 83 // Energy limits 84 85 G4DNAGenericIonsManager *instance; 86 instance = G4DNAGenericIonsManager::Instance 87 protonDef = G4Proton::ProtonDefinition(); 88 alphaPlusPlusDef = G4Alpha::Alpha(); 89 alphaPlusDef = instance->GetIon("alpha+"); 90 hydrogenDef = instance->GetIon("hydrogen"); 91 heliumDef = instance->GetIon("helium"); 92 93 G4String proton; 94 G4String alphaPlusPlus; 95 G4String alphaPlus; 96 97 // Limits 98 99 proton = protonDef->GetParticleName(); 100 lowEnergyLimit[proton] = 100. * eV; 101 highEnergyLimit[proton] = 100. * MeV; 102 103 alphaPlusPlus = alphaPlusPlusDef->GetParticl 104 lowEnergyLimit[alphaPlusPlus] = 1. * keV; 105 highEnergyLimit[alphaPlusPlus] = 400. * MeV; 106 107 alphaPlus = alphaPlusDef->GetParticleName(); 108 lowEnergyLimit[alphaPlus] = 1. * keV; 109 highEnergyLimit[alphaPlus] = 400. * MeV; 110 111 // 112 113 if (particle==protonDef) 114 { 115 SetLowEnergyLimit(lowEnergyLimit[proton]); 116 SetHighEnergyLimit(highEnergyLimit[proton] 117 } 118 119 if (particle==alphaPlusPlusDef) 120 { 121 SetLowEnergyLimit(lowEnergyLimit[alphaPlus 122 SetHighEnergyLimit(highEnergyLimit[alphaPl 123 } 124 125 if (particle==alphaPlusDef) 126 { 127 SetLowEnergyLimit(lowEnergyLimit[alphaPlus 128 SetHighEnergyLimit(highEnergyLimit[alphaPl 129 } 130 131 // Final state 132 133 //PROTON 134 f0[0][0]=1.; 135 a0[0][0]=-0.180; 136 a1[0][0]=-3.600; 137 b0[0][0]=-18.22; 138 b1[0][0]=-1.997; 139 c0[0][0]=0.215; 140 d0[0][0]=3.550; 141 x0[0][0]=3.450; 142 x1[0][0]=5.251; 143 144 numberOfPartialCrossSections[0] = 1; 145 146 //ALPHA++ 147 f0[0][1]=1.; a0[0][1]=0.95; 148 a1[0][1]=-2.75; 149 b0[0][1]=-23.00; 150 c0[0][1]=0.215; 151 d0[0][1]=2.95; 152 x0[0][1]=3.50; 153 154 f0[1][1]=1.; 155 a0[1][1]=0.95; 156 a1[1][1]=-2.75; 157 b0[1][1]=-23.73; 158 c0[1][1]=0.250; 159 d0[1][1]=3.55; 160 x0[1][1]=3.72; 161 162 x1[0][1]=-1.; 163 b1[0][1]=-1.; 164 165 x1[1][1]=-1.; 166 b1[1][1]=-1.; 167 168 numberOfPartialCrossSections[1] = 2; 169 170 // ALPHA+ 171 f0[0][2]=1.; 172 a0[0][2]=0.65; 173 a1[0][2]=-2.75; 174 b0[0][2]=-21.81; 175 c0[0][2]=0.232; 176 d0[0][2]=2.95; 177 x0[0][2]=3.53; 178 179 x1[0][2]=-1.; 180 b1[0][2]=-1.; 181 182 numberOfPartialCrossSections[2] = 1; 183 184 // 185 186 if( verboseLevel>0 ) 187 { 188 G4cout << "Dingfelder charge decrease mode 189 << "Energy range: " 190 << LowEnergyLimit() / keV << " keV - " 191 << HighEnergyLimit() / MeV << " MeV for " 192 << particle->GetParticleName() 193 << G4endl; 194 } 195 196 // Initialize water density pointer 197 fpMolWaterDensity = G4DNAMolecularMaterial:: 198 199 if (isInitialised) 200 { return;} 201 fParticleChangeForGamma = GetParticleChangeF 202 isInitialised = true; 203 204 } 205 206 //....oooOO0OOooo........oooOO0OOooo........oo 207 208 G4double G4DNADingfelderChargeDecreaseModel::C 209 210 211 212 213 { 214 if (verboseLevel > 3) 215 { 216 G4cout 217 << "Calling CrossSectionPerVolume() of 218 << G4endl; 219 } 220 221 // Calculate total cross section for model 222 if ( 223 particleDefinition != protonDef 224 && 225 particleDefinition != alphaPlusPlusDef 226 && 227 particleDefinition != alphaPlusDef 228 ) 229 230 return 0; 231 232 G4double lowLim = 0; 233 G4double highLim = 0; 234 G4double crossSection = 0.; 235 236 G4double waterDensity = (*fpMolWaterDensity) 237 238 const G4String& particleName = particleDefin 239 240 std::map< G4String,G4double,std::less<G4Stri 241 pos1 = lowEnergyLimit.find(particleName); 242 243 if (pos1 != lowEnergyLimit.end()) 244 { 245 lowLim = pos1->second; 246 } 247 248 std::map< G4String,G4double,std::less<G4Stri 249 pos2 = highEnergyLimit.find(particleName); 250 251 if (pos2 != highEnergyLimit.end()) 252 { 253 highLim = pos2->second; 254 } 255 256 if (k >= lowLim && k <= highLim) 257 { 258 crossSection = Sum(k,particleDefinition); 259 } 260 261 if (verboseLevel > 2) 262 { 263 G4cout << "_______________________________ 264 G4cout << "G4DNADingfelderChargeDecreaeMod 265 G4cout << "Kinetic energy(eV)=" << k/eV << 266 G4cout << "Cross section per water molecul 267 G4cout << "Cross section per water molecul 268 waterDensity/(1./cm) << G4endl; 269 // material->GetAtomicNumDensityVector()[1 270 } 271 272 return crossSection*waterDensity; 273 //return crossSection*material->GetAtomicNum 274 275 } 276 277 //....oooOO0OOooo........oooOO0OOooo........oo 278 279 void G4DNADingfelderChargeDecreaseModel::Sampl 280 281 282 283 284 285 { 286 if (verboseLevel > 3) 287 { 288 G4cout 289 << "Calling SampleSecondaries() of G4D 290 << G4endl; 291 } 292 293 G4double inK = aDynamicParticle->GetKineticE 294 295 G4ParticleDefinition* definition = aDynamicP 296 297 G4double particleMass = definition->GetPDGMa 298 299 G4int finalStateIndex = RandomSelect(inK,def 300 301 G4int n = NumberOfFinalStates(definition, fi 302 G4double waterBindingEnergy = WaterBindingEn 303 G4double outgoingParticleBindingEnergy = Out 304 305 G4double outK = 0.; 306 307 if (definition==G4Proton::Proton()) 308 { 309 if (!statCode) outK = inK - n*(inK*electron 310 - waterBindingEne 311 else outK = inK; 312 } 313 314 else 315 { 316 if (!statCode) outK = inK - n*(inK*electron 317 - waterBindingEne 318 else outK = inK; 319 } 320 321 if (outK<0) 322 { 323 G4Exception("G4DNADingfelderChargeDecrease 324 FatalException,"Final kinetic energy i 325 } 326 327 fParticleChangeForGamma->ProposeTrackStatus( 328 329 if (!statCode) fParticleChangeForGamma->Prop 330 331 else 332 333 { 334 if (definition==G4Proton::Proton()) 335 fParticleChangeForGamma->ProposeLocalEner 336 + waterBindingEnergy - outgoingParticleBi 337 else 338 fParticleChangeForGamma->ProposeLocalEner 339 + waterBindingEnergy - outgoingParticleBi 340 } 341 342 auto dp = new G4DynamicParticle (OutgoingPa 343 aDynamicParticle->GetMomentumDirection() 344 outK); 345 fvect->push_back(dp); 346 347 const G4Track * theIncomingTrack = fParticle 348 G4DNAChemistryManager::Instance()->Creat 349 1, 350 theIncomingTrack); 351 352 } 353 354 //....oooOO0OOooo........oooOO0OOooo........oo 355 356 G4int G4DNADingfelderChargeDecreaseModel::Numb 357 358 359 { 360 if (particleDefinition == G4Proton::Proton() 361 return 1; 362 363 if (particleDefinition == alphaPlusPlusDef) 364 { 365 if (finalStateIndex == 0) 366 return 1; 367 return 2; 368 } 369 370 if (particleDefinition == alphaPlusDef) 371 return 1; 372 373 return 0; 374 } 375 376 //....oooOO0OOooo........oooOO0OOooo........oo 377 378 G4ParticleDefinition* G4DNADingfelderChargeDec 379 380 { 381 if (particleDefinition == G4Proton::Proton() 382 return hydrogenDef; 383 384 if (particleDefinition == alphaPlusPlusDef) 385 { 386 if (finalStateIndex == 0) 387 return alphaPlusDef; 388 return heliumDef; 389 } 390 391 if (particleDefinition == alphaPlusDef) 392 return heliumDef; 393 394 return nullptr; 395 } 396 397 //....oooOO0OOooo........oooOO0OOooo........oo 398 399 G4double G4DNADingfelderChargeDecreaseModel::W 400 401 { 402 // Ionization energy of first water shell 403 // Rad. Phys. Chem. 59 p.267 by Dingf. et al 404 // W + 10.79 eV -> W+ + e- 405 406 if (particleDefinition == G4Proton::Proton() 407 return 10.79 * eV; 408 409 if (particleDefinition == alphaPlusPlusDef) 410 { 411 // Binding energy for W+ -> W++ + e- 412 // Binding energy for W -> W+ + e- 413 // 414 // Ionization energy of first water shell 415 // Rad. Phys. Chem. 59 p.267 by Dingf. et 416 417 if (finalStateIndex == 0) 418 return 10.79 * eV; 419 420 return 10.79 * 2 * eV; 421 } 422 423 if (particleDefinition == alphaPlusDef) 424 { 425 // Binding energy for W+ -> W++ + e- 426 // Binding energy for W -> W+ + e- 427 // 428 // Ionization energy of first water shell 429 // Rad. Phys. Chem. 59 p.267 by Dingf. et 430 431 return 10.79 * eV; 432 } 433 434 return 0; 435 } 436 437 //....oooOO0OOooo........oooOO0OOooo........oo 438 439 G4double G4DNADingfelderChargeDecreaseModel::O 440 441 { 442 if (particleDefinition == G4Proton::Proton() 443 return 13.6 * eV; 444 445 if (particleDefinition == alphaPlusPlusDef) 446 { 447 // Binding energy for He+ -> He++ + e- 448 // Binding energy for He -> He+ + e- 449 450 if (finalStateIndex == 0) 451 return 54.509 * eV; 452 453 return (54.509 + 24.587) * eV; 454 } 455 456 if (particleDefinition == alphaPlusDef) 457 { 458 // Binding energy for He+ -> He++ + e- 459 // Binding energy for He -> He+ + e- 460 461 return 24.587 * eV; 462 } 463 464 return 0; 465 } 466 467 //....oooOO0OOooo........oooOO0OOooo........oo 468 469 G4double G4DNADingfelderChargeDecreaseModel::P 470 471 472 { 473 G4int particleTypeIndex = 0; 474 475 if (particleDefinition == protonDef) 476 particleTypeIndex = 0; 477 478 if (particleDefinition == alphaPlusPlusDef) 479 particleTypeIndex = 1; 480 481 if (particleDefinition == alphaPlusDef) 482 particleTypeIndex = 2; 483 484 // 485 // sigma(T) = f0 10 ^ y(log10(T/eV)) 486 // 487 // / a0 x + b0 x 488 // | 489 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x 490 // | 491 // \ a1 x + b1 x 492 // 493 // 494 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are pa 495 // 496 // f0 has been added to the code in order to 497 // 498 // From Rad. Phys. and Chem. 59 (2000) 255-2 499 // Inelastic-collision cross sections of liq 500 // 501 502 if (x1[index][particleTypeIndex] < x0[index] 503 { 504 // 505 // if x1 < x0 means that x1 and b1 will be 506 // 507 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / 508 // 509 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x 510 // 511 512 x1[index][particleTypeIndex] = x0[index][p 513 + gpow->powA((a0[index][particleTypeIn 514 / (c0[index][particleTy 515 * d0[index][particl 516 1. / (d0[index][particleTyp 517 b1[index][particleTypeIndex] = (a0[index][ 518 - a1[index][particleTypeIndex]) * x1[i 519 + b0[index][particleTypeIndex] 520 - c0[index][particleTypeIndex] 521 * gpow->powA(x1[index][particleTyp 522 - x0[index][particl 523 d0[index][particleTypeI 524 } 525 526 G4double x(G4Log(k / eV)/gpow->logZ(10)); 527 G4double y; 528 529 if (x < x0[index][particleTypeIndex]) 530 y = a0[index][particleTypeIndex] * x + b0[ 531 else if (x < x1[index][particleTypeIndex]) 532 y = a0[index][particleTypeIndex] * x + b0[ 533 - c0[index][particleTypeIndex] 534 * gpow->powA(x - x0[index][particl 535 d0[index][particleTypeI 536 else 537 y = a1[index][particleTypeIndex] * x + b1[ 538 539 return f0[index][particleTypeIndex] * gpow-> 540 541 } 542 543 G4int G4DNADingfelderChargeDecreaseModel::Rand 544 545 { 546 G4int particleTypeIndex = 0; 547 548 if (particleDefinition == protonDef) 549 particleTypeIndex = 0; 550 551 if (particleDefinition == alphaPlusPlusDef) 552 particleTypeIndex = 1; 553 554 if (particleDefinition == alphaPlusDef) 555 particleTypeIndex = 2; 556 557 const G4int n = numberOfPartialCrossSections 558 auto values(new G4double[n]); 559 G4double value(0); 560 G4int i = n; 561 562 while (i > 0) 563 { 564 i--; 565 values[i] = PartialCrossSection(k, i, part 566 value += values[i]; 567 } 568 569 value *= G4UniformRand(); 570 571 i = n; 572 while (i > 0) 573 { 574 i--; 575 576 if (values[i] > value) 577 break; 578 579 value -= values[i]; 580 } 581 582 delete[] values; 583 584 return i; 585 } 586 587 //....oooOO0OOooo........oooOO0OOooo........oo 588 589 G4double G4DNADingfelderChargeDecreaseModel::S 590 591 { 592 G4int particleTypeIndex = 0; 593 594 if (particleDefinition == protonDef) 595 particleTypeIndex = 0; 596 597 if (particleDefinition == alphaPlusPlusDef) 598 particleTypeIndex = 1; 599 600 if (particleDefinition == alphaPlusDef) 601 particleTypeIndex = 2; 602 603 G4double totalCrossSection = 0.; 604 605 for (G4int i = 0; i < numberOfPartialCrossSe 606 { 607 totalCrossSection += PartialCrossSection(k 608 } 609 610 return totalCrossSection; 611 } 612 613