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 "G4DNABornIonisationModel1.hh" 29 #include "G4PhysicalConstants.hh" 30 #include "G4SystemOfUnits.hh" 31 #include "G4UAtomicDeexcitation.hh" 32 #include "G4LossTableManager.hh" 33 #include "G4DNAChemistryManager.hh" 34 #include "G4DNAMolecularMaterial.hh" 35 #include "G4DNABornAngle.hh" 36 #include "G4DeltaAngle.hh" 37 #include "G4Exp.hh" 38 39 //....oooOO0OOooo........oooOO0OOooo........oo 40 41 using namespace std; 42 43 //....oooOO0OOooo........oooOO0OOooo........oo 44 45 G4DNABornIonisationModel1::G4DNABornIonisation 46 47 G4VEmModel(nam) 48 { 49 verboseLevel = 0; 50 // Verbosity scale: 51 // 0 = nothing 52 // 1 = warning for energy non-conservation 53 // 2 = details of energy budget 54 // 3 = calculation of cross sections, file o 55 // 4 = entering in methods 56 57 if (verboseLevel > 0) 58 { 59 G4cout << "Born ionisation model is constr 60 } 61 62 // Mark this model as "applicable" for atomi 63 SetDeexcitationFlag(true); 64 fAtomDeexcitation = nullptr; 65 fParticleChangeForGamma = nullptr; 66 fpMolWaterDensity = nullptr; 67 68 // Define default angular generator 69 SetAngularDistribution(new G4DNABornAngle()) 70 71 // Selection of computation method 72 73 fasterCode = false; 74 75 // Selection of stationary mode 76 77 statCode = false; 78 79 // Selection of SP scaling 80 81 spScaling = true; 82 } 83 84 //....oooOO0OOooo........oooOO0OOooo........oo 85 86 G4DNABornIonisationModel1::~G4DNABornIonisatio 87 { 88 // Cross section 89 90 std::map<G4String, G4DNACrossSectionDataSet* 91 for (pos = tableData.begin(); pos != tableDa 92 { 93 G4DNACrossSectionDataSet* table = pos->sec 94 delete table; 95 } 96 97 // Final state 98 99 eVecm.clear(); 100 pVecm.clear(); 101 } 102 103 //....oooOO0OOooo........oooOO0OOooo........oo 104 105 void G4DNABornIonisationModel1::Initialise(con 106 con 107 { 108 109 if (verboseLevel > 3) 110 { 111 G4cout << "Calling G4DNABornIonisationMode 112 } 113 114 // Energy limits 115 116 G4String fileElectron("dna/sigma_ionisation_ 117 G4String fileProton("dna/sigma_ionisation_p_ 118 119 G4ParticleDefinition* electronDef = G4Electr 120 G4ParticleDefinition* protonDef = G4Proton:: 121 122 G4String electron; 123 G4String proton; 124 125 G4double scaleFactor = (1.e-22 / 3.343) * m* 126 127 const char *path = G4FindDataDir("G4LEDATA") 128 129 // *** ELECTRON 130 131 electron = electronDef->GetParticleName(); 132 133 tableFile[electron] = fileElectron; 134 135 lowEnergyLimit[electron] = 11. * eV; 136 highEnergyLimit[electron] = 1. * MeV; 137 138 // Cross section 139 140 auto tableE = new G4DNACrossSectionDataSet( 141 tableE->LoadData(fileElectron); 142 143 tableData[electron] = tableE; 144 145 // Final state 146 147 std::ostringstream eFullFileName; 148 149 if (fasterCode) eFullFileName << path << "/d 150 if (!fasterCode) eFullFileName << path << "/ 151 152 std::ifstream eDiffCrossSection(eFullFileNam 153 154 if (!eDiffCrossSection) 155 { 156 if (fasterCode) G4Exception("G4DNABornIoni 157 FatalException,"Missing data file:/dna 158 159 if (!fasterCode) G4Exception("G4DNABornIon 160 FatalException,"Missing data file:/dna 161 } 162 163 // Clear the arrays for re-initialization ca 164 // March 25th, 2014 - Vaclav Stepan, Sebasti 165 166 eTdummyVec.clear(); 167 pTdummyVec.clear(); 168 169 eVecm.clear(); 170 pVecm.clear(); 171 172 for (G4int j=0; j<5; j++) 173 { 174 eProbaShellMap[j].clear(); 175 pProbaShellMap[j].clear(); 176 177 eDiffCrossSectionData[j].clear(); 178 pDiffCrossSectionData[j].clear(); 179 180 eNrjTransfData[j].clear(); 181 pNrjTransfData[j].clear(); 182 } 183 184 // 185 186 eTdummyVec.push_back(0.); 187 while(!eDiffCrossSection.eof()) 188 { 189 G4double tDummy; 190 G4double eDummy; 191 eDiffCrossSection>>tDummy>>eDummy; 192 if (tDummy != eTdummyVec.back()) eTdummyVe 193 194 G4double tmp; 195 for (G4int j=0; j<5; j++) 196 { 197 eDiffCrossSection>> tmp; 198 199 eDiffCrossSectionData[j][tDummy][eDummy] 200 201 if (fasterCode) 202 { 203 eNrjTransfData[j][tDummy][eDiffCrossSe 204 eProbaShellMap[j][tDummy].push_back(eD 205 } 206 207 // SI - only if eof is not reached 208 if (!eDiffCrossSection.eof() && !fasterC 209 210 if (!fasterCode) eVecm[tDummy].push_back 211 212 } 213 } 214 215 // *** PROTON 216 217 proton = protonDef->GetParticleName(); 218 219 tableFile[proton] = fileProton; 220 221 lowEnergyLimit[proton] = 500. * keV; 222 highEnergyLimit[proton] = 100. * MeV; 223 224 // Cross section 225 226 auto tableP = new G4DNACrossSectionDataSet( 227 tableP->LoadData(fileProton); 228 229 tableData[proton] = tableP; 230 231 // Final state 232 233 std::ostringstream pFullFileName; 234 235 if (fasterCode) pFullFileName << path << "/d 236 237 if (!fasterCode) pFullFileName << path << "/ 238 239 std::ifstream pDiffCrossSection(pFullFileNam 240 241 if (!pDiffCrossSection) 242 { 243 if (fasterCode) G4Exception("G4DNABornIoni 244 FatalException,"Missing data file:/dna 245 246 if (!fasterCode) G4Exception("G4DNABornIon 247 FatalException,"Missing data file:/dna 248 } 249 250 pTdummyVec.push_back(0.); 251 while(!pDiffCrossSection.eof()) 252 { 253 G4double tDummy; 254 G4double eDummy; 255 pDiffCrossSection>>tDummy>>eDummy; 256 if (tDummy != pTdummyVec.back()) pTdummyVe 257 for (G4int j=0; j<5; j++) 258 { 259 pDiffCrossSection>>pDiffCrossSectionData 260 261 if (fasterCode) 262 { 263 pNrjTransfData[j][tDummy][pDiffCrossSe 264 pProbaShellMap[j][tDummy].push_back(pD 265 } 266 267 // SI - only if eof is not reached ! 268 if (!pDiffCrossSection.eof() && !fasterC 269 270 if (!fasterCode) pVecm[tDummy].push_back 271 } 272 } 273 274 // 275 276 if (particle==electronDef) 277 { 278 SetLowEnergyLimit(lowEnergyLimit[electron] 279 SetHighEnergyLimit(highEnergyLimit[electro 280 } 281 282 if (particle==protonDef) 283 { 284 SetLowEnergyLimit(lowEnergyLimit[proton]); 285 SetHighEnergyLimit(highEnergyLimit[proton] 286 } 287 288 if( verboseLevel>0 ) 289 { 290 G4cout << "Born ionisation model is initia 291 << "Energy range: " 292 << LowEnergyLimit() / eV << " eV - " 293 << HighEnergyLimit() / keV << " keV for " 294 << particle->GetParticleName() 295 << G4endl; 296 } 297 298 // Initialize water density pointer 299 300 fpMolWaterDensity = G4DNAMolecularMaterial:: 301 GetNumMolPerVolTableFor(G4Material::GetMater 302 303 // AD 304 305 fAtomDeexcitation = G4LossTableManager::Inst 306 307 // 308 309 if (isInitialised) 310 { return;} 311 fParticleChangeForGamma = GetParticleChangeF 312 isInitialised = true; 313 } 314 315 //....oooOO0OOooo........oooOO0OOooo........oo 316 317 G4double G4DNABornIonisationModel1::CrossSecti 318 319 320 321 322 { 323 if (verboseLevel > 3) 324 { 325 G4cout << "Calling CrossSectionPerVolume() 326 << G4endl; 327 } 328 329 if ( 330 particleDefinition != G4Proton::ProtonDe 331 && 332 particleDefinition != G4Electron::Electr 333 ) 334 335 return 0; 336 337 // Calculate total cross section for model 338 339 G4double lowLim = 0; 340 G4double highLim = 0; 341 G4double sigma=0; 342 343 G4double waterDensity = (*fpMolWaterDensity) 344 345 const G4String& particleName = particleDefin 346 347 std::map< G4String,G4double,std::less<G4Stri 348 pos1 = lowEnergyLimit.find(particleName); 349 if (pos1 != lowEnergyLimit.end()) 350 { 351 lowLim = pos1->second; 352 } 353 354 std::map< G4String,G4double,std::less<G4Stri 355 pos2 = highEnergyLimit.find(particleName); 356 if (pos2 != highEnergyLimit.end()) 357 { 358 highLim = pos2->second; 359 } 360 361 if (ekin >= lowLim && ekin <= highLim) 362 { 363 std::map< G4String,G4DNACrossSectionDataSe 364 pos = tableData.find(particleName); 365 366 if (pos != tableData.end()) 367 { 368 G4DNACrossSectionDataSet* table = pos->s 369 if (table != nullptr) 370 { 371 sigma = table->FindValue(ekin); 372 373 // ICRU49 electronic SP scaling - ZF, 374 375 if (particleDefinition == G4Proton::Pr 376 { 377 G4double A = 1.39241700556072800000E- 378 G4double B = -8.52610412942622630000E 379 sigma = sigma * G4Exp(A*(ekin/eV)+B); 380 } 381 // 382 383 } 384 } 385 else 386 { 387 G4Exception("G4DNABornIonisationModel1:: 388 FatalException,"Model not applicable 389 } 390 } 391 392 if (verboseLevel > 2) 393 { 394 G4cout << "_______________________________ 395 G4cout << "G4DNABornIonisationModel1 - XS 396 G4cout << "Kinetic energy(eV)=" << ekin/eV 397 G4cout << "Cross section per water molecul 398 G4cout << "Cross section per water molecul 399 G4cout << "G4DNABornIonisationModel1 - XS 400 } 401 402 return sigma*waterDensity; 403 } 404 405 //....oooOO0OOooo........oooOO0OOooo........oo 406 407 void G4DNABornIonisationModel1::SampleSecondar 408 409 410 411 412 { 413 414 if (verboseLevel > 3) 415 { 416 G4cout << "Calling SampleSecondaries() of 417 << G4endl; 418 } 419 420 G4double lowLim = 0; 421 G4double highLim = 0; 422 423 G4double k = particle->GetKineticEnergy(); 424 425 const G4String& particleName = particle->Get 426 427 std::map< G4String,G4double,std::less<G4Stri 428 pos1 = lowEnergyLimit.find(particleName); 429 430 if (pos1 != lowEnergyLimit.end()) 431 { 432 lowLim = pos1->second; 433 } 434 435 std::map< G4String,G4double,std::less<G4Stri 436 pos2 = highEnergyLimit.find(particleName); 437 438 if (pos2 != highEnergyLimit.end()) 439 { 440 highLim = pos2->second; 441 } 442 443 if (k >= lowLim && k <= highLim) 444 { 445 G4ParticleMomentum primaryDirection = part 446 G4double particleMass = particle->GetDefin 447 G4double totalEnergy = k + particleMass; 448 G4double pSquare = k * (totalEnergy + part 449 G4double totalMomentum = std::sqrt(pSquare 450 451 G4int ionizationShell = 0; 452 453 if (!fasterCode) ionizationShell = RandomS 454 455 // SI: The following protection is necessa 456 // sigmadiff_ionisation_e_born.dat has no 457 // sigmadiff_cumulated_ionisation_e_born. 458 // this is due to the fact that the max a 459 // strictly above this value have non zer 460 461 if (fasterCode) 462 do 463 { 464 ionizationShell = RandomSelect(k,particl 465 } while (k<19*eV && ionizationShell==2 && 466 467 G4double bindingEnergy = 0; 468 bindingEnergy = waterStructure.IonisationE 469 470 // SI: additional protection if tcs interp 471 if (k<bindingEnergy) return; 472 // 473 474 G4double secondaryKinetic=-1000*eV; 475 476 if (!fasterCode) 477 { 478 secondaryKinetic = RandomizeEjectedElect 479 } 480 else 481 { 482 secondaryKinetic = RandomizeEjectedElect 483 } 484 // 485 486 G4int Z = 8; 487 488 G4ThreeVector deltaDirection = 489 GetAngularDistribution()->SampleDirectionF 490 Z, ionizationShell, 491 couple->GetMaterial()); 492 493 if (secondaryKinetic>0) 494 { 495 auto dp = new G4DynamicParticle (G4Elec 496 fvect->push_back(dp); 497 } 498 499 if (particle->GetDefinition() == G4Electro 500 { 501 G4double deltaTotalMomentum = std::sqrt( 502 503 G4double finalPx = totalMomentum*primary 504 G4double finalPy = totalMomentum*primary 505 G4double finalPz = totalMomentum*primary 506 G4double finalMomentum = std::sqrt(final 507 finalPx /= finalMomentum; 508 finalPy /= finalMomentum; 509 finalPz /= finalMomentum; 510 511 G4ThreeVector direction; 512 direction.set(finalPx,finalPy,finalPz); 513 514 fParticleChangeForGamma->ProposeMomentum 515 } 516 517 else fParticleChangeForGamma->ProposeMomen 518 519 // AM: sample deexcitation 520 // here we assume that H_{2}O electronic l 521 // this can be considered true with a roug 522 523 std::size_t secNumberInit = 0;// need to k 524 std::size_t secNumberFinal = 0;// So I'll 525 526 G4double scatteredEnergy = k-bindingEnergy 527 528 // SI: only atomic deexcitation from K she 529 if((fAtomDeexcitation != nullptr) && ioniz 530 { 531 const G4AtomicShell* shell = 532 fAtomDeexcitation->GetAtomicShell(Z, G 533 secNumberInit = fvect->size(); 534 fAtomDeexcitation->GenerateParticles(fve 535 secNumberFinal = fvect->size(); 536 537 //TEST 538 //G4cout << "ionizationShell=" << ioniza 539 //G4cout << "bindingEnergy=" << bindingE 540 541 if(secNumberFinal > secNumberInit) 542 { 543 for (std::size_t i=secNumberInit; i<secNumbe 544 { 545 //Check if there is enough residual 546 if (bindingEnergy >= ((*fvect)[i])-> 547 { 548 //Ok, this is a valid secondary: 549 bindingEnergy -= ((*fvect)[i])->GetKine 550 //G4cout << "--deex nrj=" << ((*fvect)[ 551 //<< G4endl; 552 } 553 else 554 { 555 //Invalid secondary: not enough energy 556 //Keep its energy in the local deposit 557 delete (*fvect)[i]; 558 (*fvect)[i]=nullptr; 559 } 560 } 561 } 562 563 //TEST 564 //G4cout << "k=" << k/eV<< G4endl; 565 //G4cout << "secondaryKinetic=" << seconda 566 //G4cout << "scatteredEnergy=" << scattere 567 //G4cout << "deposited energy=" << binding 568 // 569 570 } 571 572 //This should never happen 573 if(bindingEnergy < 0.0) 574 G4Exception("G4DNABornIonisatioModel1::Sa 575 "em2050",FatalException,"Nega 576 577 //bindingEnergy has been decreased 578 //by the amount of energy taken away by de 579 if (!statCode) 580 { 581 fParticleChangeForGamma->SetProposedKine 582 fParticleChangeForGamma->ProposeLocalEne 583 } 584 else 585 { 586 fParticleChangeForGamma->SetProposedKine 587 fParticleChangeForGamma->ProposeLocalEne 588 } 589 590 // TEST ////////////////////////// 591 // if (secondaryKinetic<0) abort(); 592 // if (scatteredEnergy<0) abort(); 593 // if (k-scatteredEnergy-secondaryKinetic- 594 // if (k-scatteredEnergy<0) abort(); 595 ///////////////////////////////// 596 597 const G4Track * theIncomingTrack = fPartic 598 G4DNAChemistryManager::Instance()->CreateW 599 ionizationShell, 600 theIncomingTrack); 601 } 602 } 603 604 //....oooOO0OOooo........oooOO0OOooo........oo 605 606 G4double G4DNABornIonisationModel1::RandomizeE 607 608 609 { 610 // G4cout << "*** SLOW computation for " << 611 612 if (particleDefinition == G4Electron::Electr 613 { 614 G4double maximumEnergyTransfer = 0.; 615 if ((k + waterStructure.IonisationEnergy(s 616 maximumEnergyTransfer = k; 617 else 618 maximumEnergyTransfer = (k + waterStruct 619 620 // SI : original method 621 /* 622 G4double crossSectionMaximum = 0.; 623 for(G4double value=waterStructure.Ionisati 624 { 625 G4double differentialCrossSection = Diffe 626 if(differentialCrossSection >= crossSecti 627 } 628 */ 629 630 // SI : alternative method 631 G4double crossSectionMaximum = 0.; 632 633 G4double minEnergy = waterStructure.Ionisa 634 G4double maxEnergy = maximumEnergyTransfer 635 G4int nEnergySteps = 50; 636 637 G4double value(minEnergy); 638 G4double stpEnergy(std::pow(maxEnergy / va 639 1. / static_ca 640 G4int step(nEnergySteps); 641 while (step > 0) 642 { 643 step--; 644 G4double differentialCrossSection = 645 DifferentialCrossSection(particleDef 646 k / eV, 647 value / eV, 648 shell); 649 if (differentialCrossSection >= crossSec 650 crossSectionMaximum = differentialCros 651 value *= stpEnergy; 652 } 653 // 654 655 G4double secondaryElectronKineticEnergy = 656 do 657 { 658 secondaryElectronKineticEnergy = G4Unifo 659 } while(G4UniformRand()*crossSectionMaximu 660 DifferentialCrossSection(particleDefin 661 (secondaryElectronKineticEnergy+wa 662 663 return secondaryElectronKineticEnergy; 664 665 } 666 667 if (particleDefinition == G4Proton::ProtonDe 668 { 669 G4double maximumKineticEnergyTransfer = 4. 670 * (electron_mass_c2 / proton_mass_c2) 671 672 G4double crossSectionMaximum = 0.; 673 for (G4double value = waterStructure.Ionis 674 value <= 4. * waterStructure.Ionisatio 675 { 676 G4double differentialCrossSection = 677 DifferentialCrossSection(particleDef 678 k / eV, 679 value / eV, 680 shell); 681 if (differentialCrossSection >= crossSec 682 crossSectionMaximum = differentialCros 683 } 684 685 G4double secondaryElectronKineticEnergy = 686 do 687 { 688 secondaryElectronKineticEnergy = G4Unifo 689 } while(G4UniformRand()*crossSectionMaximu 690 DifferentialCrossSection(particleDefin 691 (secondaryElectronKineticEnergy+wa 692 693 return secondaryElectronKineticEnergy; 694 } 695 696 return 0; 697 } 698 699 //....oooOO0OOooo........oooOO0OOooo........oo 700 701 // The following section is not used anymore b 702 // GetAngularDistribution()->SampleDirectionFo 703 704 /* 705 void G4DNABornIonisationModel1::RandomizeEjec 706 G4double k, 707 G4double secKinetic, 708 G4double & cosTheta, 709 G4double & phi ) 710 { 711 if (particleDefinition == G4Electron::Electro 712 { 713 phi = twopi * G4UniformRand(); 714 if (secKinetic < 50.*eV) cosTheta = (2.*G4Uni 715 else if (secKinetic <= 200.*eV) 716 { 717 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4 718 else cosTheta = G4UniformRand()*(std::sqrt(2. 719 } 720 else 721 { 722 G4double sin2O = (1.-secKinetic/k) / (1.+secK 723 cosTheta = std::sqrt(1.-sin2O); 724 } 725 } 726 727 else if (particleDefinition == G4Proton::Prot 728 { 729 G4double maxSecKinetic = 4.* (electron_mass_c 730 phi = twopi * G4UniformRand(); 731 732 // cosTheta = std::sqrt(secKinetic / maxSecKi 733 734 // Restriction below 100 eV from Emfietzoglou 735 736 if (secKinetic>100*eV) cosTheta = std::sqrt(s 737 else cosTheta = (2.*G4UniformRand())-1.; 738 739 } 740 } 741 */ 742 743 //....oooOO0OOooo........oooOO0OOooo........oo 744 G4double G4DNABornIonisationModel1::Differenti 745 746 747 748 { 749 G4double sigma = 0.; 750 751 if (energyTransfer >= waterStructure.Ionisat 752 { 753 G4double valueT1 = 0; 754 G4double valueT2 = 0; 755 G4double valueE21 = 0; 756 G4double valueE22 = 0; 757 G4double valueE12 = 0; 758 G4double valueE11 = 0; 759 760 G4double xs11 = 0; 761 G4double xs12 = 0; 762 G4double xs21 = 0; 763 G4double xs22 = 0; 764 765 if (particleDefinition == G4Electron::Elec 766 { 767 768 // Protection against out of boundary ac 769 if (k==eTdummyVec.back()) k=k*(1.-1e-12) 770 // 771 772 // k should be in eV and energy transfer 773 774 auto t2 = std::upper_bound(eTdummyVec.be 775 776 777 778 auto t1 = t2 - 1; 779 780 // SI : the following condition avoids s 781 if (energyTransfer <= eVecm[(*t1)].back( 782 && energyTransfer <= eVecm[(*t2)].ba 783 { 784 auto e12 = 785 std::upper_bound(eVecm[(*t1)].begi 786 eVecm[(*t1)].end( 787 energyTransfer); 788 auto e11 = e12 - 1; 789 790 auto e22 = 791 std::upper_bound(eVecm[(*t2)].begi 792 eVecm[(*t2)].end( 793 energyTransfer); 794 auto e21 = e22 - 1; 795 796 valueT1 = *t1; 797 valueT2 = *t2; 798 valueE21 = *e21; 799 valueE22 = *e22; 800 valueE12 = *e12; 801 valueE11 = *e11; 802 803 xs11 = eDiffCrossSectionData[ionizatio 804 xs12 = eDiffCrossSectionData[ionizatio 805 xs21 = eDiffCrossSectionData[ionizatio 806 xs22 = eDiffCrossSectionData[ionizatio 807 808 } 809 810 } 811 812 if (particleDefinition == G4Proton::Proton 813 { 814 // Protection against out of boundary ac 815 if (k==pTdummyVec.back()) k=k*(1.-1e-12) 816 // 817 818 // k should be in eV and energy transfer 819 auto t2 = std::upper_bound(pTdummyVec.be 820 821 822 auto t1 = t2 - 1; 823 824 auto e12 = std::upper_bound(pVecm[(*t1)] 825 826 827 auto e11 = e12 - 1; 828 829 auto e22 = std::upper_bound(pVecm[(*t2)] 830 831 832 auto e21 = e22 - 1; 833 834 valueT1 = *t1; 835 valueT2 = *t2; 836 valueE21 = *e21; 837 valueE22 = *e22; 838 valueE12 = *e12; 839 valueE11 = *e11; 840 841 xs11 = pDiffCrossSectionData[ionizationL 842 xs12 = pDiffCrossSectionData[ionizationL 843 xs21 = pDiffCrossSectionData[ionizationL 844 xs22 = pDiffCrossSectionData[ionizationL 845 846 } 847 848 G4double xsProduct = xs11 * xs12 * xs21 * 849 if (xsProduct != 0.) 850 { 851 sigma = QuadInterpolator(valueE11, 852 valueE12, 853 valueE21, 854 valueE22, 855 xs11, 856 xs12, 857 xs21, 858 xs22, 859 valueT1, 860 valueT2, 861 k, 862 energyTransfer) 863 } 864 } 865 866 return sigma; 867 } 868 869 //....oooOO0OOooo........oooOO0OOooo........oo 870 871 G4double G4DNABornIonisationModel1::Interpolat 872 873 874 875 876 { 877 G4double value = 0.; 878 879 // Log-log interpolation by default 880 881 if (e1 != 0 && e2 != 0 && (std::log10(e2) - 882 && !fasterCode) 883 { 884 G4double a = (std::log10(xs2) - std::log10 885 / (std::log10(e2) - std::log10(e1)); 886 G4double b = std::log10(xs2) - a * std::lo 887 G4double sigma = a * std::log10(e) + b; 888 value = (std::pow(10., sigma)); 889 } 890 891 // Switch to lin-lin interpolation 892 /* 893 if ((e2-e1)!=0) 894 { 895 G4double d1 = xs1; 896 G4double d2 = xs2; 897 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1) 898 } 899 */ 900 901 // Switch to log-lin interpolation for faste 902 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 & 903 { 904 G4double d1 = std::log10(xs1); 905 G4double d2 = std::log10(xs2); 906 value = std::pow(10., (d1 + (d2 - d1) * (e 907 } 908 909 // Switch to lin-lin interpolation for faste 910 // in case one of xs1 or xs2 (=cum proba) va 911 912 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) 913 { 914 G4double d1 = xs1; 915 G4double d2 = xs2; 916 value = (d1 + (d2 - d1) * (e - e1) / (e2 - 917 } 918 919 /* 920 G4cout 921 << e1 << " " 922 << e2 << " " 923 << e << " " 924 << xs1 << " " 925 << xs2 << " " 926 << value 927 << G4endl; 928 */ 929 930 return value; 931 } 932 933 //....oooOO0OOooo........oooOO0OOooo........oo 934 935 G4double G4DNABornIonisationModel1::QuadInterp 936 937 938 939 940 941 942 943 944 945 946 947 { 948 G4double interpolatedvalue1 = Interpolate(e1 949 G4double interpolatedvalue2 = Interpolate(e2 950 G4double value = Interpolate(t1, 951 t2, 952 t, 953 interpolatedval 954 interpolatedval 955 956 return value; 957 } 958 959 G4double G4DNABornIonisationModel1::GetPartial 960 961 962 963 { 964 std::map<G4String, G4DNACrossSectionDataSet* 965 pos = tableData.find(particle->GetParticleNa 966 967 if (pos != tableData.end()) 968 { 969 G4DNACrossSectionDataSet* table = pos->sec 970 return table->GetComponent(level)->FindVal 971 } 972 973 return 0; 974 } 975 976 //....oooOO0OOooo........oooOO0OOooo........oo 977 978 G4int G4DNABornIonisationModel1::RandomSelect( 979 980 { 981 G4int level = 0; 982 983 std::map<G4String, G4DNACrossSectionDataSet* 984 pos = tableData.find(particle); 985 986 if (pos != tableData.end()) 987 { 988 G4DNACrossSectionDataSet* table = pos->sec 989 990 if (table != nullptr) 991 { 992 auto valuesBuffer = new G4double[table- 993 const auto n = (G4int)table->NumberOfCo 994 G4int i(n); 995 G4double value = 0.; 996 997 while (i > 0) 998 { 999 i--; 1000 valuesBuffer[i] = table->GetComponent 1001 value += valuesBuffer[i]; 1002 } 1003 1004 value *= G4UniformRand(); 1005 1006 i = n; 1007 1008 while (i > 0) 1009 { 1010 i--; 1011 1012 if (valuesBuffer[i] > value) 1013 { 1014 delete[] valuesBuffer; 1015 return i; 1016 } 1017 value -= valuesBuffer[i]; 1018 } 1019 1020 1021 delete[] valuesBuffer; 1022 1023 } 1024 } else 1025 { 1026 G4Exception("G4DNABornIonisationModel1::R 1027 "em0002", 1028 FatalException, 1029 "Model not applicable to part 1030 } 1031 1032 return level; 1033 } 1034 1035 //....oooOO0OOooo........oooOO0OOooo........o 1036 1037 G4double G4DNABornIonisationModel1::Randomize 1038 1039 1040 { 1041 //G4cout << "*** FAST computation for " << 1042 1043 G4double secondaryElectronKineticEnergy = 0 1044 1045 G4double random = G4UniformRand(); 1046 1047 secondaryElectronKineticEnergy = Transfered 1048 1049 1050 1051 - waterStructure.IonisationEnergy(shell 1052 1053 //G4cout << RandomTransferedEnergy(particle 1054 1055 if (secondaryElectronKineticEnergy < 0.) 1056 return 0.; 1057 // 1058 1059 return secondaryElectronKineticEnergy; 1060 } 1061 1062 //....oooOO0OOooo........oooOO0OOooo........o 1063 1064 G4double G4DNABornIonisationModel1::Transfere 1065 1066 1067 1068 { 1069 G4double nrj = 0.; 1070 1071 G4double valueK1 = 0; 1072 G4double valueK2 = 0; 1073 G4double valuePROB21 = 0; 1074 G4double valuePROB22 = 0; 1075 G4double valuePROB12 = 0; 1076 G4double valuePROB11 = 0; 1077 1078 G4double nrjTransf11 = 0; 1079 G4double nrjTransf12 = 0; 1080 G4double nrjTransf21 = 0; 1081 G4double nrjTransf22 = 0; 1082 1083 if (particleDefinition == G4Electron::Elect 1084 { 1085 // Protection against out of boundary acc 1086 if (k==eTdummyVec.back()) k=k*(1.-1e-12); 1087 // 1088 1089 // k should be in eV 1090 auto k2 = std::upper_bound(eTdummyVec.beg 1091 1092 1093 auto k1 = k2 - 1; 1094 1095 /* 1096 G4cout << "----> k=" << k 1097 << " " << *k1 1098 << " " << *k2 1099 << " " << random 1100 << " " << ionizationLevelIndex 1101 << " " << eProbaShellMap[ionizationLevel 1102 << " " << eProbaShellMap[ionizationLevel 1103 << G4endl; 1104 */ 1105 1106 // SI : the following condition avoids si 1107 if (random <= eProbaShellMap[ionizationLe 1108 && random <= eProbaShellMap[ionizatio 1109 { 1110 auto prob12 = 1111 std::upper_bound(eProbaShellMap[ion 1112 eProbaShellMap[ion 1113 random); 1114 1115 auto prob11 = prob12 - 1; 1116 1117 auto prob22 = 1118 std::upper_bound(eProbaShellMap[ion 1119 eProbaShellMap[ion 1120 random); 1121 1122 auto prob21 = prob22 - 1; 1123 1124 valueK1 = *k1; 1125 valueK2 = *k2; 1126 valuePROB21 = *prob21; 1127 valuePROB22 = *prob22; 1128 valuePROB12 = *prob12; 1129 valuePROB11 = *prob11; 1130 1131 /* 1132 G4cout << " " << random << " " 1133 << valuePROB12 << " " << valuePROB21 < 1134 */ 1135 1136 nrjTransf11 = eNrjTransfData[ionization 1137 nrjTransf12 = eNrjTransfData[ionization 1138 nrjTransf21 = eNrjTransfData[ionization 1139 nrjTransf22 = eNrjTransfData[ionization 1140 1141 /* 1142 G4cout << " " << ionizationLeve 1143 << random << " " <<valueK1 << " " << v 1144 1145 G4cout << " " << random << " " 1146 << nrjTransf12 << " " << nrjTransf21 < 1147 */ 1148 1149 } 1150 1151 // Avoids cases where cum xs is zero for 1152 if (random > eProbaShellMap[ionizationLev 1153 { 1154 auto prob22 = 1155 std::upper_bound(eProbaShellMap[ion 1156 eProbaShellMap[ion 1157 random); 1158 1159 auto prob21 = prob22 - 1; 1160 1161 valueK1 = *k1; 1162 valueK2 = *k2; 1163 valuePROB21 = *prob21; 1164 valuePROB22 = *prob22; 1165 1166 //G4cout << " " << random << " " 1167 1168 nrjTransf21 = eNrjTransfData[ionization 1169 nrjTransf22 = eNrjTransfData[ionization 1170 1171 G4double interpolatedvalue2 = Interpola 1172 1173 1174 1175 1176 1177 // zeros are explicitly set 1178 1179 G4double value = Interpolate(valueK1, v 1180 1181 /* 1182 G4cout << " " << ionizationLeve 1183 << random << " " <<valueK1 << " " << v 1184 1185 G4cout << " " << random << " " 1186 << nrjTransf12 << " " << nrjTransf21 < 1187 1188 G4cout << "ici" << " " << value << G4e 1189 */ 1190 1191 return value; 1192 } 1193 } 1194 // 1195 else if (particleDefinition == G4Proton::Pr 1196 { 1197 // Protection against out of boundary acc 1198 if (k==pTdummyVec.back()) k=k*(1.-1e-12); 1199 // 1200 1201 // k should be in eV 1202 1203 auto k2 = std::upper_bound(pTdummyVec.beg 1204 1205 1206 1207 auto k1 = k2 - 1; 1208 1209 /* 1210 G4cout << "----> k=" << k 1211 << " " << *k1 1212 << " " << *k2 1213 << " " << random 1214 << " " << ionizationLevelIndex 1215 << " " << pProbaShellMap[ionizationLevel 1216 << " " << pProbaShellMap[ionizationLevel 1217 << G4endl; 1218 */ 1219 1220 // SI : the following condition avoids si 1221 // for eg. when the last element is 1222 if (random <= pProbaShellMap[ionizationLe 1223 && random <= pProbaShellMap[ionizatio 1224 { 1225 auto prob12 = 1226 std::upper_bound(pProbaShellMap[ion 1227 pProbaShellMap[ion 1228 random); 1229 1230 auto prob11 = prob12 - 1; 1231 1232 auto prob22 = 1233 std::upper_bound(pProbaShellMap[ion 1234 pProbaShellMap[ion 1235 random); 1236 1237 auto prob21 = prob22 - 1; 1238 1239 valueK1 = *k1; 1240 valueK2 = *k2; 1241 valuePROB21 = *prob21; 1242 valuePROB22 = *prob22; 1243 valuePROB12 = *prob12; 1244 valuePROB11 = *prob11; 1245 1246 /* 1247 G4cout << " " << random << " " 1248 << valuePROB12 << " " << valuePROB21 < 1249 */ 1250 1251 nrjTransf11 = pNrjTransfData[ionization 1252 nrjTransf12 = pNrjTransfData[ionization 1253 nrjTransf21 = pNrjTransfData[ionization 1254 nrjTransf22 = pNrjTransfData[ionization 1255 1256 /* 1257 G4cout << " " << ionizationLeve 1258 << random << " " <<valueK1 << " " << v 1259 1260 G4cout << " " << random << " " 1261 << nrjTransf12 << " " << nrjTransf21 < 1262 */ 1263 } 1264 1265 // Avoids cases where cum xs is zero for 1266 1267 if (random > pProbaShellMap[ionizationLev 1268 { 1269 auto prob22 = 1270 std::upper_bound(pProbaShellMap[ion 1271 pProbaShellMap[ion 1272 random); 1273 1274 auto prob21 = prob22 - 1; 1275 1276 valueK1 = *k1; 1277 valueK2 = *k2; 1278 valuePROB21 = *prob21; 1279 valuePROB22 = *prob22; 1280 1281 //G4cout << " " << random << " " 1282 1283 nrjTransf21 = pNrjTransfData[ionization 1284 nrjTransf22 = pNrjTransfData[ionization 1285 1286 G4double interpolatedvalue2 = Interpola 1287 1288 1289 1290 1291 1292 // zeros are explicitly set 1293 1294 G4double value = Interpolate(valueK1, v 1295 1296 /* 1297 G4cout << " " << ionizationLeve 1298 << random << " " <<valueK1 << " " << v 1299 1300 G4cout << " " << random << " " 1301 << nrjTransf12 << " " << nrjTransf21 < 1302 1303 G4cout << "ici" << " " << value << G4e 1304 */ 1305 1306 return value; 1307 } 1308 } 1309 // End electron and proton cases 1310 1311 G4double nrjTransfProduct = nrjTransf11 * n 1312 * nrjTransf22; 1313 //G4cout << "nrjTransfProduct=" << nrjTrans 1314 1315 if (nrjTransfProduct != 0.) 1316 { 1317 nrj = QuadInterpolator(valuePROB11, 1318 valuePROB12, 1319 valuePROB21, 1320 valuePROB22, 1321 nrjTransf11, 1322 nrjTransf12, 1323 nrjTransf21, 1324 nrjTransf22, 1325 valueK1, 1326 valueK2, 1327 k, 1328 random); 1329 } 1330 //G4cout << nrj << endl; 1331 1332 return nrj; 1333 } 1334