00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "G4DNABornIonisationModel.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "G4UAtomicDeexcitation.hh"
00033 #include "G4LossTableManager.hh"
00034 #include "G4DNAChemistryManager.hh"
00035 #include "G4DNAMolecularMaterial.hh"
00036
00037
00038
00039 using namespace std;
00040
00041
00042
00043 G4DNABornIonisationModel::G4DNABornIonisationModel(const G4ParticleDefinition*,
00044 const G4String& nam)
00045 :G4VEmModel(nam),isInitialised(false)
00046 {
00047
00048
00049 verboseLevel= 0;
00050
00051
00052
00053
00054
00055
00056
00057 if( verboseLevel>0 )
00058 {
00059 G4cout << "Born ionisation model is constructed " << G4endl;
00060 }
00061
00062
00063 SetDeexcitationFlag(true);
00064 fAtomDeexcitation = 0;
00065 fParticleChangeForGamma = 0;
00066 fpMolWaterDensity = 0;
00067 }
00068
00069
00070
00071 G4DNABornIonisationModel::~G4DNABornIonisationModel()
00072 {
00073
00074
00075 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00076 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
00077 {
00078 G4DNACrossSectionDataSet* table = pos->second;
00079 delete table;
00080 }
00081
00082
00083
00084 eVecm.clear();
00085 pVecm.clear();
00086
00087 }
00088
00089
00090
00091 void G4DNABornIonisationModel::Initialise(const G4ParticleDefinition* particle,
00092 const G4DataVector& )
00093 {
00094
00095 if (verboseLevel > 3)
00096 G4cout << "Calling G4DNABornIonisationModel::Initialise()" << G4endl;
00097
00098
00099
00100 G4String fileElectron("dna/sigma_ionisation_e_born");
00101 G4String fileProton("dna/sigma_ionisation_p_born");
00102
00103 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
00104 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
00105
00106 G4String electron;
00107 G4String proton;
00108
00109 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
00110
00111 char *path = getenv("G4LEDATA");
00112
00113
00114
00115 electron = electronDef->GetParticleName();
00116
00117 tableFile[electron] = fileElectron;
00118
00119 lowEnergyLimit[electron] = 11. * eV;
00120 highEnergyLimit[electron] = 1. * MeV;
00121
00122
00123
00124 G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
00125 tableE->LoadData(fileElectron);
00126
00127 tableData[electron] = tableE;
00128
00129
00130
00131 std::ostringstream eFullFileName;
00132 eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
00133 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
00134
00135 if (!eDiffCrossSection)
00136 {
00137 G4Exception("G4DNABornIonisationModel::Initialise","em0003",
00138 FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
00139 }
00140
00141 eTdummyVec.push_back(0.);
00142 while(!eDiffCrossSection.eof())
00143 {
00144 double tDummy;
00145 double eDummy;
00146 eDiffCrossSection>>tDummy>>eDummy;
00147 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
00148 for (int j=0; j<5; j++)
00149 {
00150 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
00151
00152
00153 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
00154
00155 eVecm[tDummy].push_back(eDummy);
00156
00157 }
00158 }
00159
00160
00161
00162 proton = protonDef->GetParticleName();
00163
00164 tableFile[proton] = fileProton;
00165
00166 lowEnergyLimit[proton] = 500. * keV;
00167 highEnergyLimit[proton] = 100. * MeV;
00168
00169
00170
00171 G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
00172 tableP->LoadData(fileProton);
00173
00174 tableData[proton] = tableP;
00175
00176
00177
00178 std::ostringstream pFullFileName;
00179 pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
00180 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
00181
00182 if (!pDiffCrossSection)
00183 {
00184 G4Exception("G4DNABornIonisationModel::Initialise","em0003",
00185 FatalException,"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
00186 }
00187
00188 pTdummyVec.push_back(0.);
00189 while(!pDiffCrossSection.eof())
00190 {
00191 double tDummy;
00192 double eDummy;
00193 pDiffCrossSection>>tDummy>>eDummy;
00194 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
00195 for (int j=0; j<5; j++)
00196 {
00197 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
00198
00199
00200 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
00201
00202 pVecm[tDummy].push_back(eDummy);
00203 }
00204 }
00205
00206
00207
00208 if (particle==electronDef)
00209 {
00210 SetLowEnergyLimit(lowEnergyLimit[electron]);
00211 SetHighEnergyLimit(highEnergyLimit[electron]);
00212 }
00213
00214 if (particle==protonDef)
00215 {
00216 SetLowEnergyLimit(lowEnergyLimit[proton]);
00217 SetHighEnergyLimit(highEnergyLimit[proton]);
00218 }
00219
00220 if( verboseLevel>0 )
00221 {
00222 G4cout << "Born ionisation model is initialized " << G4endl
00223 << "Energy range: "
00224 << LowEnergyLimit() / eV << " eV - "
00225 << HighEnergyLimit() / keV << " keV for "
00226 << particle->GetParticleName()
00227 << G4endl;
00228 }
00229
00230
00231 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00232
00233
00234 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00235
00236 if (isInitialised) { return; }
00237 fParticleChangeForGamma = GetParticleChangeForGamma();
00238 isInitialised = true;
00239 }
00240
00241
00242
00243 G4double G4DNABornIonisationModel::CrossSectionPerVolume(const G4Material* material,
00244 const G4ParticleDefinition* particleDefinition,
00245 G4double ekin,
00246 G4double,
00247 G4double)
00248 {
00249 if (verboseLevel > 3)
00250 G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel" << G4endl;
00251
00252 if (
00253 particleDefinition != G4Proton::ProtonDefinition()
00254 &&
00255 particleDefinition != G4Electron::ElectronDefinition()
00256 )
00257
00258 return 0;
00259
00260
00261
00262 G4double lowLim = 0;
00263 G4double highLim = 0;
00264 G4double sigma=0;
00265
00266 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
00267
00268 if(waterDensity!= 0.0)
00269
00270 {
00271 const G4String& particleName = particleDefinition->GetParticleName();
00272
00273 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
00274 pos1 = lowEnergyLimit.find(particleName);
00275 if (pos1 != lowEnergyLimit.end())
00276 {
00277 lowLim = pos1->second;
00278 }
00279
00280 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00281 pos2 = highEnergyLimit.find(particleName);
00282 if (pos2 != highEnergyLimit.end())
00283 {
00284 highLim = pos2->second;
00285 }
00286
00287 if (ekin >= lowLim && ekin < highLim)
00288 {
00289 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00290 pos = tableData.find(particleName);
00291
00292 if (pos != tableData.end())
00293 {
00294 G4DNACrossSectionDataSet* table = pos->second;
00295 if (table != 0)
00296 {
00297 sigma = table->FindValue(ekin);
00298 }
00299 }
00300 else
00301 {
00302 G4Exception("G4DNABornIonisationModel::CrossSectionPerVolume","em0002",
00303 FatalException,"Model not applicable to particle type.");
00304 }
00305 }
00306
00307 if (verboseLevel > 2)
00308 {
00309 G4cout << "__________________________________" << G4endl;
00310 G4cout << "°°° G4DNABornIonisationModel - XS INFO START" << G4endl;
00311 G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
00312 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
00313 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
00314
00315 G4cout << "°°° G4DNABornIonisationModel - XS INFO END" << G4endl;
00316 }
00317
00318 }
00319
00320
00321 return sigma*waterDensity;
00322 }
00323
00324
00325
00326 void G4DNABornIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00327 const G4MaterialCutsCouple* ,
00328 const G4DynamicParticle* particle,
00329 G4double,
00330 G4double)
00331 {
00332
00333 if (verboseLevel > 3)
00334 G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel" << G4endl;
00335
00336 G4double lowLim = 0;
00337 G4double highLim = 0;
00338
00339 G4double k = particle->GetKineticEnergy();
00340
00341 const G4String& particleName = particle->GetDefinition()->GetParticleName();
00342
00343 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
00344 pos1 = lowEnergyLimit.find(particleName);
00345
00346 if (pos1 != lowEnergyLimit.end())
00347 {
00348 lowLim = pos1->second;
00349 }
00350
00351 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00352 pos2 = highEnergyLimit.find(particleName);
00353
00354 if (pos2 != highEnergyLimit.end())
00355 {
00356 highLim = pos2->second;
00357 }
00358
00359 if (k >= lowLim && k < highLim)
00360 {
00361 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
00362 G4double particleMass = particle->GetDefinition()->GetPDGMass();
00363 G4double totalEnergy = k + particleMass;
00364 G4double pSquare = k * (totalEnergy + particleMass);
00365 G4double totalMomentum = std::sqrt(pSquare);
00366
00367 G4int ionizationShell = RandomSelect(k,particleName);
00368
00369
00370
00371
00372
00373 G4int secNumberInit = 0;
00374 G4int secNumberFinal = 0;
00375
00376 G4double bindingEnergy = 0;
00377 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
00378
00379 if(fAtomDeexcitation) {
00380 G4int Z = 8;
00381 G4AtomicShellEnumerator as = fKShell;
00382
00383 if (ionizationShell <5 && ionizationShell >1)
00384 {
00385 as = G4AtomicShellEnumerator(4-ionizationShell);
00386 }
00387 else if (ionizationShell <2)
00388 {
00389 as = G4AtomicShellEnumerator(3);
00390 }
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00402 secNumberInit = fvect->size();
00403 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
00404 secNumberFinal = fvect->size();
00405 }
00406
00407 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
00408
00409 G4double cosTheta = 0.;
00410 G4double phi = 0.;
00411 RandomizeEjectedElectronDirection(particle->GetDefinition(), k,secondaryKinetic, cosTheta, phi);
00412
00413 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
00414 G4double dirX = sinTheta*std::cos(phi);
00415 G4double dirY = sinTheta*std::sin(phi);
00416 G4double dirZ = cosTheta;
00417 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
00418 deltaDirection.rotateUz(primaryDirection);
00419
00420 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
00421 {
00422 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
00423
00424 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
00425 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
00426 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
00427 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
00428 finalPx /= finalMomentum;
00429 finalPy /= finalMomentum;
00430 finalPz /= finalMomentum;
00431
00432 G4ThreeVector direction;
00433 direction.set(finalPx,finalPy,finalPz);
00434
00435 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
00436 }
00437
00438 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
00439
00440
00441 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
00442 G4double deexSecEnergy = 0;
00443 for (G4int j=secNumberInit; j < secNumberFinal; j++) {
00444
00445 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
00446
00447 }
00448
00449 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
00450 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
00451
00452 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
00453 fvect->push_back(dp);
00454
00455
00456 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
00457 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
00458 ionizationShell,
00459 theIncomingTrack);
00460 }
00461
00462 }
00463
00464
00465
00466 G4double G4DNABornIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
00467 G4double k, G4int shell)
00468 {
00469 if (particleDefinition == G4Electron::ElectronDefinition())
00470 {
00471 G4double maximumEnergyTransfer=0.;
00472 if ((k+waterStructure.IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
00473 else maximumEnergyTransfer = (k+waterStructure.IonisationEnergy(shell))/2.;
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488 G4double crossSectionMaximum = 0.;
00489
00490 G4double minEnergy = waterStructure.IonisationEnergy(shell);
00491 G4double maxEnergy = maximumEnergyTransfer;
00492 G4int nEnergySteps = 50;
00493
00494 G4double value(minEnergy);
00495 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
00496 G4int step(nEnergySteps);
00497 while (step>0)
00498 {
00499 step--;
00500 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
00501 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
00502 value*=stpEnergy;
00503 }
00504
00505
00506 G4double secondaryElectronKineticEnergy=0.;
00507 do
00508 {
00509 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
00510 } while(G4UniformRand()*crossSectionMaximum >
00511 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
00512
00513 return secondaryElectronKineticEnergy;
00514
00515 }
00516
00517 else if (particleDefinition == G4Proton::ProtonDefinition())
00518 {
00519 G4double maximumKineticEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
00520
00521 G4double crossSectionMaximum = 0.;
00522 for (G4double value = waterStructure.IonisationEnergy(shell);
00523 value<=4.*waterStructure.IonisationEnergy(shell) ;
00524 value+=0.1*eV)
00525 {
00526 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
00527 if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
00528 }
00529
00530 G4double secondaryElectronKineticEnergy = 0.;
00531 do
00532 {
00533 secondaryElectronKineticEnergy = G4UniformRand() * maximumKineticEnergyTransfer;
00534 } while(G4UniformRand()*crossSectionMaximum >=
00535 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
00536
00537 return secondaryElectronKineticEnergy;
00538 }
00539
00540 return 0;
00541 }
00542
00543
00544
00545 void G4DNABornIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
00546 G4double k,
00547 G4double secKinetic,
00548 G4double & cosTheta,
00549 G4double & phi )
00550 {
00551 if (particleDefinition == G4Electron::ElectronDefinition())
00552 {
00553 phi = twopi * G4UniformRand();
00554 if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
00555 else if (secKinetic <= 200.*eV)
00556 {
00557 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
00558 else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
00559 }
00560 else
00561 {
00562 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
00563 cosTheta = std::sqrt(1.-sin2O);
00564 }
00565 }
00566
00567 else if (particleDefinition == G4Proton::ProtonDefinition())
00568 {
00569 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
00570 phi = twopi * G4UniformRand();
00571
00572
00573
00574
00575
00576 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
00577 else cosTheta = (2.*G4UniformRand())-1.;
00578
00579 }
00580 }
00581
00582
00583
00584 double G4DNABornIonisationModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
00585 G4double k,
00586 G4double energyTransfer,
00587 G4int ionizationLevelIndex)
00588 {
00589 G4double sigma = 0.;
00590
00591 if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex))
00592 {
00593 G4double valueT1 = 0;
00594 G4double valueT2 = 0;
00595 G4double valueE21 = 0;
00596 G4double valueE22 = 0;
00597 G4double valueE12 = 0;
00598 G4double valueE11 = 0;
00599
00600 G4double xs11 = 0;
00601 G4double xs12 = 0;
00602 G4double xs21 = 0;
00603 G4double xs22 = 0;
00604
00605 if (particleDefinition == G4Electron::ElectronDefinition())
00606 {
00607
00608
00609 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
00610
00611 std::vector<double>::iterator t1 = t2-1;
00612
00613
00614 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
00615 {
00616 std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
00617 std::vector<double>::iterator e11 = e12-1;
00618
00619 std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
00620 std::vector<double>::iterator e21 = e22-1;
00621
00622 valueT1 =*t1;
00623 valueT2 =*t2;
00624 valueE21 =*e21;
00625 valueE22 =*e22;
00626 valueE12 =*e12;
00627 valueE11 =*e11;
00628
00629 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
00630 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
00631 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
00632 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
00633 }
00634
00635 }
00636
00637 if (particleDefinition == G4Proton::ProtonDefinition())
00638 {
00639
00640 std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
00641 std::vector<double>::iterator t1 = t2-1;
00642
00643 std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
00644 std::vector<double>::iterator e11 = e12-1;
00645
00646 std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
00647 std::vector<double>::iterator e21 = e22-1;
00648
00649 valueT1 =*t1;
00650 valueT2 =*t2;
00651 valueE21 =*e21;
00652 valueE22 =*e22;
00653 valueE12 =*e12;
00654 valueE11 =*e11;
00655
00656 xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
00657 xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
00658 xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
00659 xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
00660
00661 }
00662
00663 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
00664 if (xsProduct != 0.)
00665 {
00666 sigma = QuadInterpolator( valueE11, valueE12,
00667 valueE21, valueE22,
00668 xs11, xs12,
00669 xs21, xs22,
00670 valueT1, valueT2,
00671 k, energyTransfer);
00672 }
00673
00674 }
00675
00676 return sigma;
00677 }
00678
00679
00680
00681 G4double G4DNABornIonisationModel::LogLogInterpolate(G4double e1,
00682 G4double e2,
00683 G4double e,
00684 G4double xs1,
00685 G4double xs2)
00686 {
00687 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
00688 G4double b = std::log10(xs2) - a*std::log10(e2);
00689 G4double sigma = a*std::log10(e) + b;
00690 G4double value = (std::pow(10.,sigma));
00691 return value;
00692 }
00693
00694
00695
00696 G4double G4DNABornIonisationModel::QuadInterpolator(G4double e11, G4double e12,
00697 G4double e21, G4double e22,
00698 G4double xs11, G4double xs12,
00699 G4double xs21, G4double xs22,
00700 G4double t1, G4double t2,
00701 G4double t, G4double e)
00702 {
00703 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
00704 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
00705 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
00706 return value;
00707 }
00708
00709
00710
00711 G4int G4DNABornIonisationModel::RandomSelect(G4double k, const G4String& particle )
00712 {
00713 G4int level = 0;
00714
00715 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00716 pos = tableData.find(particle);
00717
00718 if (pos != tableData.end())
00719 {
00720 G4DNACrossSectionDataSet* table = pos->second;
00721
00722 if (table != 0)
00723 {
00724 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
00725 const size_t n(table->NumberOfComponents());
00726 size_t i(n);
00727 G4double value = 0.;
00728
00729 while (i>0)
00730 {
00731 i--;
00732 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
00733 value += valuesBuffer[i];
00734 }
00735
00736 value *= G4UniformRand();
00737
00738 i = n;
00739
00740 while (i > 0)
00741 {
00742 i--;
00743
00744 if (valuesBuffer[i] > value)
00745 {
00746 delete[] valuesBuffer;
00747 return i;
00748 }
00749 value -= valuesBuffer[i];
00750 }
00751
00752 if (valuesBuffer) delete[] valuesBuffer;
00753
00754 }
00755 }
00756 else
00757 {
00758 G4Exception("G4DNABornIonisationModel::RandomSelect","em0002",
00759 FatalException,"Model not applicable to particle type.");
00760 }
00761
00762 return level;
00763 }
00764