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 // Reference: 27 // A.D. Dominguez-Munoz, M.I. Gallardo, M.C 28 // Z. Francis, S. Incerti, M.A. Cortes-Gira 29 // Radiat. Phys. Chem. 199 (2022) 110363. 30 // 31 // Class authors: 32 // A.D. Dominguez-Munoz 33 // M.A. Cortes-Giraldo (miancortes -at- us. 34 // 35 // Class creation: 2022-03-03 36 // 37 // 38 39 #include "G4DNARPWBAIonisationModel.hh" 40 #include "G4PhysicalConstants.hh" 41 #include "G4SystemOfUnits.hh" 42 #include "G4UAtomicDeexcitation.hh" 43 #include "G4LossTableManager.hh" 44 #include "G4DNAChemistryManager.hh" 45 #include "G4DNAMolecularMaterial.hh" 46 #include "G4DNABornAngle.hh" 47 #include "G4Exp.hh" 48 using namespace std; 49 //....oooOO0OOooo........oooOO0OOooo........oo 50 G4DNARPWBAIonisationModel::G4DNARPWBAIonisatio 51 const G4ParticleDefinition*, const G4String& 52 : G4VEmModel(nam) 53 { 54 // Verbosity scale: 55 // 0 = nothing 56 // 1 = warning for energy non-conservation 57 // 2 = details of energy budget 58 // 3 = calculation of cross sections, file o 59 // 4 = entering in methods 60 61 if(verboseLevel > 0) 62 { 63 G4cout << "RPWBA ionisation model is const 64 } 65 SetDeexcitationFlag(true); 66 SetAngularDistribution(new G4DNABornAngle()) 67 } 68 69 //....oooOO0OOooo........oooOO0OOooo........oo 70 71 G4DNARPWBAIonisationModel::~G4DNARPWBAIonisati 72 { 73 eVecm.clear(); 74 pVecm.clear(); 75 } 76 //....oooOO0OOooo........oooOO0OOooo........oo 77 78 G4bool G4DNARPWBAIonisationModel::InEnergyLimi 79 { 80 if(lowEnergyLimit == highEnergyLimit) 81 { 82 G4Exception("G4DNARPWBAIonisationModel::In 83 FatalException, "lowEnergyLimi 84 } 85 return k >= lowEnergyLimit && k <= highEnerg 86 } 87 88 //....oooOO0OOooo........oooOO0OOooo........oo 89 void G4DNARPWBAIonisationModel::InitialiseForP 90 const G4ParticleDefinition* part) 91 { 92 if(part != fProtonDef) 93 { 94 G4Exception("G4DNARPWBAIonisationModel::In 95 FatalException, "Model not app 96 } 97 // Energy limits 98 G4String fileProton("dna/sigma_ionisation_p_ 99 G4double scaleFactor = 1 * cm * cm; 100 const char *path = G4FindDataDir("G4LEDAT 101 lowEnergyLimit = 100. * MeV; 102 highEnergyLimit = 300. * MeV; 103 104 if(LowEnergyLimit() < lowEnergyLimit || High 105 { 106 G4ExceptionDescription ed; 107 ed << "Model is applicable from "<<lowEner 108 G4Exception("G4DNARPWBAIonisationModel::In 109 FatalException, ed); 110 } 111 112 fpTotalCrossSection = make_unique<G4DNACross 113 new G4LogLogInterpolation, eV, scaleFactor 114 fpTotalCrossSection->LoadData(fileProton); 115 116 // Final state 117 118 std::ostringstream pFullFileName; 119 fasterCode ? pFullFileName 120 << path << "/dna/sigmadiff_cu 121 : pFullFileName << path << "/dna/ 122 std::ifstream pDiffCrossSection(pFullFileNam 123 if(!pDiffCrossSection) 124 { 125 G4ExceptionDescription exceptionDescriptio 126 exceptionDescription << "Missing data file 127 G4Exception("G4DNARPWBAIonisationModel::In 128 FatalException, exceptionDescr 129 } 130 131 pTdummyVec.push_back(0.); 132 while(!pDiffCrossSection.eof()) 133 { 134 G4double tDummy; 135 G4double eDummy; 136 pDiffCrossSection >> tDummy >> eDummy; 137 if(tDummy != pTdummyVec.back()) 138 { 139 pTdummyVec.push_back(tDummy); 140 } 141 142 for(G4int j = 0; j < 5; j++) 143 { 144 pDiffCrossSection >> pDiffCrossSectionDa 145 146 if(fasterCode) 147 { 148 pNrjTransfData[j][tDummy][pDiffCrossSe 149 eDummy; 150 pProbaShellMap[j][tDummy].push_back( 151 pDiffCrossSectionData[j][tDummy][eDu 152 } 153 154 // SI - only if eof is not reached ! 155 if(!pDiffCrossSection.eof() && !fasterCo 156 { 157 pDiffCrossSectionData[j][tDummy][eDumm 158 } 159 160 if(!fasterCode) 161 { 162 pVecm[tDummy].push_back(eDummy); 163 } 164 } 165 } 166 167 // be careful about this 168 // SetLowEnergyLimit(lowEnergyLimit); 169 // SetHighEnergyLimit(highEnergyLimit); 170 } 171 172 void G4DNARPWBAIonisationModel::Initialise(con 173 con 174 { 175 if(isInitialised) 176 { 177 return; 178 } 179 if(verboseLevel > 3) 180 { 181 G4cout << "Calling G4DNARPWBAIonisationMod 182 << particle->GetParticleName() << G 183 } 184 185 InitialiseForProton(particle); 186 187 if(verboseLevel > 0) 188 { 189 G4cout << "RPWBA ionisation model is initi 190 << "Energy range: " << LowEnergyLim 191 << HighEnergyLimit() / MeV << " MeV 192 << particle->GetParticleName() << G 193 } 194 195 // Initialize water density pointer 196 if(G4Material::GetMaterial("G4_WATER") != nu 197 { 198 fpMolWaterDensity = 199 G4DNAMolecularMaterial::Instance()->GetN 200 G4Material::GetMaterial("G4_WATER")); 201 } 202 else 203 { 204 G4ExceptionDescription exceptionDescriptio 205 exceptionDescription << "G4_WATER does not 206 G4Exception("G4DNARPWBAIonisationModel::In 207 FatalException, exceptionDescr 208 } 209 fAtomDeexcitation = G4LossTableManager 210 fParticleChangeForGamma = GetParticleChangeF 211 isInitialised = true; 212 } 213 214 //....oooOO0OOooo........oooOO0OOooo........oo 215 216 G4double G4DNARPWBAIonisationModel::CrossSecti 217 const G4Material* material, const G4Particle 218 G4double ekin, G4double, G4double) 219 { 220 if(particleDefinition != fProtonDef) 221 { 222 G4Exception("G4DNARPWBAIonisationModel::Cr 223 FatalException, "Model not app 224 } 225 if(verboseLevel > 3) 226 { 227 G4cout << "Calling CrossSectionPerVolume() 228 << G4endl; 229 } 230 G4double sigma; 231 G4double waterDensity = (*fpMolWaterDensity) 232 233 if(InEnergyLimit(ekin)) 234 { 235 sigma = fpTotalCrossSection->FindValue(eki 236 } 237 else 238 { 239 // nput energy is outside this interval th 240 // should add a warning or exception ? 241 return 0; 242 } 243 244 if(verboseLevel > 2) 245 { 246 G4cout << "_______________________________ 247 G4cout << "G4DNARPWBAIonisationModel - XS 248 G4cout << "Kinetic energy(eV)=" << ekin / 249 << " particle : " << fProtonDef->Ge 250 G4cout << "Cross section per water molecul 251 << G4endl; 252 G4cout << "Cross section per water molecul 253 << sigma * waterDensity / (1. / cm) 254 G4cout << "G4DNARPWBAIonisationModel - XS 255 } 256 257 return sigma * waterDensity; 258 } 259 260 //....oooOO0OOooo........oooOO0OOooo........oo 261 262 void G4DNARPWBAIonisationModel::SampleSecondar 263 std::vector<G4DynamicParticle*>* fvect, cons 264 const G4DynamicParticle* particle, G4double, 265 { 266 if(verboseLevel > 3) 267 { 268 G4cout << "Calling SampleSecondaries() of 269 << G4endl; 270 } 271 G4double k = particle->GetKineticEnergy(); 272 if(InEnergyLimit(k)) 273 { 274 G4ParticleMomentum primaryDirection = part 275 G4double particleMass = particle->GetDefi 276 G4double totalEnergy = k + particleMass; 277 G4double pSquare = k * (totalEnergy 278 G4double totalMomentum = std::sqrt(pSquare 279 G4int ionizationShell; 280 281 if(!fasterCode) 282 { 283 ionizationShell = RandomSelect(k); 284 } 285 else 286 { 287 // fasterCode = true 288 do 289 { 290 ionizationShell = RandomSelect(k); 291 } while(k < 19 * eV && ionizationShell = 292 particle->GetDefinition() == G4E 293 } 294 295 G4double bindingEnergy = waterStructure.Io 296 297 // SI: additional protection if tcs interp 298 if(k < bindingEnergy) 299 { 300 return; 301 } 302 // 303 G4double secondaryKinetic; 304 if(!fasterCode) 305 { 306 secondaryKinetic = RandomizeEjectedElect 307 } 308 else 309 { 310 secondaryKinetic = 311 RandomizeEjectedElectronEnergyFromCumu 312 } 313 314 G4int Z = 8; // water Z (6 Oxygen + 2 hyd 315 G4ThreeVector deltaDirection = 316 GetAngularDistribution()->SampleDirectio 317 particle, secondaryKinetic, Z, ionizat 318 319 if(secondaryKinetic > 0){ 320 auto dp = new G4DynamicParticle(G4Electr 321 secondar 322 fvect->push_back(dp); 323 } 324 325 if(particle->GetDefinition() == G4Electron 326 G4double deltaTotalMomentum = std::sqrt( 327 secondaryKinetic * (secondaryKinetic + 328 329 G4double finalPx = totalMomentum * prima 330 deltaTotalMomentum * 331 G4double finalPy = totalMomentum * prima 332 deltaTotalMomentum * 333 G4double finalPz = totalMomentum * prima 334 deltaTotalMomentum * 335 G4double finalMomentum = 336 std::sqrt(finalPx * finalPx + finalPy 337 finalPx /= finalMomentum; 338 finalPy /= finalMomentum; 339 finalPz /= finalMomentum; 340 G4ThreeVector direction; 341 direction.set(finalPx, finalPy, finalPz) 342 fParticleChangeForGamma->ProposeMomentum 343 } 344 else 345 { 346 fParticleChangeForGamma->ProposeMomentum 347 } 348 349 // AM: sample deexcitation 350 // here we assume that H_{2}O electronic l 351 // this can be considered true with a roug 352 353 size_t secNumberInit; // need to know at 354 // secondaries 355 size_t 356 secNumberFinal; // So I'll make the dif 357 358 G4double scatteredEnergy = k - bindingEner 359 360 // SI: only atomic deexcitation from K she 361 if((fAtomDeexcitation != nullptr) && ioniz 362 { 363 const G4AtomicShell* shell = 364 fAtomDeexcitation->GetAtomicShell(Z, G 365 secNumberInit = fvect->size(); 366 fAtomDeexcitation->GenerateParticles(fve 367 secNumberFinal = fvect->size(); 368 369 if(secNumberFinal > secNumberInit){ 370 for(size_t i = secNumberInit; i < secN 371 if(bindingEnergy >= ((*fvect)[i])->G 372 { 373 bindingEnergy -= ((*fvect)[i])->Ge 374 }else{ 375 delete(*fvect)[i]; 376 (*fvect)[i] = nullptr; 377 } 378 } 379 } 380 } 381 382 // This should never happen 383 if(bindingEnergy < 0.0) 384 { 385 G4Exception("G4DNARPWBAIonisatioModel::S 386 FatalException, "Negative lo 387 } 388 389 // bindingEnergy has been decreased 390 // by the amount of energy taken away by d 391 if(!statCode){ 392 fParticleChangeForGamma->SetProposedKine 393 fParticleChangeForGamma->ProposeLocalEne 394 }else{ 395 fParticleChangeForGamma->SetProposedKine 396 fParticleChangeForGamma->ProposeLocalEne 397 } 398 const G4Track* theIncomingTrack = 399 fParticleChangeForGamma->GetCurrentTrack 400 G4DNAChemistryManager::Instance()->CreateW 401 eIonizedMolecule, ionizationShell, theIn 402 } 403 } 404 405 //....oooOO0OOooo........oooOO0OOooo........oo 406 407 G4double G4DNARPWBAIonisationModel::RandomizeE 408 const G4double& k, const G4int& shell) 409 { 410 G4double maximumKineticEnergyTransfer = 411 4. * (electron_mass_c2 / proton_mass_c2) * 412 G4double IonisationEnergyInShell = waterStru 413 G4double kIneV = k / eV; 414 415 G4double crossSectionMaximum = 0.; 416 for(G4double value = IonisationEnergyInShell 417 value <= 4. * IonisationEnergyInShell; v 418 { 419 G4double differentialCrossSection = 420 DifferentialCrossSection(kIneV, value / 421 if(differentialCrossSection >= crossSectio 422 { 423 crossSectionMaximum = differentialCrossS 424 } 425 } 426 427 G4double secondaryElectronKineticEnergy; 428 do 429 { 430 secondaryElectronKineticEnergy = 431 G4UniformRand() * maximumKineticEnergyTr 432 } while(G4UniformRand() * crossSectionMaximu 433 DifferentialCrossSection(kIneV, 434 (secondaryElectronKineticEnergy + 435 IonisationEnergyInShell) / eV, 436 437 return secondaryElectronKineticEnergy; 438 } 439 440 //....oooOO0OOooo........oooOO0OOooo........oo 441 G4double G4DNARPWBAIonisationModel::Differenti 442 const G4double& kine, const G4double& energy 443 const G4int& ionizationLevelIndex) 444 { 445 G4double k = kine; 446 G4double sigma = 0.; 447 if(energyTransfer >= 448 waterStructure.IonisationEnergy(ionizatio 449 { 450 G4double valueT1 = 0; 451 G4double valueT2 = 0; 452 G4double valueE21 = 0; 453 G4double valueE22 = 0; 454 G4double valueE12 = 0; 455 G4double valueE11 = 0; 456 457 G4double xs11 = 0; 458 G4double xs12 = 0; 459 G4double xs21 = 0; 460 G4double xs22 = 0; 461 462 // Protection against out of boundary acce 463 464 if(k == pTdummyVec.back()) 465 { 466 k = k * (1. - 1e-12); 467 } 468 // k should be in eV and energy transfer e 469 auto t2 = std::upper_bound(pTdummyVec.begi 470 auto t1 = t2 - 1; 471 472 auto e12 = std::upper_bound(pVecm[(*t1)].b 473 energyTransfer 474 auto e11 = e12 - 1; 475 476 auto e22 = std::upper_bound(pVecm[(*t2)].b 477 energyTransfer 478 auto e21 = e22 - 1; 479 480 valueT1 = *t1; 481 valueT2 = *t2; 482 valueE21 = *e21; 483 valueE22 = *e22; 484 valueE12 = *e12; 485 valueE11 = *e11; 486 487 xs11 = pDiffCrossSectionData[ionizationLev 488 xs12 = pDiffCrossSectionData[ionizationLev 489 xs21 = pDiffCrossSectionData[ionizationLev 490 xs22 = pDiffCrossSectionData[ionizationLev 491 492 G4double xsProduct = xs11 * xs12 * xs21 * 493 if(xsProduct != 0.) 494 { 495 sigma = 496 QuadInterpolator(valueE11, valueE12, v 497 xs21, xs22, valueT1, 498 } 499 } 500 return sigma; 501 } 502 503 //....oooOO0OOooo........oooOO0OOooo........oo 504 505 G4double G4DNARPWBAIonisationModel::Interpolat 506 507 508 509 510 { 511 G4double value = 0.; 512 513 // Log-log interpolation by default 514 515 if(e1 != 0 && e2 != 0 && (std::log10(e2) - s 516 !fasterCode) 517 { 518 G4double a = 519 (std::log10(xs2) - std::log10(xs1)) / (s 520 G4double b = std::log10(xs2) - a * std 521 G4double sigma = a * std::log10(e) + b; 522 value = (std::pow(10., sigma)); 523 } 524 // Switch to lin-lin interpolation 525 /* 526 if ((e2-e1)!=0) 527 { 528 G4double d1 = xs1; 529 G4double d2 = xs2; 530 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1) 531 } 532 */ 533 534 // Switch to log-lin interpolation for faste 535 if((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && 536 { 537 G4double d1 = std::log10(xs1); 538 G4double d2 = std::log10(xs2); 539 value = std::pow(10., (d1 + (d2 - d1 540 } 541 542 // Switch to lin-lin interpolation for faste 543 // in case one of xs1 or xs2 (=cum proba) va 544 545 if((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) 546 { 547 G4double d1 = xs1; 548 G4double d2 = xs2; 549 value = (d1 + (d2 - d1) * (e - e1) / 550 } 551 return value; 552 } 553 554 //....oooOO0OOooo........oooOO0OOooo........oo 555 556 G4double G4DNARPWBAIonisationModel::QuadInterp 557 const G4double& e11, const G4double& e12, co 558 const G4double& e22, const G4double& xs11, c 559 const G4double& xs21, const G4double& xs22, 560 const G4double& t2, const G4double& t, const 561 { 562 G4double interpolatedvalue1 = Interpolate(e1 563 G4double interpolatedvalue2 = Interpolate(e2 564 G4double value = 565 Interpolate(t1, t2, t, interpolatedvalue1, 566 567 return value; 568 } 569 570 G4double G4DNARPWBAIonisationModel::GetPartial 571 const G4Material* /*material*/, G4int level, 572 const G4ParticleDefinition* particle, G4doub 573 { 574 if(fpTotalCrossSection != nullptr && particl 575 { 576 G4Exception("G4DNARPWBAIonisationModel::Ge 577 FatalException, "Model not app 578 } 579 return fpTotalCrossSection->GetComponent(lev 580 } 581 582 //....oooOO0OOooo........oooOO0OOooo........oo 583 584 G4int G4DNARPWBAIonisationModel::RandomSelect( 585 { 586 if(fpTotalCrossSection == nullptr) 587 { 588 G4Exception("G4DNARPWBAIonisationModel::Ra 589 FatalException, "Model not app 590 } 591 else 592 { 593 auto valuesBuffer = new G4double[fpTotalCr 594 const auto n = (G4int)fpTotalCrossSection 595 G4int i(n); 596 G4double value = 0.; 597 while(i > 0) 598 { 599 --i; 600 valuesBuffer[i] = fpTotalCrossSection->G 601 value += valuesBuffer[i]; 602 } 603 value *= G4UniformRand(); 604 i = n; 605 606 while(i > 0) 607 { 608 --i; 609 if(valuesBuffer[i] > value) 610 { 611 delete[] valuesBuffer; 612 return i; 613 } 614 value -= valuesBuffer[i]; 615 } 616 delete[] valuesBuffer; 617 } 618 return 0; 619 } 620 621 //....oooOO0OOooo........oooOO0OOooo........oo 622 623 G4double 624 G4DNARPWBAIonisationModel::RandomizeEjectedEle 625 const G4double& k, const G4int& shell) 626 { 627 G4double random = G4 628 G4double secondaryKineticEnergy = 629 TransferedEnergy(k / eV, shell, random) * 630 waterStructure.IonisationEnergy(shell); 631 if(secondaryKineticEnergy < 0.) 632 { 633 return 0.; 634 } 635 return secondaryKineticEnergy; 636 } 637 638 //....oooOO0OOooo........oooOO0OOooo........oo 639 640 G4double G4DNARPWBAIonisationModel::Transfered 641 642 643 { 644 G4double nrj = 0.; 645 G4double valueK1 = 0; 646 G4double valueK2 = 0; 647 G4double valuePROB21 = 0; 648 G4double valuePROB22 = 0; 649 G4double valuePROB12 = 0; 650 G4double valuePROB11 = 0; 651 G4double nrjTransf11 = 0; 652 G4double nrjTransf12 = 0; 653 G4double nrjTransf21 = 0; 654 G4double nrjTransf22 = 0; 655 // Protection against out of boundary access 656 if(k == pTdummyVec.back()) 657 { 658 k = k * (1. - 1e-12); 659 } 660 // k should be in eV 661 662 auto k2 = std::upper_bound(pTdummyVec.begin( 663 auto k1 = k2 - 1; 664 665 // SI : the following condition avoids situa 666 // element, 667 // for eg. when the last element is zer 668 if(random <= pProbaShellMap[ionizationLevelI 669 random <= pProbaShellMap[ionizationLevelI 670 { 671 auto prob12 = std::upper_bound( 672 pProbaShellMap[ionizationLevelIndex][(*k 673 pProbaShellMap[ionizationLevelIndex][(*k 674 auto prob11 = prob12 - 1; 675 auto prob22 = std::upper_bound( 676 pProbaShellMap[ionizationLevelIndex][(*k 677 pProbaShellMap[ionizationLevelIndex][(*k 678 679 auto prob21 = prob22 - 1; 680 681 valueK1 = *k1; 682 valueK2 = *k2; 683 valuePROB21 = *prob21; 684 valuePROB22 = *prob22; 685 valuePROB12 = *prob12; 686 valuePROB11 = *prob11; 687 688 nrjTransf11 = pNrjTransfData[ionizationLev 689 nrjTransf12 = pNrjTransfData[ionizationLev 690 nrjTransf21 = pNrjTransfData[ionizationLev 691 nrjTransf22 = pNrjTransfData[ionizationLev 692 } 693 694 // Avoids cases where cum xs is zero for k1 695 // k1<k2) 696 697 if(random > pProbaShellMap[ionizationLevelIn 698 { 699 auto prob22 = std::upper_bound( 700 pProbaShellMap[ionizationLevelIndex][(*k 701 pProbaShellMap[ionizationLevelIndex][(*k 702 auto prob21 = prob22 - 1; 703 704 valueK1 = *k1; 705 valueK2 = *k2; 706 valuePROB21 = *prob21; 707 valuePROB22 = *prob22; 708 nrjTransf21 = pNrjTransfData[ionizationLev 709 nrjTransf22 = pNrjTransfData[ionizationLev 710 711 G4double interpolatedvalue2 = 712 Interpolate(valuePROB21, valuePROB22, ra 713 714 G4double value = Interpolate(valueK1, valu 715 return value; 716 } 717 G4double nrjTransfProduct = 718 nrjTransf11 * nrjTransf12 * nrjTransf21 * 719 720 if(nrjTransfProduct != 0.) 721 { 722 nrj = QuadInterpolator(valuePROB11, valueP 723 nrjTransf11, nrjTra 724 valueK1, valueK2, k 725 } 726 return nrj; 727 } 728