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
00030 #include "G4DNARuddIonisationModel.hh"
00031 #include "G4PhysicalConstants.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "G4UAtomicDeexcitation.hh"
00034 #include "G4LossTableManager.hh"
00035 #include "G4DNAChemistryManager.hh"
00036 #include "G4DNAMolecularMaterial.hh"
00037
00038
00039
00040 using namespace std;
00041
00042
00043
00044 G4DNARuddIonisationModel::G4DNARuddIonisationModel(const G4ParticleDefinition*,
00045 const G4String& nam)
00046 :G4VEmModel(nam),isInitialised(false)
00047 {
00048
00049 fpWaterDensity = 0;
00050
00051 slaterEffectiveCharge[0]=0.;
00052 slaterEffectiveCharge[1]=0.;
00053 slaterEffectiveCharge[2]=0.;
00054 sCoefficient[0]=0.;
00055 sCoefficient[1]=0.;
00056 sCoefficient[2]=0.;
00057
00058 lowEnergyLimitForZ1 = 0 * eV;
00059 lowEnergyLimitForZ2 = 0 * eV;
00060 lowEnergyLimitOfModelForZ1 = 100 * eV;
00061 lowEnergyLimitOfModelForZ2 = 1 * keV;
00062 killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
00063 killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
00064
00065 verboseLevel= 0;
00066
00067
00068
00069
00070
00071
00072
00073 if( verboseLevel>0 )
00074 {
00075 G4cout << "Rudd ionisation model is constructed " << G4endl;
00076 }
00077
00078
00079 SetDeexcitationFlag(true);
00080 fAtomDeexcitation = 0;
00081 fParticleChangeForGamma = 0;
00082 }
00083
00084
00085
00086 G4DNARuddIonisationModel::~G4DNARuddIonisationModel()
00087 {
00088
00089
00090 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00091 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
00092 {
00093 G4DNACrossSectionDataSet* table = pos->second;
00094 delete table;
00095 }
00096
00097
00098
00099
00100
00101
00102 }
00103
00104
00105
00106 void G4DNARuddIonisationModel::Initialise(const G4ParticleDefinition* particle,
00107 const G4DataVector& )
00108 {
00109
00110 if (verboseLevel > 3)
00111 G4cout << "Calling G4DNARuddIonisationModel::Initialise()" << G4endl;
00112
00113
00114
00115 G4String fileProton("dna/sigma_ionisation_p_rudd");
00116 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
00117 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
00118 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
00119 G4String fileHelium("dna/sigma_ionisation_he_rudd");
00120
00121 G4DNAGenericIonsManager *instance;
00122 instance = G4DNAGenericIonsManager::Instance();
00123 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
00124 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
00125 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
00126 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
00127 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
00128
00129 G4String proton;
00130 G4String hydrogen;
00131 G4String alphaPlusPlus;
00132 G4String alphaPlus;
00133 G4String helium;
00134
00135 G4double scaleFactor = 1 * m*m;
00136
00137
00138
00139
00140
00141 proton = protonDef->GetParticleName();
00142 tableFile[proton] = fileProton;
00143
00144 lowEnergyLimit[proton] = lowEnergyLimitForZ1;
00145 highEnergyLimit[proton] = 500. * keV;
00146
00147
00148
00149 G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00150 eV,
00151 scaleFactor );
00152 tableProton->LoadData(fileProton);
00153 tableData[proton] = tableProton;
00154
00155
00156
00157 hydrogen = hydrogenDef->GetParticleName();
00158 tableFile[hydrogen] = fileHydrogen;
00159
00160 lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
00161 highEnergyLimit[hydrogen] = 100. * MeV;
00162
00163
00164
00165 G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00166 eV,
00167 scaleFactor );
00168 tableHydrogen->LoadData(fileHydrogen);
00169
00170 tableData[hydrogen] = tableHydrogen;
00171
00172
00173
00174 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
00175 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
00176
00177 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
00178 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
00179
00180
00181
00182 G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00183 eV,
00184 scaleFactor );
00185 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
00186
00187 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
00188
00189
00190
00191 alphaPlus = alphaPlusDef->GetParticleName();
00192 tableFile[alphaPlus] = fileAlphaPlus;
00193
00194 lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
00195 highEnergyLimit[alphaPlus] = 400. * MeV;
00196
00197
00198
00199 G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00200 eV,
00201 scaleFactor );
00202 tableAlphaPlus->LoadData(fileAlphaPlus);
00203 tableData[alphaPlus] = tableAlphaPlus;
00204
00205
00206
00207 helium = heliumDef->GetParticleName();
00208 tableFile[helium] = fileHelium;
00209
00210 lowEnergyLimit[helium] = lowEnergyLimitForZ2;
00211 highEnergyLimit[helium] = 400. * MeV;
00212
00213
00214
00215 G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
00216 eV,
00217 scaleFactor );
00218 tableHelium->LoadData(fileHelium);
00219 tableData[helium] = tableHelium;
00220
00221
00222
00223 if (particle==protonDef)
00224 {
00225 SetLowEnergyLimit(lowEnergyLimit[proton]);
00226 SetHighEnergyLimit(highEnergyLimit[proton]);
00227 }
00228
00229 if (particle==hydrogenDef)
00230 {
00231 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
00232 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
00233 }
00234
00235 if (particle==heliumDef)
00236 {
00237 SetLowEnergyLimit(lowEnergyLimit[helium]);
00238 SetHighEnergyLimit(highEnergyLimit[helium]);
00239 }
00240
00241 if (particle==alphaPlusDef)
00242 {
00243 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
00244 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
00245 }
00246
00247 if (particle==alphaPlusPlusDef)
00248 {
00249 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
00250 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
00251 }
00252
00253 if( verboseLevel>0 )
00254 {
00255 G4cout << "Rudd ionisation model is initialized " << G4endl
00256 << "Energy range: "
00257 << LowEnergyLimit() / eV << " eV - "
00258 << HighEnergyLimit() / keV << " keV for "
00259 << particle->GetParticleName()
00260 << G4endl;
00261 }
00262
00263
00264 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00265
00266
00267
00268 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00269
00270 if (isInitialised) { return; }
00271 fParticleChangeForGamma = GetParticleChangeForGamma();
00272 isInitialised = true;
00273
00274 }
00275
00276
00277
00278 G4double G4DNARuddIonisationModel::CrossSectionPerVolume(const G4Material* material,
00279 const G4ParticleDefinition* particleDefinition,
00280 G4double k,
00281 G4double,
00282 G4double)
00283 {
00284 if (verboseLevel > 3)
00285 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationModel" << G4endl;
00286
00287
00288
00289 G4DNAGenericIonsManager *instance;
00290 instance = G4DNAGenericIonsManager::Instance();
00291
00292 if (
00293 particleDefinition != G4Proton::ProtonDefinition()
00294 &&
00295 particleDefinition != instance->GetIon("hydrogen")
00296 &&
00297 particleDefinition != instance->GetIon("alpha++")
00298 &&
00299 particleDefinition != instance->GetIon("alpha+")
00300 &&
00301 particleDefinition != instance->GetIon("helium")
00302 )
00303
00304 return 0;
00305
00306 G4double lowLim = 0;
00307
00308 if ( particleDefinition == G4Proton::ProtonDefinition()
00309 || particleDefinition == instance->GetIon("hydrogen")
00310 )
00311
00312 lowLim = lowEnergyLimitOfModelForZ1;
00313
00314 if ( particleDefinition == instance->GetIon("alpha++")
00315 || particleDefinition == instance->GetIon("alpha+")
00316 || particleDefinition == instance->GetIon("helium")
00317 )
00318
00319 lowLim = lowEnergyLimitOfModelForZ2;
00320
00321 G4double highLim = 0;
00322 G4double sigma=0;
00323
00324 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
00325
00326 if(waterDensity!= 0.0)
00327
00328 {
00329 const G4String& particleName = particleDefinition->GetParticleName();
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00343 pos2 = highEnergyLimit.find(particleName);
00344
00345 if (pos2 != highEnergyLimit.end())
00346 {
00347 highLim = pos2->second;
00348 }
00349
00350 if (k <= highLim)
00351 {
00352
00353
00354 if (k < lowLim) k = lowLim;
00355
00356
00357
00358 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00359 pos = tableData.find(particleName);
00360
00361 if (pos != tableData.end())
00362 {
00363 G4DNACrossSectionDataSet* table = pos->second;
00364 if (table != 0)
00365 {
00366 sigma = table->FindValue(k);
00367 }
00368 }
00369 else
00370 {
00371 G4Exception("G4DNARuddIonisationModel::CrossSectionPerVolume","em0002",
00372 FatalException,"Model not applicable to particle type.");
00373 }
00374
00375 }
00376
00377 if (verboseLevel > 2)
00378 {
00379 G4cout << "__________________________________" << G4endl;
00380 G4cout << "°°° G4DNARuddIonisationModel - XS INFO START" << G4endl;
00381 G4cout << "°°° Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
00382 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
00383 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
00384
00385 G4cout << "°°° G4DNARuddIonisationModel - XS INFO END" << G4endl;
00386 }
00387
00388 }
00389 else
00390 {
00391 if (verboseLevel > 2)
00392 {
00393 G4cout << "Warning : RuddIonisationModel: WATER DENSITY IS NULL" << G4endl;
00394 }
00395 }
00396
00397 return sigma*waterDensity;
00398
00399
00400 }
00401
00402
00403
00404 void G4DNARuddIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00405 const G4MaterialCutsCouple* ,
00406 const G4DynamicParticle* particle,
00407 G4double,
00408 G4double)
00409 {
00410 if (verboseLevel > 3)
00411 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationModel" << G4endl;
00412
00413 G4double lowLim = 0;
00414 G4double highLim = 0;
00415
00416 G4DNAGenericIonsManager *instance;
00417 instance = G4DNAGenericIonsManager::Instance();
00418
00419 if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
00420 || particle->GetDefinition() == instance->GetIon("hydrogen")
00421 )
00422
00423 lowLim = killBelowEnergyForZ1;
00424
00425 if ( particle->GetDefinition() == instance->GetIon("alpha++")
00426 || particle->GetDefinition() == instance->GetIon("alpha+")
00427 || particle->GetDefinition() == instance->GetIon("helium")
00428 )
00429
00430 lowLim = killBelowEnergyForZ2;
00431
00432 G4double k = particle->GetKineticEnergy();
00433
00434 const G4String& particleName = particle->GetDefinition()->GetParticleName();
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00448 pos2 = highEnergyLimit.find(particleName);
00449
00450 if (pos2 != highEnergyLimit.end())
00451 {
00452 highLim = pos2->second;
00453 }
00454
00455 if (k >= lowLim && k <= highLim)
00456 {
00457 G4ParticleDefinition* definition = particle->GetDefinition();
00458 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
00459
00460
00461
00462
00463
00464
00465
00466 G4int ionizationShell = RandomSelect(k,particleName);
00467
00468
00469
00470
00471
00472 G4int secNumberInit = 0;
00473 G4int secNumberFinal = 0;
00474
00475 G4double bindingEnergy = 0;
00476 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
00477
00478 if(fAtomDeexcitation) {
00479 G4int Z = 8;
00480 G4AtomicShellEnumerator as = fKShell;
00481
00482 if (ionizationShell <5 && ionizationShell >1)
00483 {
00484 as = G4AtomicShellEnumerator(4-ionizationShell);
00485 }
00486 else if (ionizationShell <2)
00487 {
00488 as = G4AtomicShellEnumerator(3);
00489 }
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00503 secNumberInit = fvect->size();
00504 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
00505 secNumberFinal = fvect->size();
00506 }
00507
00508
00509 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
00510
00511 G4double cosTheta = 0.;
00512 G4double phi = 0.;
00513 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi);
00514
00515 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
00516 G4double dirX = sinTheta*std::cos(phi);
00517 G4double dirY = sinTheta*std::sin(phi);
00518 G4double dirZ = cosTheta;
00519 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
00520 deltaDirection.rotateUz(primaryDirection);
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
00540
00541 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
00542 G4double deexSecEnergy = 0;
00543 for (G4int j=secNumberInit; j < secNumberFinal; j++) {
00544
00545 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
00546
00547 }
00548
00549 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
00550 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
00551
00552
00553
00554
00555
00556
00557
00558 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
00559 fvect->push_back(dp);
00560
00561 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
00562 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
00563 ionizationShell,
00564 theIncomingTrack);
00565 }
00566
00567
00568
00569 if (k < lowLim)
00570 {
00571 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
00572 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
00573 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
00574 }
00575
00576 }
00577
00578
00579
00580 G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
00581 G4double k,
00582 G4int shell)
00583 {
00584 G4double maximumKineticEnergyTransfer = 0.;
00585
00586 G4DNAGenericIonsManager *instance;
00587 instance = G4DNAGenericIonsManager::Instance();
00588
00589 if (particleDefinition == G4Proton::ProtonDefinition()
00590 || particleDefinition == instance->GetIon("hydrogen"))
00591 {
00592 maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
00593 }
00594
00595 else if (particleDefinition == instance->GetIon("helium")
00596 || particleDefinition == instance->GetIon("alpha+")
00597 || particleDefinition == instance->GetIon("alpha++"))
00598 {
00599 maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
00600 }
00601
00602 G4double crossSectionMaximum = 0.;
00603
00604 for(G4double value=waterStructure.IonisationEnergy(shell); value<=5.*waterStructure.IonisationEnergy(shell) && k>=value ; value+=0.1*eV)
00605 {
00606 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
00607 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
00608 }
00609
00610
00611 G4double secElecKinetic = 0.;
00612
00613 do
00614 {
00615 secElecKinetic = G4UniformRand() * maximumKineticEnergyTransfer;
00616 } while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
00617 k,
00618 secElecKinetic+waterStructure.IonisationEnergy(shell),
00619 shell));
00620
00621 return(secElecKinetic);
00622 }
00623
00624
00625
00626
00627 void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
00628 G4double k,
00629 G4double secKinetic,
00630 G4double & cosTheta,
00631 G4double & phi )
00632 {
00633 G4DNAGenericIonsManager *instance;
00634 instance = G4DNAGenericIonsManager::Instance();
00635
00636 G4double maxSecKinetic = 0.;
00637
00638 if (particleDefinition == G4Proton::ProtonDefinition()
00639 || particleDefinition == instance->GetIon("hydrogen"))
00640 {
00641 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
00642 }
00643
00644 else if (particleDefinition == instance->GetIon("helium")
00645 || particleDefinition == instance->GetIon("alpha+")
00646 || particleDefinition == instance->GetIon("alpha++"))
00647 {
00648 maxSecKinetic = 4.* (0.511 / 3728) * k;
00649 }
00650
00651 phi = twopi * G4UniformRand();
00652
00653
00654
00655
00656
00657 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
00658 else cosTheta = (2.*G4UniformRand())-1.;
00659
00660 }
00661
00662
00663
00664
00665 G4double G4DNARuddIonisationModel::DifferentialCrossSection(G4ParticleDefinition* particleDefinition,
00666 G4double k,
00667 G4double energyTransfer,
00668 G4int ionizationLevelIndex)
00669 {
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686 const G4int j=ionizationLevelIndex;
00687
00688 G4double A1 ;
00689 G4double B1 ;
00690 G4double C1 ;
00691 G4double D1 ;
00692 G4double E1 ;
00693 G4double A2 ;
00694 G4double B2 ;
00695 G4double C2 ;
00696 G4double D2 ;
00697 G4double alphaConst ;
00698
00699
00700
00701 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
00702
00703 if (j == 4)
00704 {
00705
00706 A1 = 1.25;
00707 B1 = 0.5;
00708 C1 = 1.00;
00709 D1 = 1.00;
00710 E1 = 3.00;
00711 A2 = 1.10;
00712 B2 = 1.30;
00713 C2 = 1.00;
00714 D2 = 0.00;
00715 alphaConst = 0.66;
00716 }
00717 else
00718 {
00719
00720 A1 = 1.02;
00721 B1 = 82.0;
00722 C1 = 0.45;
00723 D1 = -0.80;
00724 E1 = 0.38;
00725 A2 = 1.07;
00726
00727 B2 = 11.6;
00728
00729 C2 = 0.60;
00730 D2 = 0.04;
00731 alphaConst = 0.64;
00732 }
00733
00734 const G4double n = 2.;
00735 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
00736
00737 G4DNAGenericIonsManager* instance;
00738 instance = G4DNAGenericIonsManager::Instance();
00739
00740 G4double wBig = (energyTransfer - waterStructure.IonisationEnergy(ionizationLevelIndex));
00741 if (wBig<0) return 0.;
00742
00743 G4double w = wBig / Bj[ionizationLevelIndex];
00744
00745 if (j==4) w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
00746
00747 G4double Ry = 13.6*eV;
00748
00749 G4double tau = 0.;
00750
00751 G4bool isProtonOrHydrogen = false;
00752 G4bool isHelium = false;
00753
00754 if (particleDefinition == G4Proton::ProtonDefinition()
00755 || particleDefinition == instance->GetIon("hydrogen"))
00756 {
00757 isProtonOrHydrogen = true;
00758 tau = (electron_mass_c2/proton_mass_c2) * k ;
00759 }
00760
00761 else if ( particleDefinition == instance->GetIon("helium")
00762 || particleDefinition == instance->GetIon("alpha+")
00763 || particleDefinition == instance->GetIon("alpha++"))
00764 {
00765 isHelium = true;
00766 tau = (0.511/3728.) * k ;
00767 }
00768
00769 G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/Bj[ionizationLevelIndex]),2);
00770 if (j==4) S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
00771
00772 G4double v2 = tau / Bj[ionizationLevelIndex];
00773 if (j==4) v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
00774
00775 G4double v = std::sqrt(v2);
00776 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj[ionizationLevelIndex]));
00777 if (j==4) wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.IonisationEnergy(ionizationLevelIndex)));
00778
00779 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
00780 G4double L2 = C2*std::pow(v,(D2));
00781 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
00782 G4double H2 = (A2/v2) + (B2/(v2*v2));
00783
00784 G4double F1 = L1+H1;
00785 G4double F2 = (L2*H2)/(L2+H2);
00786
00787 G4double sigma = CorrectionFactor(particleDefinition, k)
00788 * Gj[j] * (S/Bj[ionizationLevelIndex])
00789 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
00790
00791 if (j==4) sigma = CorrectionFactor(particleDefinition, k)
00792 * Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
00793 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
00794
00795 if ( (particleDefinition == instance->GetIon("hydrogen")) && (ionizationLevelIndex==4))
00796 {
00797
00798 sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
00799 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
00800 }
00801
00802
00803
00804 if(isProtonOrHydrogen)
00805 {
00806 return(sigma);
00807 }
00808
00809 if (particleDefinition == instance->GetIon("alpha++") )
00810 {
00811 slaterEffectiveCharge[0]=0.;
00812 slaterEffectiveCharge[1]=0.;
00813 slaterEffectiveCharge[2]=0.;
00814 sCoefficient[0]=0.;
00815 sCoefficient[1]=0.;
00816 sCoefficient[2]=0.;
00817 }
00818
00819 else if (particleDefinition == instance->GetIon("alpha+") )
00820 {
00821 slaterEffectiveCharge[0]=2.0;
00822
00823 slaterEffectiveCharge[1]=2.0;
00824 slaterEffectiveCharge[2]=2.0;
00825
00826 sCoefficient[0]=0.7;
00827 sCoefficient[1]=0.15;
00828 sCoefficient[2]=0.15;
00829 }
00830
00831 else if (particleDefinition == instance->GetIon("helium") )
00832 {
00833 slaterEffectiveCharge[0]=1.7;
00834 slaterEffectiveCharge[1]=1.15;
00835 slaterEffectiveCharge[2]=1.15;
00836 sCoefficient[0]=0.5;
00837 sCoefficient[1]=0.25;
00838 sCoefficient[2]=0.25;
00839 }
00840
00841
00842
00843
00844
00845 if(isHelium)
00846 {
00847 sigma = Gj[j] * (S/Bj[ionizationLevelIndex]) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
00848
00849 if (j==4) sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
00850 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
00851
00852 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
00853
00854 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
00855 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
00856 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
00857
00858 return zEff * zEff * sigma ;
00859 }
00860
00861 return 0;
00862 }
00863
00864
00865
00866 G4double G4DNARuddIonisationModel::S_1s(G4double t,
00867 G4double energyTransferred,
00868 G4double slaterEffectiveChg,
00869 G4double shellNumber)
00870 {
00871
00872
00873
00874 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
00875 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
00876
00877 return value;
00878 }
00879
00880
00881
00882 G4double G4DNARuddIonisationModel::S_2s(G4double t,
00883 G4double energyTransferred,
00884 G4double slaterEffectiveChg,
00885 G4double shellNumber)
00886 {
00887
00888
00889
00890 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
00891 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
00892
00893 return value;
00894
00895 }
00896
00897
00898
00899 G4double G4DNARuddIonisationModel::S_2p(G4double t,
00900 G4double energyTransferred,
00901 G4double slaterEffectiveChg,
00902 G4double shellNumber)
00903 {
00904
00905
00906
00907 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
00908 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
00909
00910 return value;
00911 }
00912
00913
00914
00915 G4double G4DNARuddIonisationModel::R(G4double t,
00916 G4double energyTransferred,
00917 G4double slaterEffectiveChg,
00918 G4double shellNumber)
00919 {
00920
00921
00922
00923 G4double tElectron = 0.511/3728. * t;
00924
00925 G4double H = 2.*13.60569172 * eV;
00926 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
00927
00928 return value;
00929 }
00930
00931
00932
00933 G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k)
00934 {
00935 G4DNAGenericIonsManager *instance;
00936 instance = G4DNAGenericIonsManager::Instance();
00937
00938 if (particleDefinition == G4Proton::Proton())
00939 {
00940 return(1.);
00941 }
00942 else
00943 if (particleDefinition == instance->GetIon("hydrogen"))
00944 {
00945 G4double value = (std::log10(k/eV)-4.2)/0.5;
00946
00947 return((0.6/(1+std::exp(value))) + 0.9);
00948 }
00949 else
00950 {
00951 return(1.);
00952 }
00953 }
00954
00955
00956
00957 G4int G4DNARuddIonisationModel::RandomSelect(G4double k, const G4String& particle )
00958 {
00959
00960
00961
00962
00963
00964 G4int level = 0;
00965
00966
00967
00968 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00969 pos = tableData.find(particle);
00970
00971 if (pos != tableData.end())
00972 {
00973 G4DNACrossSectionDataSet* table = pos->second;
00974
00975 if (table != 0)
00976 {
00977 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
00978
00979 const size_t n(table->NumberOfComponents());
00980 size_t i(n);
00981 G4double value = 0.;
00982
00983 while (i>0)
00984 {
00985 i--;
00986 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
00987 value += valuesBuffer[i];
00988 }
00989
00990 value *= G4UniformRand();
00991
00992 i = n;
00993
00994 while (i > 0)
00995 {
00996 i--;
00997
00998
00999 if (valuesBuffer[i] > value)
01000 {
01001 delete[] valuesBuffer;
01002 return i;
01003 }
01004 value -= valuesBuffer[i];
01005 }
01006
01007 if (valuesBuffer) delete[] valuesBuffer;
01008
01009 }
01010 }
01011 else
01012 {
01013 G4Exception("G4DNARuddIonisationModel::RandomSelect","em0002",
01014 FatalException,"Model not applicable to particle type.");
01015 }
01016
01017 return level;
01018 }
01019
01020
01021
01022 G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track )
01023 {
01024 G4double sigma = 0.;
01025
01026 const G4DynamicParticle* particle = track.GetDynamicParticle();
01027 G4double k = particle->GetKineticEnergy();
01028
01029 G4double lowLim = 0;
01030 G4double highLim = 0;
01031
01032 const G4String& particleName = particle->GetDefinition()->GetParticleName();
01033
01034 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
01035 pos1 = lowEnergyLimit.find(particleName);
01036
01037 if (pos1 != lowEnergyLimit.end())
01038 {
01039 lowLim = pos1->second;
01040 }
01041
01042 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
01043 pos2 = highEnergyLimit.find(particleName);
01044
01045 if (pos2 != highEnergyLimit.end())
01046 {
01047 highLim = pos2->second;
01048 }
01049
01050 if (k >= lowLim && k <= highLim)
01051 {
01052 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
01053 pos = tableData.find(particleName);
01054
01055 if (pos != tableData.end())
01056 {
01057 G4DNACrossSectionDataSet* table = pos->second;
01058 if (table != 0)
01059 {
01060 sigma = table->FindValue(k);
01061 }
01062 }
01063 else
01064 {
01065 G4Exception("G4DNARuddIonisationModel::PartialCrossSection","em0002",
01066 FatalException,"Model not applicable to particle type.");
01067 }
01068 }
01069
01070 return sigma;
01071 }
01072
01073
01074
01075 G4double G4DNARuddIonisationModel::Sum(G4double , const G4String& )
01076 {
01077 return 0;
01078 }
01079