Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4DNARuddIonisationExtendedModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// Modified by Z. Francis, S. Incerti to handle HZE
28// && inverse rudd function sampling 26-10-2010
29
30#include "G4DNARuddIonisationExtendedModel.hh"
31#include "G4PhysicalConstants.hh"
32#include "G4SystemOfUnits.hh"
33#include "G4UAtomicDeexcitation.hh"
34#include "G4LossTableManager.hh"
35#include "G4DNAChemistryManager.hh"
36#include "G4DNAMolecularMaterial.hh"
37
38#include "G4IonTable.hh"
39#include "G4DNARuddAngle.hh"
40#include "G4DeltaAngle.hh"
41#include "G4Exp.hh"
42#include "G4Log.hh"
43#include "G4Pow.hh"
44#include "G4Alpha.hh"
45
46static G4Pow * gpow = G4Pow::GetInstance();
47
48
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50
51using namespace std;
52
53//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54
55G4DNARuddIonisationExtendedModel::G4DNARuddIonisationExtendedModel(const G4ParticleDefinition*,
56 const G4String& nam)
57:G4VEmModel(nam),isInitialised(false)
58{
59 slaterEffectiveCharge[0]=0.;
60 slaterEffectiveCharge[1]=0.;
61 slaterEffectiveCharge[2]=0.;
62 sCoefficient[0]=0.;
63 sCoefficient[1]=0.;
64 sCoefficient[2]=0.;
65
66 lowEnergyLimitForA[1] = 0 * eV;
67 lowEnergyLimitForA[2] = 0 * eV;
68 lowEnergyLimitForA[3] = 0 * eV;
69 lowEnergyLimitOfModelForA[1] = 100 * eV;
70 lowEnergyLimitOfModelForA[4] = 1 * keV;
71 lowEnergyLimitOfModelForA[5] = 0.5 * MeV; // For A = 3 or above, limit is MeV/uma
72 killBelowEnergyForA[1] = lowEnergyLimitOfModelForA[1];
73 killBelowEnergyForA[4] = lowEnergyLimitOfModelForA[4];
74 killBelowEnergyForA[5] = lowEnergyLimitOfModelForA[5];
75
76 verboseLevel= 0;
77 // Verbosity scale:
78 // 0 = nothing
79 // 1 = warning for energy non-conservation
80 // 2 = details of energy budget
81 // 3 = calculation of cross sections, file openings, sampling of atoms
82 // 4 = entering in methods
83
84 if( verboseLevel>0 )
85 {
86 G4cout << "Rudd ionisation model is constructed " << G4endl;
87 }
88
89 // Define default angular generator
90 SetAngularDistribution(new G4DNARuddAngle());
91
92 // Mark this model as "applicable" for atomic deexcitation
93 SetDeexcitationFlag(true);
94
95 // Selection of stationary mode
96
97 statCode = false;
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
102G4DNARuddIonisationExtendedModel::~G4DNARuddIonisationExtendedModel()
103{
104 if(isIon) {
105 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
106 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
107 {
108 G4DNACrossSectionDataSet* table = pos->second;
109 delete table;
110 }
111 } else {
112 delete mainTable;
113 }
114}
115
116//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117
118void G4DNARuddIonisationExtendedModel::Initialise(const G4ParticleDefinition* particle,
119 const G4DataVector& /*cuts*/)
120{
121 if (isInitialised) { return; }
122 if (verboseLevel > 3)
123 G4cout << "Calling G4DNARuddIonisationExtendedModel::Initialise()" << G4endl;
124
125 // Energy limits
126 G4String fileProton("dna/sigma_ionisation_p_rudd");
127 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
128 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
129 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
130 G4String fileHelium("dna/sigma_ionisation_he_rudd");
131 G4String fileCarbon("dna/sigma_ionisation_c_rudd");
132 G4String fileNitrogen("dna/sigma_ionisation_n_rudd");
133 G4String fileOxygen("dna/sigma_ionisation_o_rudd");
134 G4String fileSilicon("dna/sigma_ionisation_si_rudd");
135 G4String fileIron("dna/sigma_ionisation_fe_rudd");
136
137 G4String pname = particle->GetParticleName();
138
139 G4DNAGenericIonsManager *instance;
140 instance = G4DNAGenericIonsManager::Instance();
141 protonDef = G4Proton::ProtonDefinition();
142 hydrogenDef = instance->GetIon("hydrogen");
143 alphaPlusPlusDef = G4Alpha::Alpha();
144 alphaPlusDef = instance->GetIon("alpha+");
145 heliumDef = instance->GetIon("helium");
146
147 carbonDef = instance->GetIon("carbon");
148 nitrogenDef = instance->GetIon("nitrogen");
149 oxygenDef = instance->GetIon("oxygen");
150 siliconDef = instance->GetIon("silicon");
151 ironDef = instance->GetIon("iron");
152
153 G4String carbon;
154 G4String nitrogen;
155 G4String oxygen;
156 G4String silicon;
157 G4String iron;
158
159 G4double scaleFactor = 1 * m*m;
160 massC12 = carbonDef->GetPDGMass();
161
162 // LIMITS AND DATA
163
164 // **********************************************************************************************
165
166 if(pname == "proton") {
167 localMinEnergy = lowEnergyLimitForA[1];
168
169 // Cross section
170 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
171 mainTable->LoadData(fileProton);
172
173 // **********************************************************************************************
174
175 } else if(pname == "hydrogen") {
176
177 localMinEnergy = lowEnergyLimitForA[1];
178
179 // Cross section
180 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
181 mainTable->LoadData(fileHydrogen);
182
183 // **********************************************************************************************
184
185 } else if(pname == "alpha") {
186
187 localMinEnergy = lowEnergyLimitForA[4];
188
189 // Cross section
190 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
191 mainTable->LoadData(fileAlphaPlusPlus);
192
193 // **********************************************************************************************
194
195 } else if(pname == "alpha+") {
196
197 localMinEnergy = lowEnergyLimitForA[4];
198
199 // Cross section
200 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
201 mainTable->LoadData(fileAlphaPlus);
202
203 // **********************************************************************************************
204
205 } else if(pname == "helium") {
206
207 localMinEnergy = lowEnergyLimitForA[4];
208
209 // Cross section
210 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
211 mainTable->LoadData(fileHelium);
212
213 // **********************************************************************************************
214
215 } else if(pname == "GenericIon") {
216
217 isIon = true;
218 carbon = carbonDef->GetParticleName();
219 localMinEnergy = lowEnergyLimitForA[5]*massC12/proton_mass_c2;
220
221 // Cross section
222 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
223 mainTable->LoadData(fileCarbon);
224
225 tableData[carbon] = mainTable;
226
227 // **********************************************************************************************
228
229 oxygen = oxygenDef->GetParticleName();
230 tableFile[oxygen] = fileOxygen;
231
232 // Cross section
233 G4DNACrossSectionDataSet* tableOxygen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
234 eV, scaleFactor );
235 tableOxygen->LoadData(fileOxygen);
236 tableData[oxygen] = tableOxygen;
237
238 // **********************************************************************************************
239
240 nitrogen = nitrogenDef->GetParticleName();
241 tableFile[nitrogen] = fileNitrogen;
242
243 // Cross section
244 G4DNACrossSectionDataSet* tableNitrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
245 eV, scaleFactor );
246 tableNitrogen->LoadData(fileNitrogen);
247 tableData[nitrogen] = tableNitrogen;
248
249 // **********************************************************************************************
250
251 silicon = siliconDef->GetParticleName();
252 tableFile[silicon] = fileSilicon;
253
254 // Cross section
255 G4DNACrossSectionDataSet* tableSilicon = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
256 eV, scaleFactor );
257 tableSilicon->LoadData(fileSilicon);
258 tableData[silicon] = tableSilicon;
259
260 // **********************************************************************************************
261
262 iron = ironDef->GetParticleName();
263 tableFile[iron] = fileIron;
264
265 // Cross section
266
267 G4DNACrossSectionDataSet* tableIron = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
268 eV, scaleFactor );
269 tableIron->LoadData(fileIron);
270 tableData[iron] = tableIron;
271 }
272 // **********************************************************************************************
273
274 if( verboseLevel>0 )
275 {
276 G4cout << "Rudd ionisation model is initialized " << G4endl
277 << "Energy range for model: "
278 << LowEnergyLimit() / keV << " keV - "
279 << HighEnergyLimit() / keV << " keV for "
280 << pname
281 << " internal low energy limit E(keV)=" << localMinEnergy / keV
282 << G4endl;
283 }
284
285 // Initialize water density pointer
286 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
287
288 //
289 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
290
291 fParticleChangeForGamma = GetParticleChangeForGamma();
292 isInitialised = true;
293}
294
295//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296
297G4double G4DNARuddIonisationExtendedModel::CrossSectionPerVolume(const G4Material* material,
298 const G4ParticleDefinition* partDef,
299 G4double k,
300 G4double,
301 G4double)
302{
303 if (verboseLevel > 3)
304 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" << G4endl;
305
306 currParticle = GetDNAIonParticleDefinition(partDef);
307 currentScaledEnergy = k;
308 G4double e = k;
309 G4double q2 = 1.0;
310 currentTable = mainTable;
311
312 if (isIon){
313 if (currParticle == nullptr) {//not DNA particle
314 currentScaledEnergy *= massC12/partDef->GetPDGMass();
315 G4double q = partDef->GetPDGCharge()/(eplus*6);
316 q2 *= q*q;
317 e = currentScaledEnergy;
318 currParticle = carbonDef;
319 }
320 G4String pname = currParticle->GetParticleName();
321 auto goodTable = tableData.find(pname);
322 currentTable = goodTable->second;
323 }
324 // below low the ion should be stopped
325 if (currentScaledEnergy < localMinEnergy) { return DBL_MAX; }
326
327 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
328 G4double sigma = 0.0;
329 if (nullptr != currentTable) {
330 sigma = currentTable->FindValue(e)*q2;
331 } else {
332 G4cout << "G4DNARuddIonisationExtendedModel - no data table for "
333 << partDef->GetParticleName() << G4endl;
334 G4Exception("G4DNARuddIonisationExtendedModel::CrossSectionPerVolume(...)","em0002",
335 FatalException,"Data table is not available for the model.");
336 }
337 if (verboseLevel > 2)
338 {
339 G4cout << "__________________________________" << G4endl;
340 G4cout << "G4DNARuddIonisationExtendedModel - XS INFO START" << G4endl;
341 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << partDef->GetParticleName() << G4endl;
342 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
343 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
344 G4cout << "G4DNARuddIonisationExtendedModel - XS INFO END" << G4endl;
345 }
346 return sigma*waterDensity;
347}
348
349//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
350
351void G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
352 const G4MaterialCutsCouple* couple,
353 const G4DynamicParticle* particle,
354 G4double,
355 G4double)
356{
357 if (verboseLevel > 3)
358 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" << G4endl;
359
360 // stop ion with energy below low energy limit
361 G4double k = particle->GetKineticEnergy();
362 if (currentScaledEnergy < localMinEnergy) {
363 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
364 fParticleChangeForGamma->ProposeTrackStatus(fStopButAlive);
365 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
366 return;
367 }
368
369 // sampling of final state
370 G4ParticleDefinition* definition = particle->GetDefinition();
371 const G4ThreeVector& primaryDirection = particle->GetMomentumDirection();
372
373 G4int ionizationShell = RandomSelect(currentScaledEnergy);
374
375 // sample deexcitation
376 // here we assume that H_{2}O electronic levels are the same as Oxygen.
377 // this can be considered true with a rough 10% error in energy on K-shell,
378
379 G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
380
381 //SI: additional protection if tcs interpolation method is modified
382 if (k < bindingEnergy) return;
383 //
384
385 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
386
387 // is ionisation possible?
388 G4double scatteredEnergy = k - bindingEnergy - secondaryKinetic;
389 if(scatteredEnergy < 0.0) { return; }
390
391 G4int Z = 8;
392
393 G4ThreeVector deltaDirection =
394 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
395 Z, ionizationShell,
396 couple->GetMaterial());
397
398 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
399 fvect->push_back(dp);
400
401 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
402
403 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
404 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
405
406 // SI: only atomic deexcitation from K shell is considered
407 if(fAtomDeexcitation != nullptr && ionizationShell == 4)
408 {
409 const G4AtomicShell* shell
410 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
411 secNumberInit = fvect->size();
412 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
413 secNumberFinal = fvect->size();
414
415 if(secNumberFinal > secNumberInit)
416 {
417 for (size_t i=secNumberInit; i<secNumberFinal; ++i)
418 {
419 //Check if there is enough residual energy
420 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
421 {
422 //Ok, this is a valid secondary: keep it
423 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
424 }
425 else
426 {
427 //Invalid secondary: not enough energy to create it!
428 //Keep its energy in the local deposit
429 delete (*fvect)[i];
430 (*fvect)[i]=0;
431 }
432 }
433 }
434 }
435
436 //bindingEnergy has been decreased
437 //by the amount of energy taken away by deexc. products
438 if (!statCode)
439 {
440 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
441 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
442 }
443 else
444 {
445 fParticleChangeForGamma->SetProposedKineticEnergy(k);
446 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
447 }
448
449 // TEST //////////////////////////
450 // if (secondaryKinetic<0) abort();
451 // if (scatteredEnergy<0) abort();
452 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
453 // if (k-scatteredEnergy<0) abort();
454 /////////////////////////////////
455
456 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
457 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
458 ionizationShell,
459 theIncomingTrack);
460 // SI - not useful since low energy of model is 0 eV
461}
462
463//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
464
465G4double G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
466 G4double k,
467 G4int shell)
468{
469 //-- Fast sampling method -----
470 G4double proposed_energy;
471 G4double random1;
472 G4double value_sampling;
473 G4double max1;
474
475 do
476 {
477 proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell); // Proposed energy by inverse function sampling
478
479 max1=0.;
480
481 for(G4double en=0.; en<20.; en+=1.) if(RejectionFunction(particleDefinition, k, en, shell) > max1)
482 max1=RejectionFunction(particleDefinition, k, en, shell);
483
484 random1 = G4UniformRand()*max1;
485
486 value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell);
487
488 } while(random1 > value_sampling);
489
490 return(proposed_energy);
491}
492
493//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
494
495G4double G4DNARuddIonisationExtendedModel::RejectionFunction(G4ParticleDefinition* particleDefinition,
496 G4double k,
497 G4double proposed_ws,
498 G4int ionizationLevelIndex)
499{
500 const G4int j=ionizationLevelIndex;
501 G4double Bj_energy, alphaConst;
502 G4double Ry = 13.6*eV;
503 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
504
505 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
506
507 // Following values provided by M. Dingfelder (priv. comm)
508 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
509
510 if (j == 4)
511 {
512 alphaConst = 0.66;
513 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
514 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
515 //---
516 }
517 else
518 {
519 alphaConst = 0.64;
520 Bj_energy = Bj[ionizationLevelIndex];
521 }
522
523 G4double energyTransfer = proposed_ws + Bj_energy;
524 proposed_ws/=Bj_energy;
525
526 G4double tau = 0.;
527 G4double A_ion = 0.;
528 tau = (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
529 A_ion = particleDefinition->GetAtomicMass();
530
531 G4double v2;
532 G4double beta2;
533
534 if((tau/MeV)<5.447761194e-2)
535 {
536 v2 = tau / Bj_energy;
537 beta2 = 2.*tau / electron_mass_c2;
538 }
539 // Relativistic
540 else
541 {
542 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
543 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
544 }
545
546 G4double v = std::sqrt(v2);
547 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
548 G4double rejection_term = 1.+G4Exp(alphaConst*(proposed_ws - wc) / v);
549 rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
550 //* (S/Bj_energy) ; Not needed anymore
551
552 G4bool isHelium = false;
553
554 if ( particleDefinition == protonDef
555 || particleDefinition == hydrogenDef
556 )
557 {
558 return(rejection_term);
559 }
560
561 else if(particleDefinition->GetAtomicMass() > 4) // anything above Helium
562 {
563 G4double Z = particleDefinition->GetAtomicNumber();
564
565 G4double x = 100.*std::sqrt(beta2)/gpow->powA(Z, 2./3.);
566 G4double Zeffion = Z*(1.-G4Exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
567 rejection_term*=Zeffion*Zeffion;
568 }
569
570 else if (particleDefinition == alphaPlusPlusDef )
571 {
572 isHelium = true;
573 slaterEffectiveCharge[0]=0.;
574 slaterEffectiveCharge[1]=0.;
575 slaterEffectiveCharge[2]=0.;
576 sCoefficient[0]=0.;
577 sCoefficient[1]=0.;
578 sCoefficient[2]=0.;
579 }
580
581 else if (particleDefinition == alphaPlusDef )
582 {
583 isHelium = true;
584 slaterEffectiveCharge[0]=2.0;
585 // The following values are provided by M. Dingfelder (priv. comm)
586 slaterEffectiveCharge[1]=2.0;
587 slaterEffectiveCharge[2]=2.0;
588 //
589 sCoefficient[0]=0.7;
590 sCoefficient[1]=0.15;
591 sCoefficient[2]=0.15;
592 }
593
594 else if (particleDefinition == heliumDef )
595 {
596 isHelium = true;
597 slaterEffectiveCharge[0]=1.7;
598 slaterEffectiveCharge[1]=1.15;
599 slaterEffectiveCharge[2]=1.15;
600 sCoefficient[0]=0.5;
601 sCoefficient[1]=0.25;
602 sCoefficient[2]=0.25;
603 }
604
605 if (isHelium)
606 {
607
608 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
609
610 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
611 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
612 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
613
614 rejection_term*= zEff * zEff;
615 }
616
617 return (rejection_term);
618}
619
620
621//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
622
623
624G4double G4DNARuddIonisationExtendedModel::ProposedSampledEnergy(G4ParticleDefinition* particle,
625 G4double k,
626 G4int ionizationLevelIndex)
627{
628
629 const G4int j=ionizationLevelIndex;
630
631 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2;
632 //G4double alphaConst ;
633 G4double Bj_energy;
634
635 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
636 // Following values provided by M. Dingfelder (priv. comm)
637
638 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
639
640 if (j == 4)
641 {
642 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
643 A1 = 1.25;
644 B1 = 0.5;
645 C1 = 1.00;
646 D1 = 1.00;
647 E1 = 3.00;
648 A2 = 1.10;
649 B2 = 1.30;
650 C2 = 1.00;
651 D2 = 0.00;
652 //alphaConst = 0.66;
653 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
654 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
655 //---
656 }
657 else
658 {
659 //Data For Liquid Water from Dingfelder (Protons in Water)
660 A1 = 1.02;
661 B1 = 82.0;
662 C1 = 0.45;
663 D1 = -0.80;
664 E1 = 0.38;
665 A2 = 1.07;
666 //B2 = 14.6; From Ding Paper
667 // Value provided by M. Dingfelder (priv. comm)
668 B2 = 11.6;
669 //
670 C2 = 0.60;
671 D2 = 0.04;
672 //alphaConst = 0.64;
673
674 Bj_energy = Bj[ionizationLevelIndex];
675 }
676
677 G4double tau = 0.;
678 G4double A_ion = 0.;
679 tau = (electron_mass_c2 / particle->GetPDGMass()) * k;
680
681 A_ion = particle->GetAtomicMass();
682
683 G4double v2;
684 G4double beta2;
685 if((tau/MeV)<5.447761194e-2)
686 {
687 v2 = tau / Bj_energy;
688 beta2 = 2.*tau / electron_mass_c2;
689 }
690 // Relativistic
691 else
692 {
693 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
694 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
695 }
696
697 G4double v = std::sqrt(v2);
698 //G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
699 G4double L1 = (C1* gpow->powA(v,(D1))) / (1.+ E1*gpow->powA(v, (D1+4.)));
700 G4double L2 = C2*gpow->powA(v,(D2));
701 G4double H1 = (A1*G4Log(1.+v2)) / (v2+(B1/v2));
702 G4double H2 = (A2/v2) + (B2/(v2*v2));
703 G4double F1 = L1+H1;
704 G4double F2 = (L2*H2)/(L2+H2);
705
706 // ZF. generalized & relativistic version
707 G4double maximumEnergy;
708
709 //---- maximum kinetic energy , non relativistic ------
710 if( (k/MeV)/(particle->GetPDGMass()/MeV) <= 0.1 )
711 {
712 maximumEnergy = 4.* (electron_mass_c2 / particle->GetPDGMass()) * k;
713 }
714 //---- relativistic -----------------------------------
715 else
716 {
717 G4double gamma = 1./sqrt(1.-beta2);
718 maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
719 (1.+2.*gamma*(electron_mass_c2/particle->GetPDGMass())+pow(electron_mass_c2/particle->GetPDGMass(), 2.) );
720 }
721
722 //either it is transfered energy or secondary electron energy ...
723 //maximumEnergy-=Bj_energy;
724
725 //-----------------------------------------------------
726 G4double wmax = maximumEnergy/Bj_energy;
727 G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
728 c=1./c; //!!!!!!!!!!! manual calculus leads to c=1/c
729 G4double randVal = G4UniformRand();
730 G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
731 proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
732 // proposed_ws = -F1*c+2.*randVal-std::sqrt(proposed_ws);
733 proposed_ws/= ( F1*c + F2*c - 2.*randVal );
734 proposed_ws*=Bj_energy;
735
736 return(proposed_ws);
737
738}
739
740//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
741
742G4double G4DNARuddIonisationExtendedModel::S_1s(G4double t,
743 G4double energyTransferred,
744 G4double slaterEffectiveChg,
745 G4double shellNumber)
746{
747 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
748 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
749
750 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
751 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
752
753 return value;
754}
755
756//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
757
758G4double G4DNARuddIonisationExtendedModel::S_2s(G4double t,
759 G4double energyTransferred,
760 G4double slaterEffectiveChg,
761 G4double shellNumber)
762{
763 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
764 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
765
766 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
767 G4double value = 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
768
769 return value;
770
771}
772
773//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
774
775G4double G4DNARuddIonisationExtendedModel::S_2p(G4double t,
776 G4double energyTransferred,
777 G4double slaterEffectiveChg,
778 G4double shellNumber)
779{
780 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
781 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
782
783 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
784 G4double value = 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
785
786 return value;
787}
788
789//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
790
791G4double G4DNARuddIonisationExtendedModel::R(G4double t,
792 G4double energyTransferred,
793 G4double slaterEffectiveChg,
794 G4double shellNumber)
795{
796 // tElectron = m_electron / m_alpha * t
797 // Dingfelder, in Chattanooga 2005 proceedings, p 4
798
799 G4double tElectron = 0.511/3728. * t;
800 // The following values are provided by M. Dingfelder (priv. comm)
801 G4double H = 2.*13.60569172 * eV;
802 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
803
804 return value;
805}
806
807//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
808
809G4double G4DNARuddIonisationExtendedModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
810{
811 // ZF Shortened
812
813 if (particleDefinition == hydrogenDef && shell < 4)
814 {
815 G4double value = ((G4Log(k/eV)/gpow->logZ(10))-4.2)/0.5;
816 // The following values are provided by M. Dingfelder (priv. comm)
817 return((0.6/(1+G4Exp(value))) + 0.9);
818 }
819 else
820 {
821 return(1.);
822 }
823}
824
825//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
826
827G4int G4DNARuddIonisationExtendedModel::RandomSelect(G4double k)
828{
829 G4int level = 0;
830
831 // Retrieve data table corresponding to the current particle type
832
833 G4DNACrossSectionDataSet* table = currentTable;
834
835 if (table != nullptr)
836 {
837 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
838
839 const G4int n = (G4int)table->NumberOfComponents();
840 G4int i(n);
841 G4double value = 0.;
842
843 while (i>0)
844 {
845 --i;
846 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
847
848 value += valuesBuffer[i];
849 }
850
851 value *= G4UniformRand();
852
853 i = n;
854
855 while (i > 0)
856 {
857 --i;
858
859 if (valuesBuffer[i] > value)
860 {
861 delete[] valuesBuffer;
862 return i;
863 }
864 value -= valuesBuffer[i];
865 }
866
867 if (valuesBuffer) delete[] valuesBuffer;
868
869 }
870 return level;
871}
872
873//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
874
875G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(const G4Track& track )
876{
877 G4double sigma = 0.;
878
879 const G4DynamicParticle* particle = track.GetDynamicParticle();
880 G4double k = particle->GetKineticEnergy();
881
882 G4double lowLim = 0;
883 G4double highLim = 0;
884
885 const G4String& particleName = particle->GetDefinition()->GetParticleName();
886
887 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
888 pos1 = lowEnergyLimit.find(particleName);
889
890 if (pos1 != lowEnergyLimit.end())
891 {
892 lowLim = pos1->second;
893 }
894
895 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
896 pos2 = highEnergyLimit.find(particleName);
897
898 if (pos2 != highEnergyLimit.end())
899 {
900 highLim = pos2->second;
901 }
902
903 if (k >= lowLim && k <= highLim)
904 {
905 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
906 pos = tableData.find(particleName);
907
908 if (pos != tableData.end())
909 {
910 G4DNACrossSectionDataSet* table = pos->second;
911 if (table != 0)
912 {
913 sigma = table->FindValue(k);
914 }
915 }
916 else
917 {
918 G4Exception("G4DNARuddIonisationExtendedModel::PartialCrossSection","em0002",
919 FatalException,"Model not applicable to particle type.");
920 }
921 }
922
923 return sigma;
924}
925
926//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
927
928G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */)
929{
930 return 0;
931}
932
933//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
934
935G4ParticleDefinition* G4DNARuddIonisationExtendedModel::GetDNAIonParticleDefinition(const G4ParticleDefinition* particleDefinition)
936{
937 //for proton, hydrogen, alphas
938 if(particleDefinition->GetAtomicMass() <= 4)
939 {
940 return const_cast<G4ParticleDefinition*>(particleDefinition);
941 }
942 else{
943 auto instance = G4DNAGenericIonsManager::Instance();
944
945 auto PDGEncoding = particleDefinition->GetPDGEncoding();
946 if(PDGEncoding == 1000140280){
947 return instance->GetIon("silicon");
948 }else if(PDGEncoding == 1000260560){
949 return instance->GetIon("iron");
950 }else if(PDGEncoding == 1000080160){
951 return instance->GetIon("oxygen");
952 }else if(PDGEncoding == 1000070140){
953 return instance->GetIon("nitrogen");
954 }else if(PDGEncoding == 1000060120){
955 return instance->GetIon("carbon");
956 }
957 //if there is no DNA particle, get nullptr
958 return nullptr;
959 }
960}
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:50
eIonizedMolecule
@ eIonizedMolecule
Definition: G4DNAChemistryManager.hh:70
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
fStopButAlive
@ fStopButAlive
Definition: G4TrackStatus.hh:44
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
Z
const G4int Z[17]
Definition: G4WaterStopping.cc:51
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
C2
const double C2
Definition: MCGIDI_KalbachMann.cc:15
C1
#define C1
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:88
G4DNAChemistryManager::Instance
static G4DNAChemistryManager * Instance()
Definition: G4DNAChemistryManager.cc:132
G4DNAChemistryManager::CreateWaterMolecule
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Definition: G4DNAChemistryManager.cc:607
G4DNACrossSectionDataSet::FindValue
virtual G4double FindValue(G4double e, G4int componentId=0) const
Definition: G4DNACrossSectionDataSet.cc:483
G4DNACrossSectionDataSet::NumberOfComponents
virtual size_t NumberOfComponents(void) const
Definition: G4DNACrossSectionDataSet.hh:74
G4DNACrossSectionDataSet::GetComponent
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
Definition: G4DNACrossSectionDataSet.hh:68
G4DNACrossSectionDataSet::LoadData
virtual G4bool LoadData(const G4String &argFileName)
Definition: G4DNACrossSectionDataSet.cc:83
G4DNAGenericIonsManager::Instance
static G4DNAGenericIonsManager * Instance(void)
Definition: G4DNAGenericIonsManager.cc:35
G4DNAGenericIonsManager::GetIon
G4ParticleDefinition * GetIon(const G4String &name)
Definition: G4DNAGenericIonsManager.cc:45
G4DNAMolecularMaterial::GetNumMolPerVolTableFor
const std::vector< G4double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
Definition: G4DNAMolecularMaterial.cc:418
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition: G4DNAMolecularMaterial.cc:80
G4DNARuddIonisationExtendedModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4DNARuddIonisationExtendedModel.cc:351
G4DNARuddIonisationExtendedModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
Definition: G4DNARuddIonisationExtendedModel.cc:297
G4DNARuddIonisationExtendedModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4DNARuddIonisationExtendedModel.cc:118
G4DNARuddIonisationExtendedModel::G4DNARuddIonisationExtendedModel
G4DNARuddIonisationExtendedModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationExtendedModel")
Definition: G4DNARuddIonisationExtendedModel.cc:55
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:87
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.hh:389
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:166
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:255
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:691
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:106
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
Definition: G4ParticleChangeForGamma.hh:120
G4ParticleDefinition::GetAtomicNumber
G4int GetAtomicNumber() const
G4ParticleDefinition::GetAtomicMass
G4int GetAtomicMass() const
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:113
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:109
G4Pow
Definition: G4Pow.hh:49
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4Pow::logZ
G4double logZ(G4int Z) const
Definition: G4Pow.hh:137
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:87
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4VAtomDeexcitation::GetAtomicShell
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4VAtomDeexcitation::GenerateParticles
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
Definition: G4VAtomDeexcitation.cc:234
G4VEMDataSet::FindValue
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4VEmAngularDistribution::SampleDirectionForShell
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
Definition: G4VEmAngularDistribution.cc:68
G4VEmModel::GetAngularDistribution
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:600
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:124
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:802
G4VEmModel::SetAngularDistribution
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:607
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VParticleChange::GetCurrentTrack
const G4Track * GetCurrentTrack() const
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62