| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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........oo
44
45 using namespace std;
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 G4DNARuddIonisationModel::G4DNARuddIonisationM
50
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 = lowEnergyLimitOfModel
65 killBelowEnergyForZ2 = lowEnergyLimitOfModel
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 o
73 // 4 = entering in methods
74
75 if (verboseLevel > 0)
76 {
77 G4cout << "Rudd ionisation model is constr
78 }
79
80 // Define default angular generator
81 SetAngularDistribution(new G4DNARuddAngle())
82
83 // Mark this model as "applicable" for atomi
84 SetDeexcitationFlag(true);
85
86 // Selection of stationary mode
87
88 statCode = false;
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oo
92
93 G4DNARuddIonisationModel::~G4DNARuddIonisation
94 {
95 // Cross section
96
97 std::map<G4String, G4DNACrossSectionDataSet*
98 for (pos = tableData.begin(); pos != tableDa
99 {
100 G4DNACrossSectionDataSet* table = pos->sec
101 delete table;
102 }
103
104 // The following removal is forbidden since
105 // Coverity however will signal this as an e
106 // if (fAtomDeexcitation) {delete fAtomDeex
107
108 }
109
110 //....oooOO0OOooo........oooOO0OOooo........oo
111
112 void G4DNARuddIonisationModel::Initialise(cons
113 cons
114 {
115
116 if (verboseLevel > 3)
117 {
118 G4cout << "Calling G4DNARuddIonisationMode
119 }
120
121 // Energy limits
122
123 G4String fileProton("dna/sigma_ionisation_p_
124 G4String fileHydrogen("dna/sigma_ionisation_
125 G4String fileAlphaPlusPlus("dna/sigma_ionisa
126 G4String fileAlphaPlus("dna/sigma_ionisation
127 G4String fileHelium("dna/sigma_ionisation_he
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 G4DNACrossSectionDat
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] = lowEnergyLimitFor
169 highEnergyLimit[hydrogen] = 100. * MeV;
170
171 // Cross section
172
173 auto tableHydrogen = new G4DNACrossSectionD
174 eV,
175 scaleFactor );
176 tableHydrogen->LoadData(fileHydrogen);
177
178 tableData[hydrogen] = tableHydrogen;
179
180 // *****************************************
181
182 alphaPlusPlus = alphaPlusPlusDef->GetParticl
183 tableFile[alphaPlusPlus] = fileAlphaPlusPlus
184
185 lowEnergyLimit[alphaPlusPlus] = lowEnergyLim
186 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
187
188 // Cross section
189
190 auto tableAlphaPlusPlus = new G4DNACrossSec
191 eV,
192 scaleFactor );
193 tableAlphaPlusPlus->LoadData(fileAlphaPlusPl
194
195 tableData[alphaPlusPlus] = tableAlphaPlusPlu
196
197 // *****************************************
198
199 alphaPlus = alphaPlusDef->GetParticleName();
200 tableFile[alphaPlus] = fileAlphaPlus;
201
202 lowEnergyLimit[alphaPlus] = lowEnergyLimitFo
203 highEnergyLimit[alphaPlus] = 400. * MeV;
204
205 // Cross section
206
207 auto tableAlphaPlus = new G4DNACrossSection
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 G4DNACrossSectionDat
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[hydroge
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[alphaPl
253 }
254
255 if (particle==alphaPlusPlusDef)
256 {
257 SetLowEnergyLimit(lowEnergyLimit[alphaPlus
258 SetHighEnergyLimit(highEnergyLimit[alphaPl
259 }
260
261 if( verboseLevel>0 )
262 {
263 G4cout << "Rudd ionisation model is initia
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::Ins
273
274 //
275
276 fAtomDeexcitation = G4LossTableManager::Inst
277
278 if (isInitialised)
279 { return;}
280 fParticleChangeForGamma = GetParticleChangeF
281 isInitialised = true;
282
283 }
284
285 //....oooOO0OOooo........oooOO0OOooo........oo
286
287 G4double G4DNARuddIonisationModel::CrossSectio
288
289
290
291
292 {
293 if (verboseLevel > 3)
294 {
295 G4cout << "Calling CrossSectionPerVolume()
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)[ma
334
335 const G4String& particleName = particleDefin
336
337 // SI - the following is useless since lowLi
338 /*
339 std::map< G4String,G4double,std::less<G4Stri
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<G4Stri
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 sampl
359
360 if (k < lowLim) k = lowLim;
361
362 //
363
364 std::map< G4String,G4DNACrossSectionDataSe
365 pos = tableData.find(particleName);
366
367 if (pos != tableData.end())
368 {
369 G4DNACrossSectionDataSet* table = pos->s
370 if (table != nullptr)
371 {
372 sigma = table->FindValue(k);
373 }
374 }
375 else
376 {
377 G4Exception("G4DNARuddIonisationModel::C
378 FatalException,"Model not applicable
379 }
380
381 }
382
383 if (verboseLevel > 2)
384 {
385 G4cout << "_______________________________
386 G4cout << "G4DNARuddIonisationModel - XS I
387 G4cout << "Kinetic energy(eV)=" << k/eV <<
388 G4cout << "Cross section per water molecul
389 G4cout << "Cross section per water molecul
390 //G4cout << " - Cross section per water mo
391 //<< sigma*material->GetAtomicNumDensityVe
392 G4cout << "G4DNARuddIonisationModel - XS I
393 }
394
395 return sigma*waterDensity;
396
397 }
398
399 //....oooOO0OOooo........oooOO0OOooo........oo
400
401 void G4DNARuddIonisationModel::SampleSecondari
402
403
404
405
406 {
407 if (verboseLevel > 3)
408 {
409 G4cout << "Calling SampleSecondaries() of
410 << G4endl;
411 }
412
413 G4double lowLim = 0;
414 G4double highLim = 0;
415
416 if ( particle->GetDefinition() == protonDef
417 || particle->GetDefinition() == hydrogen
418 )
419
420 lowLim = killBelowEnergyForZ1;
421
422 if ( particle->GetDefinition() == alphaPlusP
423 || particle->GetDefinition() == alphaPlu
424 || particle->GetDefinition() == heliumDe
425 )
426
427 lowLim = killBelowEnergyForZ2;
428
429 G4double k = particle->GetKineticEnergy();
430
431 const G4String& particleName = particle->Get
432
433 // SI - the following is useless since lowLi
434 /*
435 std::map< G4String,G4double,std::less<G4Str
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<G4Stri
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 = particl
455 G4ParticleMomentum primaryDirection = part
456 /*
457 G4double particleMass = definition->GetPD
458 G4double totalEnergy = k + particleMass;
459 G4double pSquare = k*(totalEnergy+particl
460 G4double totalMomentum = std::sqrt(pSquar
461 */
462
463 G4int ionizationShell = RandomSelect(k,par
464
465 G4double bindingEnergy = 0;
466 bindingEnergy = waterStructure.IonisationE
467
468 //SI: additional protection if tcs interpo
469 if (k<bindingEnergy) return;
470 //
471
472 // SI - For atom. deexc. tagging - 23/05/2
473 G4int Z = 8;
474 //
475
476 G4double secondaryKinetic = RandomizeEject
477
478 G4ThreeVector deltaDirection =
479 GetAngularDistribution()->SampleDirectionF
480 Z, ionizationShell,
481 couple->GetMaterial());
482
483 auto dp = new G4DynamicParticle (G4Electr
484 fvect->push_back(dp);
485
486 // Ignored for ions on electrons
487 /*
488 G4double deltaTotalMomentum = std::sqrt(s
489
490 G4double finalPx = totalMomentum*primaryD
491 G4double finalPy = totalMomentum*primaryD
492 G4double finalPz = totalMomentum*primaryD
493 G4double finalMomentum = std::sqrt(finalP
494 finalPx /= finalMomentum;
495 finalPy /= finalMomentum;
496 finalPz /= finalMomentum;
497
498 G4ThreeVector direction;
499 direction.set(finalPx,finalPy,finalPz);
500
501 fParticleChangeForGamma->ProposeMomentumD
502 */
503
504 fParticleChangeForGamma->ProposeMomentumDi
505
506 // sample deexcitation
507 // here we assume that H_{2}O electronic l
508 // this can be considered true with a roug
509
510 size_t secNumberInit = 0;// need to know a
511 size_t secNumberFinal = 0;// So I'll make
512
513 G4double scatteredEnergy = k-bindingEnergy
514
515 // SI: only atomic deexcitation from K she
516 if((fAtomDeexcitation != nullptr) && ioniz
517 {
518 const G4AtomicShell* shell
519 = fAtomDeexcitation->GetAtomicShell(Z,
520 secNumberInit = fvect->size();
521 fAtomDeexcitation->GenerateParticles(fve
522 secNumberFinal = fvect->size();
523
524 if(secNumberFinal > secNumberInit)
525 {
526 for (size_t i=secNumberInit; i<secNumberFina
527 {
528 //Check if there is enough residual
529 if (bindingEnergy >= ((*fvect)[i])->
530 {
531 //Ok, this is a valid secondary:
532 bindingEnergy -= ((*fvect)[i])->GetKine
533 }
534 else
535 {
536 //Invalid secondary: not enough energy
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("G4DNAEmfietzoglouIonisatioMo
549 "em2050",FatalException,"Nega
550
551 //bindingEnergy has been decreased
552 //by the amount of energy taken away by de
553 if (!statCode)
554 {
555 fParticleChangeForGamma->SetProposedKine
556 fParticleChangeForGamma->ProposeLocalEne
557 }
558 else
559 {
560 fParticleChangeForGamma->SetProposedKine
561 fParticleChangeForGamma->ProposeLocalEne
562 }
563
564 // debug
565 // k-scatteredEnergy-secondaryKinetic-deex
566 // = k-k+bindingEnergy+secondaryKinetic-se
567 // = bindingEnergy-deexSecEnergy
568 // SO deexSecEnergy=0 => LocalEnergyDeposi
569
570 // TEST //////////////////////////
571 // if (secondaryKinetic<0) abort();
572 // if (scatteredEnergy<0) abort();
573 // if (k-scatteredEnergy-secondaryKinetic-
574 // if (k-scatteredEnergy<0) abort();
575 /////////////////////////////////
576
577 const G4Track * theIncomingTrack = fPartic
578 G4DNAChemistryManager::Instance()->CreateW
579 ionizationShell,
580 theIncomingTrack);
581 }
582
583 // SI - not useful since low energy of model
584
585 if (k < lowLim)
586 {
587 fParticleChangeForGamma->SetProposedKineti
588 fParticleChangeForGamma->ProposeTrackStatu
589 fParticleChangeForGamma->ProposeLocalEnerg
590 }
591
592 }
593
594 //....oooOO0OOooo........oooOO0OOooo........oo
595
596 G4double G4DNARuddIonisationModel::RandomizeEj
597
598
599 {
600 G4double maximumKineticEnergyTransfer = 0.;
601
602 if (particleDefinition == protonDef
603 || particleDefinition == hydrogenDef)
604 {
605 maximumKineticEnergyTransfer = 4. * (elect
606 }
607
608 else if (particleDefinition == heliumDef
609 || particleDefinition == alphaPlusDef
610 || particleDefinition == alphaPlusPlusDe
611 {
612 maximumKineticEnergyTransfer = 4. * (0.511
613 }
614
615 G4double crossSectionMaximum = 0.;
616
617 for (G4double value = waterStructure.Ionisat
618 value <= 5. * waterStructure.IonisationE
619 value += 0.1 * eV)
620 {
621 G4double differentialCrossSection =
622 DifferentialCrossSection(particleDefin
623 if (differentialCrossSection >= crossSecti
624 crossSectionMaximum = differentialCrossS
625 }
626
627 G4double secElecKinetic = 0.;
628
629 do
630 {
631 secElecKinetic = G4UniformRand()* maximumK
632 }while(G4UniformRand()*crossSectionMaximum >
633 k,
634 secElecKinetic+waterStructure.Ionisa
635 shell));
636
637 return (secElecKinetic);
638 }
639
640 //....oooOO0OOooo........oooOO0OOooo........oo
641
642 // The following section is not used anymore b
643 // GetAngularDistribution()->SampleDirectionFo
644
645 /*
646 void G4DNARuddIonisationModel::RandomizeEject
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 / proto
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 / maxSecKi
673
674 // Restriction below 100 eV from Emfietzoglou
675
676 if (secKinetic>100*eV) cosTheta = std::sqrt(s
677 else cosTheta = (2.*G4UniformRand())-1.;
678
679 }
680 */
681 //....oooOO0OOooo........oooOO0OOooo........oo
682 G4double G4DNARuddIonisationModel::Differentia
683
684
685
686 {
687 // Shells ids are 0 1 2 3 4 (4 is k shell)
688 // !!Attention, "energyTransfer" here is the
689 // that the secondary kinetic en
690 //
691 // ds S F1(nu) +
692 // ---- = G(k) * ---- -----------------
693 // dw Bj (1+w)^3 * [1 + e
694 //
695 // w is the secondary electron kinetic Energ
696 //
697 // All the other parameters can be found in
698 //
699 // M.Eugene Rudd, 1988, User-Friendly model
700 // electrons from protons or electron collis
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*
717 // The following values are provided by M. d
718 const G4double Bj[5] = { 12.60 * eV, 14.70 *
719 * eV };
720
721 if (j == 4)
722 {
723 //Data For Liquid Water K SHELL from Dingf
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 (P
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.
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
753
754 G4double wBig = (energyTransfer
755 - waterStructure.IonisationEnergy(ioniza
756 if (wBig < 0)
757 return 0.;
758
759 G4double w = wBig / Bj[ionizationLevelIndex]
760 // Note that the following (j==4) cases are
761 if (j == 4)
762 w = wBig / waterStructure.IonisationEnergy
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)
776 }
777
778 else if (particleDefinition == heliumDef
779 || particleDefinition == alphaPlusDef
780 || particleDefinition == alphaPlusPlusDe
781 {
782 isHelium = true;
783 tau = (0.511 / 3728.) * k;
784 }
785
786 G4double S = 4. * pi * Bohr_radius * Bohr_ra
787 * gpow->powN((Ry / Bj[ionizationLevelInd
788 if (j == 4)
789 S = 4. * pi * Bohr_radius * Bohr_radius *
790 * gpow->powN((Ry / waterStructure.Ioni
791 2);
792
793 G4double v2 = tau / Bj[ionizationLevelIndex]
794 if (j == 4)
795 v2 = tau / waterStructure.IonisationEnergy
796
797 G4double v = std::sqrt(v2);
798 G4double wc = 4. * v2 - 2. * v - (Ry / (4. *
799 if (j == 4)
800 wc = 4. * v2 - 2. * v
801 - (Ry / (4. * waterStructure.Ionisatio
802
803 G4double L1 = (C1 * gpow->powA(v, (D1))) / (
804 G4double L2 = C2 * gpow->powA(v, (D2));
805 G4double H1 = (A1 * G4Log(1. + v2)) / (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)
813 * (S / Bj[ionizationLevelIndex])
814 * ((F1 + w * F2)
815 / (gpow->powN((1. + w), 3)
816 * (1. + G4Exp(alphaConst * (
817
818 if (j == 4)
819 sigma = CorrectionFactor(particleDefinitio
820 * (S / waterStructure.IonisationEnergy
821 * ((F1 + w * F2)
822 / (gpow->powN((1. + w), 3)
823 * (1. + G4Exp(alphaConst * (w
824
825 if ((particleDefinition == hydrogenDef)
826 && (ionizationLevelIndex == 4))
827 {
828 // sigma = Gj[j] * (S/Bj[ionizationLeve
829 sigma = Gj[j] * (S / waterStructure.Ionisa
830 * ((F1 + w * F2)
831 / (gpow->powN((1. + w), 3)
832 * (1. + G4Exp(alphaConst * (w
833 }
834
835 // if ( particleDefinition == protonDe
836 // || particleDefinition == hydro
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.
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 == heliumDe
877 // || particleDefinition == alpha
878 // || particleDefinition == alpha
879 // )
880 if (isHelium)
881 {
882 sigma = Gj[j] * (S / Bj[ionizationLevelInd
883 * ((F1 + w * F2)
884 / (gpow->powN((1. + w), 3)
885 * (1. + G4Exp(alphaConst * (w
886
887 if (j == 4)
888 sigma = Gj[j]
889 * (S / waterStructure.IonisationEner
890 * ((F1 + w * F2)
891 / (gpow->powN((1. + w), 3)
892 * (1. + G4Exp(alphaConst * (
893
894 G4double zEff = particleDefinition->GetPDG
895 + particleDefinition->GetLeptonNumber(
896
897 zEff -= (sCoefficient[0]
898 * S_1s(k, energyTransfer, slaterEffect
899 + sCoefficient[1]
900 * S_2s(k, energyTransfer, slaterEf
901 + sCoefficient[2]
902 * S_2p(k, energyTransfer, slaterEf
903
904 return zEff * zEff * sigma;
905 }
906
907 return 0;
908 }
909
910 //....oooOO0OOooo........oooOO0OOooo........oo
911
912 G4double G4DNARuddIonisationModel::S_1s(G4doub
913 G4doub
914 G4doub
915 G4doub
916 {
917 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
918 // Dingfelder, in Chattanooga 2005 proceedin
919
920 G4double r = R(t, energyTransferred, slaterE
921 G4double value = 1. - G4Exp(-2 * r) * ((2. *
922
923 return value;
924 }
925
926 //....oooOO0OOooo........oooOO0OOooo........oo
927
928 G4double G4DNARuddIonisationModel::S_2s(G4doub
929 G4doub
930 G4doub
931 G4doub
932 {
933 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
934 // Dingfelder, in Chattanooga 2005 proceedin
935
936 G4double r = R(t, energyTransferred, slaterE
937 G4double value = 1.
938 - G4Exp(-2 * r) * (((2. * r * r + 2.) *
939
940 return value;
941
942 }
943
944 //....oooOO0OOooo........oooOO0OOooo........oo
945
946 G4double G4DNARuddIonisationModel::S_2p(G4doub
947 G4doub
948 G4doub
949 G4doub
950 {
951 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^
952 // Dingfelder, in Chattanooga 2005 proceedin
953
954 G4double r = R(t, energyTransferred, slaterE
955 G4double value = 1.
956 - G4Exp(-2 * r)
957 * ((((2. / 3. * r + 4. / 3.) * r + 2
958
959 return value;
960 }
961
962 //....oooOO0OOooo........oooOO0OOooo........oo
963
964 G4double G4DNARuddIonisationModel::R(G4double
965 G4double
966 G4double
967 G4double
968 {
969 // tElectron = m_electron / m_alpha * t
970 // Dingfelder, in Chattanooga 2005 proceedin
971
972 G4double tElectron = 0.511 / 3728. * t;
973 // The following values are provided by M. D
974 G4double H = 2. * 13.60569172 * eV;
975 G4double value = std::sqrt(2. * tElectron /
976 * (slaterEffectiveChg / shellNumber);
977
978 return value;
979 }
980
981 //....oooOO0OOooo........oooOO0OOooo........oo
982
983 G4double G4DNARuddIonisationModel::CorrectionF
984
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
994 // The following values are provided by M.
995 return ((0.6 / (1 + G4Exp(value))) + 0.9);
996 }
997 return (1.);
998 }
999
1000 //....oooOO0OOooo........oooOO0OOooo........o
1001
1002 G4int G4DNARuddIonisationModel::RandomSelect(
1003
1004 {
1005
1006 // BEGIN PART 1/2 OF ELECTRON CORRECTION
1007
1008 // add ONE or TWO electron-water ionisation
1009
1010 G4int level = 0;
1011
1012 // Retrieve data table corresponding to the
1013
1014 std::map<G4String, G4DNACrossSectionDataSet
1015 pos = tableData.find(particle);
1016
1017 if (pos != tableData.end())
1018 {
1019 G4DNACrossSectionDataSet* table = pos->se
1020
1021 if (table != nullptr)
1022 {
1023 auto valuesBuffer = new G4double[table
1024
1025 const auto n = (G4int)table->NumberOfC
1026 G4int i(n);
1027 G4double value = 0.;
1028
1029 while (i > 0)
1030 {
1031 --i;
1032 valuesBuffer[i] = table->GetComponent
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::Ra
1059 "em0002",
1060 FatalException,
1061 "Model not applicable to part
1062 }
1063
1064 return level;
1065 }
1066
1067 //....oooOO0OOooo........oooOO0OOooo........o
1068
1069 G4double G4DNARuddIonisationModel::PartialCro
1070 {
1071 G4double sigma = 0.;
1072
1073 const G4DynamicParticle* particle = track.G
1074 G4double k = particle->GetKineticEnergy();
1075
1076 G4double lowLim = 0;
1077 G4double highLim = 0;
1078
1079 const G4String& particleName = particle->Ge
1080
1081 std::map<G4String, G4double, std::less<G4St
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<G4St
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, G4DNACrossSectionDataS
1100 pos = tableData.find(particleName);
1101
1102 if (pos != tableData.end())
1103 {
1104 G4DNACrossSectionDataSet* table = pos->
1105 if (table != nullptr)
1106 {
1107 sigma = table->FindValue(k);
1108 }
1109 } else
1110 {
1111 G4Exception("G4DNARuddIonisationModel::
1112 "em0002",
1113 FatalException,
1114 "Model not applicable to pa
1115 }
1116 }
1117
1118 return sigma;
1119 }
1120
1121 //....oooOO0OOooo........oooOO0OOooo........o
1122
1123 G4double G4DNARuddIonisationModel::Sum(G4doub
1124 const
1125 {
1126 return 0;
1127 }
1128
1129