| 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 // Reference:
27 // A.D. Dominguez-Munoz, M.I. Gallardo, M.C
28 // Z. Francis, S. Incerti, M.A. Cortes-Gira
29 // Radiat. Phys. Chem. 199 (2022) 110363.
30 //
31 // Class authors:
32 // A.D. Dominguez-Munoz
33 // M.A. Cortes-Giraldo (miancortes -at- us.
34 //
35 // Class creation: 2022-03-03
36 //
37 //
38
39 #include "G4DNARPWBAIonisationModel.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4UAtomicDeexcitation.hh"
43 #include "G4LossTableManager.hh"
44 #include "G4DNAChemistryManager.hh"
45 #include "G4DNAMolecularMaterial.hh"
46 #include "G4DNABornAngle.hh"
47 #include "G4Exp.hh"
48 using namespace std;
49 //....oooOO0OOooo........oooOO0OOooo........oo
50 G4DNARPWBAIonisationModel::G4DNARPWBAIonisatio
51 const G4ParticleDefinition*, const G4String&
52 : G4VEmModel(nam)
53 {
54 // Verbosity scale:
55 // 0 = nothing
56 // 1 = warning for energy non-conservation
57 // 2 = details of energy budget
58 // 3 = calculation of cross sections, file o
59 // 4 = entering in methods
60
61 if(verboseLevel > 0)
62 {
63 G4cout << "RPWBA ionisation model is const
64 }
65 SetDeexcitationFlag(true);
66 SetAngularDistribution(new G4DNABornAngle())
67 }
68
69 //....oooOO0OOooo........oooOO0OOooo........oo
70
71 G4DNARPWBAIonisationModel::~G4DNARPWBAIonisati
72 {
73 eVecm.clear();
74 pVecm.clear();
75 }
76 //....oooOO0OOooo........oooOO0OOooo........oo
77
78 G4bool G4DNARPWBAIonisationModel::InEnergyLimi
79 {
80 if(lowEnergyLimit == highEnergyLimit)
81 {
82 G4Exception("G4DNARPWBAIonisationModel::In
83 FatalException, "lowEnergyLimi
84 }
85 return k >= lowEnergyLimit && k <= highEnerg
86 }
87
88 //....oooOO0OOooo........oooOO0OOooo........oo
89 void G4DNARPWBAIonisationModel::InitialiseForP
90 const G4ParticleDefinition* part)
91 {
92 if(part != fProtonDef)
93 {
94 G4Exception("G4DNARPWBAIonisationModel::In
95 FatalException, "Model not app
96 }
97 // Energy limits
98 G4String fileProton("dna/sigma_ionisation_p_
99 G4double scaleFactor = 1 * cm * cm;
100 const char *path = G4FindDataDir("G4LEDAT
101 lowEnergyLimit = 100. * MeV;
102 highEnergyLimit = 300. * MeV;
103
104 if(LowEnergyLimit() < lowEnergyLimit || High
105 {
106 G4ExceptionDescription ed;
107 ed << "Model is applicable from "<<lowEner
108 G4Exception("G4DNARPWBAIonisationModel::In
109 FatalException, ed);
110 }
111
112 fpTotalCrossSection = make_unique<G4DNACross
113 new G4LogLogInterpolation, eV, scaleFactor
114 fpTotalCrossSection->LoadData(fileProton);
115
116 // Final state
117
118 std::ostringstream pFullFileName;
119 fasterCode ? pFullFileName
120 << path << "/dna/sigmadiff_cu
121 : pFullFileName << path << "/dna/
122 std::ifstream pDiffCrossSection(pFullFileNam
123 if(!pDiffCrossSection)
124 {
125 G4ExceptionDescription exceptionDescriptio
126 exceptionDescription << "Missing data file
127 G4Exception("G4DNARPWBAIonisationModel::In
128 FatalException, exceptionDescr
129 }
130
131 pTdummyVec.push_back(0.);
132 while(!pDiffCrossSection.eof())
133 {
134 G4double tDummy;
135 G4double eDummy;
136 pDiffCrossSection >> tDummy >> eDummy;
137 if(tDummy != pTdummyVec.back())
138 {
139 pTdummyVec.push_back(tDummy);
140 }
141
142 for(G4int j = 0; j < 5; j++)
143 {
144 pDiffCrossSection >> pDiffCrossSectionDa
145
146 if(fasterCode)
147 {
148 pNrjTransfData[j][tDummy][pDiffCrossSe
149 eDummy;
150 pProbaShellMap[j][tDummy].push_back(
151 pDiffCrossSectionData[j][tDummy][eDu
152 }
153
154 // SI - only if eof is not reached !
155 if(!pDiffCrossSection.eof() && !fasterCo
156 {
157 pDiffCrossSectionData[j][tDummy][eDumm
158 }
159
160 if(!fasterCode)
161 {
162 pVecm[tDummy].push_back(eDummy);
163 }
164 }
165 }
166
167 // be careful about this
168 // SetLowEnergyLimit(lowEnergyLimit);
169 // SetHighEnergyLimit(highEnergyLimit);
170 }
171
172 void G4DNARPWBAIonisationModel::Initialise(con
173 con
174 {
175 if(isInitialised)
176 {
177 return;
178 }
179 if(verboseLevel > 3)
180 {
181 G4cout << "Calling G4DNARPWBAIonisationMod
182 << particle->GetParticleName() << G
183 }
184
185 InitialiseForProton(particle);
186
187 if(verboseLevel > 0)
188 {
189 G4cout << "RPWBA ionisation model is initi
190 << "Energy range: " << LowEnergyLim
191 << HighEnergyLimit() / MeV << " MeV
192 << particle->GetParticleName() << G
193 }
194
195 // Initialize water density pointer
196 if(G4Material::GetMaterial("G4_WATER") != nu
197 {
198 fpMolWaterDensity =
199 G4DNAMolecularMaterial::Instance()->GetN
200 G4Material::GetMaterial("G4_WATER"));
201 }
202 else
203 {
204 G4ExceptionDescription exceptionDescriptio
205 exceptionDescription << "G4_WATER does not
206 G4Exception("G4DNARPWBAIonisationModel::In
207 FatalException, exceptionDescr
208 }
209 fAtomDeexcitation = G4LossTableManager
210 fParticleChangeForGamma = GetParticleChangeF
211 isInitialised = true;
212 }
213
214 //....oooOO0OOooo........oooOO0OOooo........oo
215
216 G4double G4DNARPWBAIonisationModel::CrossSecti
217 const G4Material* material, const G4Particle
218 G4double ekin, G4double, G4double)
219 {
220 if(particleDefinition != fProtonDef)
221 {
222 G4Exception("G4DNARPWBAIonisationModel::Cr
223 FatalException, "Model not app
224 }
225 if(verboseLevel > 3)
226 {
227 G4cout << "Calling CrossSectionPerVolume()
228 << G4endl;
229 }
230 G4double sigma;
231 G4double waterDensity = (*fpMolWaterDensity)
232
233 if(InEnergyLimit(ekin))
234 {
235 sigma = fpTotalCrossSection->FindValue(eki
236 }
237 else
238 {
239 // nput energy is outside this interval th
240 // should add a warning or exception ?
241 return 0;
242 }
243
244 if(verboseLevel > 2)
245 {
246 G4cout << "_______________________________
247 G4cout << "G4DNARPWBAIonisationModel - XS
248 G4cout << "Kinetic energy(eV)=" << ekin /
249 << " particle : " << fProtonDef->Ge
250 G4cout << "Cross section per water molecul
251 << G4endl;
252 G4cout << "Cross section per water molecul
253 << sigma * waterDensity / (1. / cm)
254 G4cout << "G4DNARPWBAIonisationModel - XS
255 }
256
257 return sigma * waterDensity;
258 }
259
260 //....oooOO0OOooo........oooOO0OOooo........oo
261
262 void G4DNARPWBAIonisationModel::SampleSecondar
263 std::vector<G4DynamicParticle*>* fvect, cons
264 const G4DynamicParticle* particle, G4double,
265 {
266 if(verboseLevel > 3)
267 {
268 G4cout << "Calling SampleSecondaries() of
269 << G4endl;
270 }
271 G4double k = particle->GetKineticEnergy();
272 if(InEnergyLimit(k))
273 {
274 G4ParticleMomentum primaryDirection = part
275 G4double particleMass = particle->GetDefi
276 G4double totalEnergy = k + particleMass;
277 G4double pSquare = k * (totalEnergy
278 G4double totalMomentum = std::sqrt(pSquare
279 G4int ionizationShell;
280
281 if(!fasterCode)
282 {
283 ionizationShell = RandomSelect(k);
284 }
285 else
286 {
287 // fasterCode = true
288 do
289 {
290 ionizationShell = RandomSelect(k);
291 } while(k < 19 * eV && ionizationShell =
292 particle->GetDefinition() == G4E
293 }
294
295 G4double bindingEnergy = waterStructure.Io
296
297 // SI: additional protection if tcs interp
298 if(k < bindingEnergy)
299 {
300 return;
301 }
302 //
303 G4double secondaryKinetic;
304 if(!fasterCode)
305 {
306 secondaryKinetic = RandomizeEjectedElect
307 }
308 else
309 {
310 secondaryKinetic =
311 RandomizeEjectedElectronEnergyFromCumu
312 }
313
314 G4int Z = 8; // water Z (6 Oxygen + 2 hyd
315 G4ThreeVector deltaDirection =
316 GetAngularDistribution()->SampleDirectio
317 particle, secondaryKinetic, Z, ionizat
318
319 if(secondaryKinetic > 0){
320 auto dp = new G4DynamicParticle(G4Electr
321 secondar
322 fvect->push_back(dp);
323 }
324
325 if(particle->GetDefinition() == G4Electron
326 G4double deltaTotalMomentum = std::sqrt(
327 secondaryKinetic * (secondaryKinetic +
328
329 G4double finalPx = totalMomentum * prima
330 deltaTotalMomentum *
331 G4double finalPy = totalMomentum * prima
332 deltaTotalMomentum *
333 G4double finalPz = totalMomentum * prima
334 deltaTotalMomentum *
335 G4double finalMomentum =
336 std::sqrt(finalPx * finalPx + finalPy
337 finalPx /= finalMomentum;
338 finalPy /= finalMomentum;
339 finalPz /= finalMomentum;
340 G4ThreeVector direction;
341 direction.set(finalPx, finalPy, finalPz)
342 fParticleChangeForGamma->ProposeMomentum
343 }
344 else
345 {
346 fParticleChangeForGamma->ProposeMomentum
347 }
348
349 // AM: sample deexcitation
350 // here we assume that H_{2}O electronic l
351 // this can be considered true with a roug
352
353 size_t secNumberInit; // need to know at
354 // secondaries
355 size_t
356 secNumberFinal; // So I'll make the dif
357
358 G4double scatteredEnergy = k - bindingEner
359
360 // SI: only atomic deexcitation from K she
361 if((fAtomDeexcitation != nullptr) && ioniz
362 {
363 const G4AtomicShell* shell =
364 fAtomDeexcitation->GetAtomicShell(Z, G
365 secNumberInit = fvect->size();
366 fAtomDeexcitation->GenerateParticles(fve
367 secNumberFinal = fvect->size();
368
369 if(secNumberFinal > secNumberInit){
370 for(size_t i = secNumberInit; i < secN
371 if(bindingEnergy >= ((*fvect)[i])->G
372 {
373 bindingEnergy -= ((*fvect)[i])->Ge
374 }else{
375 delete(*fvect)[i];
376 (*fvect)[i] = nullptr;
377 }
378 }
379 }
380 }
381
382 // This should never happen
383 if(bindingEnergy < 0.0)
384 {
385 G4Exception("G4DNARPWBAIonisatioModel::S
386 FatalException, "Negative lo
387 }
388
389 // bindingEnergy has been decreased
390 // by the amount of energy taken away by d
391 if(!statCode){
392 fParticleChangeForGamma->SetProposedKine
393 fParticleChangeForGamma->ProposeLocalEne
394 }else{
395 fParticleChangeForGamma->SetProposedKine
396 fParticleChangeForGamma->ProposeLocalEne
397 }
398 const G4Track* theIncomingTrack =
399 fParticleChangeForGamma->GetCurrentTrack
400 G4DNAChemistryManager::Instance()->CreateW
401 eIonizedMolecule, ionizationShell, theIn
402 }
403 }
404
405 //....oooOO0OOooo........oooOO0OOooo........oo
406
407 G4double G4DNARPWBAIonisationModel::RandomizeE
408 const G4double& k, const G4int& shell)
409 {
410 G4double maximumKineticEnergyTransfer =
411 4. * (electron_mass_c2 / proton_mass_c2) *
412 G4double IonisationEnergyInShell = waterStru
413 G4double kIneV = k / eV;
414
415 G4double crossSectionMaximum = 0.;
416 for(G4double value = IonisationEnergyInShell
417 value <= 4. * IonisationEnergyInShell; v
418 {
419 G4double differentialCrossSection =
420 DifferentialCrossSection(kIneV, value /
421 if(differentialCrossSection >= crossSectio
422 {
423 crossSectionMaximum = differentialCrossS
424 }
425 }
426
427 G4double secondaryElectronKineticEnergy;
428 do
429 {
430 secondaryElectronKineticEnergy =
431 G4UniformRand() * maximumKineticEnergyTr
432 } while(G4UniformRand() * crossSectionMaximu
433 DifferentialCrossSection(kIneV,
434 (secondaryElectronKineticEnergy +
435 IonisationEnergyInShell) / eV,
436
437 return secondaryElectronKineticEnergy;
438 }
439
440 //....oooOO0OOooo........oooOO0OOooo........oo
441 G4double G4DNARPWBAIonisationModel::Differenti
442 const G4double& kine, const G4double& energy
443 const G4int& ionizationLevelIndex)
444 {
445 G4double k = kine;
446 G4double sigma = 0.;
447 if(energyTransfer >=
448 waterStructure.IonisationEnergy(ionizatio
449 {
450 G4double valueT1 = 0;
451 G4double valueT2 = 0;
452 G4double valueE21 = 0;
453 G4double valueE22 = 0;
454 G4double valueE12 = 0;
455 G4double valueE11 = 0;
456
457 G4double xs11 = 0;
458 G4double xs12 = 0;
459 G4double xs21 = 0;
460 G4double xs22 = 0;
461
462 // Protection against out of boundary acce
463
464 if(k == pTdummyVec.back())
465 {
466 k = k * (1. - 1e-12);
467 }
468 // k should be in eV and energy transfer e
469 auto t2 = std::upper_bound(pTdummyVec.begi
470 auto t1 = t2 - 1;
471
472 auto e12 = std::upper_bound(pVecm[(*t1)].b
473 energyTransfer
474 auto e11 = e12 - 1;
475
476 auto e22 = std::upper_bound(pVecm[(*t2)].b
477 energyTransfer
478 auto e21 = e22 - 1;
479
480 valueT1 = *t1;
481 valueT2 = *t2;
482 valueE21 = *e21;
483 valueE22 = *e22;
484 valueE12 = *e12;
485 valueE11 = *e11;
486
487 xs11 = pDiffCrossSectionData[ionizationLev
488 xs12 = pDiffCrossSectionData[ionizationLev
489 xs21 = pDiffCrossSectionData[ionizationLev
490 xs22 = pDiffCrossSectionData[ionizationLev
491
492 G4double xsProduct = xs11 * xs12 * xs21 *
493 if(xsProduct != 0.)
494 {
495 sigma =
496 QuadInterpolator(valueE11, valueE12, v
497 xs21, xs22, valueT1,
498 }
499 }
500 return sigma;
501 }
502
503 //....oooOO0OOooo........oooOO0OOooo........oo
504
505 G4double G4DNARPWBAIonisationModel::Interpolat
506
507
508
509
510 {
511 G4double value = 0.;
512
513 // Log-log interpolation by default
514
515 if(e1 != 0 && e2 != 0 && (std::log10(e2) - s
516 !fasterCode)
517 {
518 G4double a =
519 (std::log10(xs2) - std::log10(xs1)) / (s
520 G4double b = std::log10(xs2) - a * std
521 G4double sigma = a * std::log10(e) + b;
522 value = (std::pow(10., sigma));
523 }
524 // Switch to lin-lin interpolation
525 /*
526 if ((e2-e1)!=0)
527 {
528 G4double d1 = xs1;
529 G4double d2 = xs2;
530 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1)
531 }
532 */
533
534 // Switch to log-lin interpolation for faste
535 if((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 &&
536 {
537 G4double d1 = std::log10(xs1);
538 G4double d2 = std::log10(xs2);
539 value = std::pow(10., (d1 + (d2 - d1
540 }
541
542 // Switch to lin-lin interpolation for faste
543 // in case one of xs1 or xs2 (=cum proba) va
544
545 if((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0)
546 {
547 G4double d1 = xs1;
548 G4double d2 = xs2;
549 value = (d1 + (d2 - d1) * (e - e1) /
550 }
551 return value;
552 }
553
554 //....oooOO0OOooo........oooOO0OOooo........oo
555
556 G4double G4DNARPWBAIonisationModel::QuadInterp
557 const G4double& e11, const G4double& e12, co
558 const G4double& e22, const G4double& xs11, c
559 const G4double& xs21, const G4double& xs22,
560 const G4double& t2, const G4double& t, const
561 {
562 G4double interpolatedvalue1 = Interpolate(e1
563 G4double interpolatedvalue2 = Interpolate(e2
564 G4double value =
565 Interpolate(t1, t2, t, interpolatedvalue1,
566
567 return value;
568 }
569
570 G4double G4DNARPWBAIonisationModel::GetPartial
571 const G4Material* /*material*/, G4int level,
572 const G4ParticleDefinition* particle, G4doub
573 {
574 if(fpTotalCrossSection != nullptr && particl
575 {
576 G4Exception("G4DNARPWBAIonisationModel::Ge
577 FatalException, "Model not app
578 }
579 return fpTotalCrossSection->GetComponent(lev
580 }
581
582 //....oooOO0OOooo........oooOO0OOooo........oo
583
584 G4int G4DNARPWBAIonisationModel::RandomSelect(
585 {
586 if(fpTotalCrossSection == nullptr)
587 {
588 G4Exception("G4DNARPWBAIonisationModel::Ra
589 FatalException, "Model not app
590 }
591 else
592 {
593 auto valuesBuffer = new G4double[fpTotalCr
594 const auto n = (G4int)fpTotalCrossSection
595 G4int i(n);
596 G4double value = 0.;
597 while(i > 0)
598 {
599 --i;
600 valuesBuffer[i] = fpTotalCrossSection->G
601 value += valuesBuffer[i];
602 }
603 value *= G4UniformRand();
604 i = n;
605
606 while(i > 0)
607 {
608 --i;
609 if(valuesBuffer[i] > value)
610 {
611 delete[] valuesBuffer;
612 return i;
613 }
614 value -= valuesBuffer[i];
615 }
616 delete[] valuesBuffer;
617 }
618 return 0;
619 }
620
621 //....oooOO0OOooo........oooOO0OOooo........oo
622
623 G4double
624 G4DNARPWBAIonisationModel::RandomizeEjectedEle
625 const G4double& k, const G4int& shell)
626 {
627 G4double random = G4
628 G4double secondaryKineticEnergy =
629 TransferedEnergy(k / eV, shell, random) *
630 waterStructure.IonisationEnergy(shell);
631 if(secondaryKineticEnergy < 0.)
632 {
633 return 0.;
634 }
635 return secondaryKineticEnergy;
636 }
637
638 //....oooOO0OOooo........oooOO0OOooo........oo
639
640 G4double G4DNARPWBAIonisationModel::Transfered
641
642
643 {
644 G4double nrj = 0.;
645 G4double valueK1 = 0;
646 G4double valueK2 = 0;
647 G4double valuePROB21 = 0;
648 G4double valuePROB22 = 0;
649 G4double valuePROB12 = 0;
650 G4double valuePROB11 = 0;
651 G4double nrjTransf11 = 0;
652 G4double nrjTransf12 = 0;
653 G4double nrjTransf21 = 0;
654 G4double nrjTransf22 = 0;
655 // Protection against out of boundary access
656 if(k == pTdummyVec.back())
657 {
658 k = k * (1. - 1e-12);
659 }
660 // k should be in eV
661
662 auto k2 = std::upper_bound(pTdummyVec.begin(
663 auto k1 = k2 - 1;
664
665 // SI : the following condition avoids situa
666 // element,
667 // for eg. when the last element is zer
668 if(random <= pProbaShellMap[ionizationLevelI
669 random <= pProbaShellMap[ionizationLevelI
670 {
671 auto prob12 = std::upper_bound(
672 pProbaShellMap[ionizationLevelIndex][(*k
673 pProbaShellMap[ionizationLevelIndex][(*k
674 auto prob11 = prob12 - 1;
675 auto prob22 = std::upper_bound(
676 pProbaShellMap[ionizationLevelIndex][(*k
677 pProbaShellMap[ionizationLevelIndex][(*k
678
679 auto prob21 = prob22 - 1;
680
681 valueK1 = *k1;
682 valueK2 = *k2;
683 valuePROB21 = *prob21;
684 valuePROB22 = *prob22;
685 valuePROB12 = *prob12;
686 valuePROB11 = *prob11;
687
688 nrjTransf11 = pNrjTransfData[ionizationLev
689 nrjTransf12 = pNrjTransfData[ionizationLev
690 nrjTransf21 = pNrjTransfData[ionizationLev
691 nrjTransf22 = pNrjTransfData[ionizationLev
692 }
693
694 // Avoids cases where cum xs is zero for k1
695 // k1<k2)
696
697 if(random > pProbaShellMap[ionizationLevelIn
698 {
699 auto prob22 = std::upper_bound(
700 pProbaShellMap[ionizationLevelIndex][(*k
701 pProbaShellMap[ionizationLevelIndex][(*k
702 auto prob21 = prob22 - 1;
703
704 valueK1 = *k1;
705 valueK2 = *k2;
706 valuePROB21 = *prob21;
707 valuePROB22 = *prob22;
708 nrjTransf21 = pNrjTransfData[ionizationLev
709 nrjTransf22 = pNrjTransfData[ionizationLev
710
711 G4double interpolatedvalue2 =
712 Interpolate(valuePROB21, valuePROB22, ra
713
714 G4double value = Interpolate(valueK1, valu
715 return value;
716 }
717 G4double nrjTransfProduct =
718 nrjTransf11 * nrjTransf12 * nrjTransf21 *
719
720 if(nrjTransfProduct != 0.)
721 {
722 nrj = QuadInterpolator(valuePROB11, valueP
723 nrjTransf11, nrjTra
724 valueK1, valueK2, k
725 }
726 return nrj;
727 }
728