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