Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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