| Geant4 Cross Reference |
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
42 static G4Pow * gpow = G4Pow::GetInstance();
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44
45 using namespace std;
46
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
49 G4DNARuddIonisationModel::G4DNARuddIonisationModel(const G4ParticleDefinition*,
50 const G4String& nam) :
51 G4VEmModel(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
93 G4DNARuddIonisationModel::~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
112 void 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( 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
272 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
273
274 //
275
276 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
277
278 if (isInitialised)
279 { return;}
280 fParticleChangeForGamma = GetParticleChangeForGamma();
281 isInitialised = true;
282
283 }
284
285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286
287 G4double G4DNARuddIonisationModel::CrossSectionPerVolume(const G4Material* material,
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
401 void 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
504 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
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 {
555 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
556 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
557 }
558 else
559 {
560 fParticleChangeForGamma->SetProposedKineticEnergy(k);
561 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
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();
578 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
579 ionizationShell,
580 theIncomingTrack);
581 }
582
583 // SI - not useful since low energy of model is 0 eV
584
585 if (k < lowLim)
586 {
587 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
588 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
589 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
590 }
591
592 }
593
594 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
595
596 G4double 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......
682 G4double 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
912 G4double 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
928 G4double 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
946 G4double 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
964 G4double 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
983 G4double 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
1002 G4int 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",
1060 FatalException,
1061 "Model not applicable to particle type.");
1062 }
1063
1064 return level;
1065 }
1066
1067 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1068
1069 G4double 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",
1113 FatalException,
1114 "Model not applicable to particle type.");
1115 }
1116 }
1117
1118 return sigma;
1119 }
1120
1121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1122
1123 G4double G4DNARuddIonisationModel::Sum(G4double /* energy */,
1124 const G4String& /* particle */)
1125 {
1126 return 0;
1127 }
1128
1129