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