Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNARuddIonisationModel.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
30#include "G4SystemOfUnits.hh"
32#include "G4LossTableManager.hh"
35#include "G4DNARuddAngle.hh"
36#include "G4DeltaAngle.hh"
37#include "G4Exp.hh"
38#include "G4Pow.hh"
39#include "G4Log.hh"
40#include "G4Alpha.hh"
41
42static G4Pow * gpow = G4Pow::GetInstance();
43//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44
45using namespace std;
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
50 const G4String& nam) :
51G4VEmModel(nam)
52{
53 slaterEffectiveCharge[0] = 0.;
54 slaterEffectiveCharge[1] = 0.;
55 slaterEffectiveCharge[2] = 0.;
56 sCoefficient[0] = 0.;
57 sCoefficient[1] = 0.;
58 sCoefficient[2] = 0.;
59
60 lowEnergyLimitForZ1 = 0 * eV;
61 lowEnergyLimitForZ2 = 0 * eV;
62 lowEnergyLimitOfModelForZ1 = 100 * eV;
63 lowEnergyLimitOfModelForZ2 = 1 * keV;
64 killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
65 killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
66
67 verboseLevel = 0;
68 // Verbosity scale:
69 // 0 = nothing
70 // 1 = warning for energy non-conservation
71 // 2 = details of energy budget
72 // 3 = calculation of cross sections, file openings, sampling of atoms
73 // 4 = entering in methods
74
75 if (verboseLevel > 0)
76 {
77 G4cout << "Rudd ionisation model is constructed " << G4endl;
78 }
79
80 // Define default angular generator
82
83 // Mark this model as "applicable" for atomic deexcitation
85
86 // Selection of stationary mode
87
88 statCode = false;
89}
90
91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92
94{
95 // Cross section
96
97 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
98 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
99 {
100 G4DNACrossSectionDataSet* table = pos->second;
101 delete table;
102 }
103
104 // The following removal is forbidden since G4VEnergyLossmodel takes care of deletion
105 // Coverity however will signal this as an error
106 // if (fAtomDeexcitation) {delete fAtomDeexcitation;}
107
108}
109
110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111
113 const G4DataVector& /*cuts*/)
114{
115
116 if (verboseLevel > 3)
117 {
118 G4cout << "Calling G4DNARuddIonisationModel::Initialise()" << G4endl;
119 }
120
121 // Energy limits
122
123 G4String fileProton("dna/sigma_ionisation_p_rudd");
124 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
125 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
126 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
127 G4String fileHelium("dna/sigma_ionisation_he_rudd");
128
129 G4DNAGenericIonsManager *instance;
131 protonDef = G4Proton::ProtonDefinition();
132 hydrogenDef = instance->GetIon("hydrogen");
133 alphaPlusPlusDef = G4Alpha::Alpha();
134 alphaPlusDef = instance->GetIon("alpha+");
135 heliumDef = instance->GetIon("helium");
136
137 G4String proton;
138 G4String hydrogen;
139 G4String alphaPlusPlus;
140 G4String alphaPlus;
141 G4String helium;
142
143 G4double scaleFactor = 1 * m*m;
144
145 // LIMITS AND DATA
146
147 // ********************************************************
148
149 proton = protonDef->GetParticleName();
150 tableFile[proton] = fileProton;
151
152 lowEnergyLimit[proton] = lowEnergyLimitForZ1;
153 highEnergyLimit[proton] = 500. * keV;
154
155 // Cross section
156
157 auto tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
158 eV,
159 scaleFactor );
160 tableProton->LoadData(fileProton);
161 tableData[proton] = tableProton;
162
163 // ********************************************************
164
165 hydrogen = hydrogenDef->GetParticleName();
166 tableFile[hydrogen] = fileHydrogen;
167
168 lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
169 highEnergyLimit[hydrogen] = 100. * MeV;
170
171 // Cross section
172
173 auto tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
174 eV,
175 scaleFactor );
176 tableHydrogen->LoadData(fileHydrogen);
177
178 tableData[hydrogen] = tableHydrogen;
179
180 // ********************************************************
181
182 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
183 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
184
185 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
186 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
187
188 // Cross section
189
190 auto tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
191 eV,
192 scaleFactor );
193 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
194
195 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
196
197 // ********************************************************
198
199 alphaPlus = alphaPlusDef->GetParticleName();
200 tableFile[alphaPlus] = fileAlphaPlus;
201
202 lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
203 highEnergyLimit[alphaPlus] = 400. * MeV;
204
205 // Cross section
206
207 auto tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
208 eV,
209 scaleFactor );
210 tableAlphaPlus->LoadData(fileAlphaPlus);
211 tableData[alphaPlus] = tableAlphaPlus;
212
213 // ********************************************************
214
215 helium = heliumDef->GetParticleName();
216 tableFile[helium] = fileHelium;
217
218 lowEnergyLimit[helium] = lowEnergyLimitForZ2;
219 highEnergyLimit[helium] = 400. * MeV;
220
221 // Cross section
222
223 auto tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
224 eV,
225 scaleFactor );
226 tableHelium->LoadData(fileHelium);
227 tableData[helium] = tableHelium;
228
229 //
230
231 if (particle==protonDef)
232 {
233 SetLowEnergyLimit(lowEnergyLimit[proton]);
234 SetHighEnergyLimit(highEnergyLimit[proton]);
235 }
236
237 if (particle==hydrogenDef)
238 {
239 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
240 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
241 }
242
243 if (particle==heliumDef)
244 {
245 SetLowEnergyLimit(lowEnergyLimit[helium]);
246 SetHighEnergyLimit(highEnergyLimit[helium]);
247 }
248
249 if (particle==alphaPlusDef)
250 {
251 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
252 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
253 }
254
255 if (particle==alphaPlusPlusDef)
256 {
257 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
258 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
259 }
260
261 if( verboseLevel>0 )
262 {
263 G4cout << "Rudd ionisation model is initialized " << G4endl
264 << "Energy range: "
265 << LowEnergyLimit() / eV << " eV - "
266 << HighEnergyLimit() / keV << " keV for "
267 << particle->GetParticleName()
268 << G4endl;
269 }
270
271 // Initialize water density pointer
273
274 //
275
276 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
277
278 if (isInitialised)
279 { return;}
281 isInitialised = true;
282
283}
284
285//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286
288 const G4ParticleDefinition* particleDefinition,
289 G4double k,
290 G4double,
291 G4double)
292{
293 if (verboseLevel > 3)
294 {
295 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationModel"
296 << G4endl;
297 }
298
299 // Calculate total cross section for model
300
301 if (
302 particleDefinition != protonDef
303 &&
304 particleDefinition != hydrogenDef
305 &&
306 particleDefinition != alphaPlusPlusDef
307 &&
308 particleDefinition != alphaPlusDef
309 &&
310 particleDefinition != heliumDef
311 )
312
313 return 0;
314
315 G4double lowLim = 0;
316
317 if ( particleDefinition == protonDef
318 || particleDefinition == hydrogenDef
319 )
320
321 lowLim = lowEnergyLimitOfModelForZ1;
322
323 if ( particleDefinition == alphaPlusPlusDef
324 || particleDefinition == alphaPlusDef
325 || particleDefinition == heliumDef
326 )
327
328 lowLim = lowEnergyLimitOfModelForZ2;
329
330 G4double highLim = 0;
331 G4double sigma=0;
332
333 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
334
335 const G4String& particleName = particleDefinition->GetParticleName();
336
337 // SI - the following is useless since lowLim is already defined
338 /*
339 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
340 pos1 = lowEnergyLimit.find(particleName);
341
342 if (pos1 != lowEnergyLimit.end())
343 {
344 lowLim = pos1->second;
345 }
346 */
347
348 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
349 pos2 = highEnergyLimit.find(particleName);
350
351 if (pos2 != highEnergyLimit.end())
352 {
353 highLim = pos2->second;
354 }
355
356 if (k <= highLim)
357 {
358 //SI : XS must not be zero otherwise sampling of secondaries method ignored
359
360 if (k < lowLim) k = lowLim;
361
362 //
363
364 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
365 pos = tableData.find(particleName);
366
367 if (pos != tableData.end())
368 {
369 G4DNACrossSectionDataSet* table = pos->second;
370 if (table != nullptr)
371 {
372 sigma = table->FindValue(k);
373 }
374 }
375 else
376 {
377 G4Exception("G4DNARuddIonisationModel::CrossSectionPerVolume","em0002",
378 FatalException,"Model not applicable to particle type.");
379 }
380
381 }
382
383 if (verboseLevel > 2)
384 {
385 G4cout << "__________________________________" << G4endl;
386 G4cout << "G4DNARuddIonisationModel - XS INFO START" << G4endl;
387 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
388 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
389 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
390 //G4cout << " - Cross section per water molecule (cm^-1)="
391 //<< sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
392 G4cout << "G4DNARuddIonisationModel - XS INFO END" << G4endl;
393 }
394
395 return sigma*waterDensity;
396
397}
398
399//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400
401void G4DNARuddIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
402 const G4MaterialCutsCouple* couple,
403 const G4DynamicParticle* particle,
404 G4double,
405 G4double)
406{
407 if (verboseLevel > 3)
408 {
409 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationModel"
410 << G4endl;
411 }
412
413 G4double lowLim = 0;
414 G4double highLim = 0;
415
416 if ( particle->GetDefinition() == protonDef
417 || particle->GetDefinition() == hydrogenDef
418 )
419
420 lowLim = killBelowEnergyForZ1;
421
422 if ( particle->GetDefinition() == alphaPlusPlusDef
423 || particle->GetDefinition() == alphaPlusDef
424 || particle->GetDefinition() == heliumDef
425 )
426
427 lowLim = killBelowEnergyForZ2;
428
429 G4double k = particle->GetKineticEnergy();
430
431 const G4String& particleName = particle->GetDefinition()->GetParticleName();
432
433 // SI - the following is useless since lowLim is already defined
434 /*
435 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
436 pos1 = lowEnergyLimit.find(particleName);
437
438 if (pos1 != lowEnergyLimit.end())
439 {
440 lowLim = pos1->second;
441 }
442 */
443
444 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
445 pos2 = highEnergyLimit.find(particleName);
446
447 if (pos2 != highEnergyLimit.end())
448 {
449 highLim = pos2->second;
450 }
451
452 if (k >= lowLim && k <= highLim)
453 {
454 G4ParticleDefinition* definition = particle->GetDefinition();
455 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
456 /*
457 G4double particleMass = definition->GetPDGMass();
458 G4double totalEnergy = k + particleMass;
459 G4double pSquare = k*(totalEnergy+particleMass);
460 G4double totalMomentum = std::sqrt(pSquare);
461 */
462
463 G4int ionizationShell = RandomSelect(k,particleName);
464
465 G4double bindingEnergy = 0;
466 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
467
468 //SI: additional protection if tcs interpolation method is modified
469 if (k<bindingEnergy) return;
470 //
471
472 // SI - For atom. deexc. tagging - 23/05/2017
473 G4int Z = 8;
474 //
475
476 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
477
478 G4ThreeVector deltaDirection =
479 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
480 Z, ionizationShell,
481 couple->GetMaterial());
482
483 auto dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
484 fvect->push_back(dp);
485
486 // Ignored for ions on electrons
487 /*
488 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
489
490 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
491 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
492 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
493 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
494 finalPx /= finalMomentum;
495 finalPy /= finalMomentum;
496 finalPz /= finalMomentum;
497
498 G4ThreeVector direction;
499 direction.set(finalPx,finalPy,finalPz);
500
501 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
502 */
503
505
506 // sample deexcitation
507 // here we assume that H_{2}O electronic levels are the same of Oxigen.
508 // this can be considered true with a rough 10% error in energy on K-shell,
509
510 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
511 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
512
513 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
514
515 // SI: only atomic deexcitation from K shell is considered
516 if((fAtomDeexcitation != nullptr) && ionizationShell == 4)
517 {
518 const G4AtomicShell* shell
519 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
520 secNumberInit = fvect->size();
521 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
522 secNumberFinal = fvect->size();
523
524 if(secNumberFinal > secNumberInit)
525 {
526 for (size_t i=secNumberInit; i<secNumberFinal; ++i)
527 {
528 //Check if there is enough residual energy
529 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
530 {
531 //Ok, this is a valid secondary: keep it
532 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
533 }
534 else
535 {
536 //Invalid secondary: not enough energy to create it!
537 //Keep its energy in the local deposit
538 delete (*fvect)[i];
539 (*fvect)[i]=nullptr;
540 }
541 }
542 }
543
544 }
545
546 //This should never happen
547 if(bindingEnergy < 0.0)
548 G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
549 "em2050",FatalException,"Negative local energy deposit");
550
551 //bindingEnergy has been decreased
552 //by the amount of energy taken away by deexc. products
553 if (!statCode)
554 {
557 }
558 else
559 {
562 }
563
564 // debug
565 // k-scatteredEnergy-secondaryKinetic-deexSecEnergy = k-(k-bindingEnergy-secondaryKinetic)-secondaryKinetic-deexSecEnergy =
566 // = k-k+bindingEnergy+secondaryKinetic-secondaryKinetic-deexSecEnergy=
567 // = bindingEnergy-deexSecEnergy
568 // SO deexSecEnergy=0 => LocalEnergyDeposit = bindingEnergy
569
570 // TEST //////////////////////////
571 // if (secondaryKinetic<0) abort();
572 // if (scatteredEnergy<0) abort();
573 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
574 // if (k-scatteredEnergy<0) abort();
575 /////////////////////////////////
576
577 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
579 ionizationShell,
580 theIncomingTrack);
581 }
582
583 // SI - not useful since low energy of model is 0 eV
584
585 if (k < lowLim)
586 {
590 }
591
592}
593
594//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
595
596G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
597 G4double k,
598 G4int shell)
599{
600 G4double maximumKineticEnergyTransfer = 0.;
601
602 if (particleDefinition == protonDef
603 || particleDefinition == hydrogenDef)
604 {
605 maximumKineticEnergyTransfer = 4. * (electron_mass_c2 / proton_mass_c2) * k;
606 }
607
608 else if (particleDefinition == heliumDef
609 || particleDefinition == alphaPlusDef
610 || particleDefinition == alphaPlusPlusDef)
611 {
612 maximumKineticEnergyTransfer = 4. * (0.511 / 3728) * k;
613 }
614
615 G4double crossSectionMaximum = 0.;
616
617 for (G4double value = waterStructure.IonisationEnergy(shell);
618 value <= 5. * waterStructure.IonisationEnergy(shell) && k >= value;
619 value += 0.1 * eV)
620 {
621 G4double differentialCrossSection =
622 DifferentialCrossSection(particleDefinition, k, value, shell);
623 if (differentialCrossSection >= crossSectionMaximum)
624 crossSectionMaximum = differentialCrossSection;
625 }
626
627 G4double secElecKinetic = 0.;
628
629 do
630 {
631 secElecKinetic = G4UniformRand()* maximumKineticEnergyTransfer;
632 }while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
633 k,
634 secElecKinetic+waterStructure.IonisationEnergy(shell),
635 shell));
636
637 return (secElecKinetic);
638}
639
640//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
641
642// The following section is not used anymore but is kept for memory
643// GetAngularDistribution()->SampleDirectionForShell is used instead
644
645/*
646 void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
647 G4double k,
648 G4double secKinetic,
649 G4double & cosTheta,
650 G4double & phi )
651 {
652 G4DNAGenericIonsManager *instance;
653 instance = G4DNAGenericIonsManager::Instance();
654
655 G4double maxSecKinetic = 0.;
656
657 if (particleDefinition == protonDef
658 || particleDefinition == hydrogenDef)
659 {
660 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
661 }
662
663 else if (particleDefinition == heliumDef
664 || particleDefinition == alphaPlusDef
665 || particleDefinition == alphaPlusPlusDef)
666 {
667 maxSecKinetic = 4.* (0.511 / 3728) * k;
668 }
669
670 phi = twopi * G4UniformRand();
671
672 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
673
674 // Restriction below 100 eV from Emfietzoglou (2000)
675
676 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
677 else cosTheta = (2.*G4UniformRand())-1.;
678
679 }
680 */
681//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
682G4double G4DNARuddIonisationModel::DifferentialCrossSection(G4ParticleDefinition* particleDefinition,
683 G4double k,
684 G4double energyTransfer,
685 G4int ionizationLevelIndex)
686{
687 // Shells ids are 0 1 2 3 4 (4 is k shell)
688 // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
689 // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
690 //
691 // ds S F1(nu) + w * F2(nu)
692 // ---- = G(k) * ---- -------------------------------------------
693 // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
694 //
695 // w is the secondary electron kinetic Energy in eV
696 //
697 // All the other parameters can be found in Rudd's Papers
698 //
699 // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
700 // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
701 //
702
703 const G4int j = ionizationLevelIndex;
704
705 G4double A1;
706 G4double B1;
707 G4double C1;
708 G4double D1;
709 G4double E1;
710 G4double A2;
711 G4double B2;
712 G4double C2;
713 G4double D2;
714 G4double alphaConst;
715
716 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
717 // The following values are provided by M. dingfelder (priv. comm)
718 const G4double Bj[5] = { 12.60 * eV, 14.70 * eV, 18.40 * eV, 32.20 * eV, 540
719 * eV };
720
721 if (j == 4)
722 {
723 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
724 A1 = 1.25;
725 B1 = 0.5;
726 C1 = 1.00;
727 D1 = 1.00;
728 E1 = 3.00;
729 A2 = 1.10;
730 B2 = 1.30;
731 C2 = 1.00;
732 D2 = 0.00;
733 alphaConst = 0.66;
734 } else
735 {
736 //Data For Liquid Water from Dingfelder (Protons in Water)
737 A1 = 1.02;
738 B1 = 82.0;
739 C1 = 0.45;
740 D1 = -0.80;
741 E1 = 0.38;
742 A2 = 1.07;
743 // Value provided by M. Dingfelder (priv. comm)
744 B2 = 11.6;
745 //
746 C2 = 0.60;
747 D2 = 0.04;
748 alphaConst = 0.64;
749 }
750
751 const G4double n = 2.;
752 const G4double Gj[5] = { 0.99, 1.11, 1.11, 0.52, 1. };
753
754 G4double wBig = (energyTransfer
755 - waterStructure.IonisationEnergy(ionizationLevelIndex));
756 if (wBig < 0)
757 return 0.;
758
759 G4double w = wBig / Bj[ionizationLevelIndex];
760 // Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
761 if (j == 4)
762 w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
763
764 G4double Ry = 13.6 * eV;
765
766 G4double tau = 0.;
767
768 G4bool isProtonOrHydrogen = false;
769 G4bool isHelium = false;
770
771 if (particleDefinition == protonDef
772 || particleDefinition == hydrogenDef)
773 {
774 isProtonOrHydrogen = true;
775 tau = (electron_mass_c2 / proton_mass_c2) * k;
776 }
777
778 else if (particleDefinition == heliumDef
779 || particleDefinition == alphaPlusDef
780 || particleDefinition == alphaPlusPlusDef)
781 {
782 isHelium = true;
783 tau = (0.511 / 3728.) * k;
784 }
785
786 G4double S = 4. * pi * Bohr_radius * Bohr_radius * n
787 * gpow->powN((Ry / Bj[ionizationLevelIndex]), 2);
788 if (j == 4)
789 S = 4. * pi * Bohr_radius * Bohr_radius * n
790 * gpow->powN((Ry / waterStructure.IonisationEnergy(ionizationLevelIndex)),
791 2);
792
793 G4double v2 = tau / Bj[ionizationLevelIndex];
794 if (j == 4)
795 v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
796
797 G4double v = std::sqrt(v2);
798 G4double wc = 4. * v2 - 2. * v - (Ry / (4. * Bj[ionizationLevelIndex]));
799 if (j == 4)
800 wc = 4. * v2 - 2. * v
801 - (Ry / (4. * waterStructure.IonisationEnergy(ionizationLevelIndex)));
802
803 G4double L1 = (C1 * gpow->powA(v, (D1))) / (1. + E1 * gpow->powA(v, (D1 + 4.)));
804 G4double L2 = C2 * gpow->powA(v, (D2));
805 G4double H1 = (A1 * G4Log(1. + v2)) / (v2 + (B1 / v2));
806 G4double H2 = (A2 / v2) + (B2 / (v2 * v2));
807
808 G4double F1 = L1 + H1;
809 G4double F2 = (L2 * H2) / (L2 + H2);
810
811 G4double sigma =
812 CorrectionFactor(particleDefinition, k) * Gj[j]
813 * (S / Bj[ionizationLevelIndex])
814 * ((F1 + w * F2)
815 / (gpow->powN((1. + w), 3)
816 * (1. + G4Exp(alphaConst * (w - wc) / v))));
817
818 if (j == 4)
819 sigma = CorrectionFactor(particleDefinition, k) * Gj[j]
820 * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
821 * ((F1 + w * F2)
822 / (gpow->powN((1. + w), 3)
823 * (1. + G4Exp(alphaConst * (w - wc) / v))));
824
825 if ((particleDefinition == hydrogenDef)
826 && (ionizationLevelIndex == 4))
827 {
828 // sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
829 sigma = Gj[j] * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
830 * ((F1 + w * F2)
831 / (gpow->powN((1. + w), 3)
832 * (1. + G4Exp(alphaConst * (w - wc) / v))));
833 }
834
835 // if ( particleDefinition == protonDef
836 // || particleDefinition == hydrogenDef
837 // )
838
839 if (isProtonOrHydrogen)
840 {
841 return (sigma);
842 }
843
844 if (particleDefinition == alphaPlusPlusDef)
845 {
846 slaterEffectiveCharge[0] = 0.;
847 slaterEffectiveCharge[1] = 0.;
848 slaterEffectiveCharge[2] = 0.;
849 sCoefficient[0] = 0.;
850 sCoefficient[1] = 0.;
851 sCoefficient[2] = 0.;
852 }
853
854 else if (particleDefinition == alphaPlusDef)
855 {
856 slaterEffectiveCharge[0] = 2.0;
857 // The following values are provided by M. Dingfelder (priv. comm)
858 slaterEffectiveCharge[1] = 2.0;
859 slaterEffectiveCharge[2] = 2.0;
860 //
861 sCoefficient[0] = 0.7;
862 sCoefficient[1] = 0.15;
863 sCoefficient[2] = 0.15;
864 }
865
866 else if (particleDefinition == heliumDef)
867 {
868 slaterEffectiveCharge[0] = 1.7;
869 slaterEffectiveCharge[1] = 1.15;
870 slaterEffectiveCharge[2] = 1.15;
871 sCoefficient[0] = 0.5;
872 sCoefficient[1] = 0.25;
873 sCoefficient[2] = 0.25;
874 }
875
876 // if ( particleDefinition == heliumDef
877 // || particleDefinition == alphaPlusDef
878 // || particleDefinition == alphaPlusPlusDef
879 // )
880 if (isHelium)
881 {
882 sigma = Gj[j] * (S / Bj[ionizationLevelIndex])
883 * ((F1 + w * F2)
884 / (gpow->powN((1. + w), 3)
885 * (1. + G4Exp(alphaConst * (w - wc) / v))));
886
887 if (j == 4)
888 sigma = Gj[j]
889 * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
890 * ((F1 + w * F2)
891 / (gpow->powN((1. + w), 3)
892 * (1. + G4Exp(alphaConst * (w - wc) / v))));
893
894 G4double zEff = particleDefinition->GetPDGCharge() / eplus
895 + particleDefinition->GetLeptonNumber();
896
897 zEff -= (sCoefficient[0]
898 * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.)
899 + sCoefficient[1]
900 * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.)
901 + sCoefficient[2]
902 * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.));
903
904 return zEff * zEff * sigma;
905 }
906
907 return 0;
908}
909
910//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
911
912G4double G4DNARuddIonisationModel::S_1s(G4double t,
913 G4double energyTransferred,
914 G4double slaterEffectiveChg,
915 G4double shellNumber)
916{
917 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
918 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
919
920 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
921 G4double value = 1. - G4Exp(-2 * r) * ((2. * r + 2.) * r + 1.);
922
923 return value;
924}
925
926//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
927
928G4double G4DNARuddIonisationModel::S_2s(G4double t,
929 G4double energyTransferred,
930 G4double slaterEffectiveChg,
931 G4double shellNumber)
932{
933 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
934 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
935
936 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
937 G4double value = 1.
938 - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
939
940 return value;
941
942}
943
944//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
945
946G4double G4DNARuddIonisationModel::S_2p(G4double t,
947 G4double energyTransferred,
948 G4double slaterEffectiveChg,
949 G4double shellNumber)
950{
951 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
952 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
953
954 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
955 G4double value = 1.
956 - G4Exp(-2 * r)
957 * ((((2. / 3. * r + 4. / 3.) * r + 2.) * r + 2.) * r + 1.);
958
959 return value;
960}
961
962//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
963
964G4double G4DNARuddIonisationModel::R(G4double t,
965 G4double energyTransferred,
966 G4double slaterEffectiveChg,
967 G4double shellNumber)
968{
969 // tElectron = m_electron / m_alpha * t
970 // Dingfelder, in Chattanooga 2005 proceedings, p 4
971
972 G4double tElectron = 0.511 / 3728. * t;
973 // The following values are provided by M. Dingfelder (priv. comm)
974 G4double H = 2. * 13.60569172 * eV;
975 G4double value = std::sqrt(2. * tElectron / H) / (energyTransferred / H)
976 * (slaterEffectiveChg / shellNumber);
977
978 return value;
979}
980
981//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
982
983G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition,
984 G4double k)
985{
986
987 if (particleDefinition == G4Proton::Proton())
988 {
989 return (1.);
990 }
991 if (particleDefinition == hydrogenDef)
992 {
993 G4double value = (G4Log(k / eV)/gpow->logZ(10) - 4.2) / 0.5;
994 // The following values are provided by M. Dingfelder (priv. comm)
995 return ((0.6 / (1 + G4Exp(value))) + 0.9);
996 }
997 return (1.);
998}
999
1000//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1001
1002G4int G4DNARuddIonisationModel::RandomSelect(G4double k,
1003 const G4String& particle)
1004{
1005
1006 // BEGIN PART 1/2 OF ELECTRON CORRECTION
1007
1008 // add ONE or TWO electron-water ionisation for alpha+ and helium
1009
1010 G4int level = 0;
1011
1012 // Retrieve data table corresponding to the current particle type
1013
1014 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1015 pos = tableData.find(particle);
1016
1017 if (pos != tableData.end())
1018 {
1019 G4DNACrossSectionDataSet* table = pos->second;
1020
1021 if (table != nullptr)
1022 {
1023 auto valuesBuffer = new G4double[table->NumberOfComponents()];
1024
1025 const auto n = (G4int)table->NumberOfComponents();
1026 G4int i(n);
1027 G4double value = 0.;
1028
1029 while (i > 0)
1030 {
1031 --i;
1032 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1033 value += valuesBuffer[i];
1034 }
1035
1036 value *= G4UniformRand();
1037
1038 i = n;
1039
1040 while (i > 0)
1041 {
1042 --i;
1043
1044 if (valuesBuffer[i] > value)
1045 {
1046 delete[] valuesBuffer;
1047 return i;
1048 }
1049 value -= valuesBuffer[i];
1050 }
1051
1052
1053 delete[] valuesBuffer;
1054
1055 }
1056 } else
1057 {
1058 G4Exception("G4DNARuddIonisationModel::RandomSelect",
1059 "em0002",
1061 "Model not applicable to particle type.");
1062 }
1063
1064 return level;
1065}
1066
1067//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1068
1069G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track)
1070{
1071 G4double sigma = 0.;
1072
1073 const G4DynamicParticle* particle = track.GetDynamicParticle();
1074 G4double k = particle->GetKineticEnergy();
1075
1076 G4double lowLim = 0;
1077 G4double highLim = 0;
1078
1079 const G4String& particleName = particle->GetDefinition()->GetParticleName();
1080
1081 std::map<G4String, G4double, std::less<G4String> >::iterator pos1;
1082 pos1 = lowEnergyLimit.find(particleName);
1083
1084 if (pos1 != lowEnergyLimit.end())
1085 {
1086 lowLim = pos1->second;
1087 }
1088
1089 std::map<G4String, G4double, std::less<G4String> >::iterator pos2;
1090 pos2 = highEnergyLimit.find(particleName);
1091
1092 if (pos2 != highEnergyLimit.end())
1093 {
1094 highLim = pos2->second;
1095 }
1096
1097 if (k >= lowLim && k <= highLim)
1098 {
1099 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1100 pos = tableData.find(particleName);
1101
1102 if (pos != tableData.end())
1103 {
1104 G4DNACrossSectionDataSet* table = pos->second;
1105 if (table != nullptr)
1106 {
1107 sigma = table->FindValue(k);
1108 }
1109 } else
1110 {
1111 G4Exception("G4DNARuddIonisationModel::PartialCrossSection",
1112 "em0002",
1114 "Model not applicable to particle type.");
1115 }
1116 }
1117
1118 return sigma;
1119}
1120
1121//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1122
1123G4double G4DNARuddIonisationModel::Sum(G4double /* energy */,
1124 const G4String& /* particle */)
1125{
1126 return 0;
1127}
1128
@ eIonizedMolecule
G4double S(G4double temp)
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:180
G4double G4Log(G4double x)
Definition G4Log.hh:227
@ fStopAndKill
double G4double
Definition G4Types.hh:83
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
const double C2
#define C1
#define G4UniformRand()
Definition Randomize.hh:52
static G4Alpha * Alpha()
Definition G4Alpha.cc:83
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
size_t NumberOfComponents() const override
const G4VEMDataSet * GetComponent(G4int componentId) const override
G4double FindValue(G4double e, G4int componentId=0) const override
static G4DNAGenericIonsManager * Instance()
G4ParticleDefinition * GetIon(const G4String &name)
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...
static G4DNAMolecularMaterial * Instance()
G4ParticleChangeForGamma * fParticleChangeForGamma
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax) override
G4DNARuddIonisationModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="DNARuddIonisationModel")
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
Definition G4Electron.cc:91
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
std::size_t GetIndex() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
const G4String & GetParticleName() const
Definition G4Pow.hh:49
static G4Pow * GetInstance()
Definition G4Pow.cc:41
G4double logZ(G4int Z) const
Definition G4Pow.hh:137
G4double powN(G4double x, G4int n) const
Definition G4Pow.cc:162
G4double powA(G4double A, G4double y) const
Definition G4Pow.hh:230
static G4Proton * ProtonDefinition()
Definition G4Proton.cc:85
static G4Proton * Proton()
Definition G4Proton.cc:90
const G4DynamicParticle * GetDynamicParticle() const
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
void SetHighEnergyLimit(G4double)
G4VEmAngularDistribution * GetAngularDistribution()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void SetLowEnergyLimit(G4double)
void SetDeexcitationFlag(G4bool val)
void SetAngularDistribution(G4VEmAngularDistribution *)
void ProposeTrackStatus(G4TrackStatus status)
const G4Track * GetCurrentTrack() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)