| 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 // Reference:
27 // A.D. Dominguez-Munoz, M.I. Gallardo, M.C. Bordage,
28 // Z. Francis, S. Incerti, M.A. Cortes-Giraldo,
29 // Radiat. Phys. Chem. 199 (2022) 110363.
30 //
31 // Class authors:
32 // A.D. Dominguez-Munoz
33 // M.A. Cortes-Giraldo (miancortes -at- us.es)
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........oooOO0OOooo........oooOO0OOooo....
50 G4DNARPWBAIonisationModel::G4DNARPWBAIonisationModel(
51 const G4ParticleDefinition*, const G4String& nam)
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 openings, sampling of atoms
59 // 4 = entering in methods
60
61 if(verboseLevel > 0)
62 {
63 G4cout << "RPWBA ionisation model is constructed " << G4endl;
64 }
65 SetDeexcitationFlag(true);
66 SetAngularDistribution(new G4DNABornAngle());
67 }
68
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70
71 G4DNARPWBAIonisationModel::~G4DNARPWBAIonisationModel()
72 {
73 eVecm.clear();
74 pVecm.clear();
75 }
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78 G4bool G4DNARPWBAIonisationModel::InEnergyLimit(const G4double& k)
79 {
80 if(lowEnergyLimit == highEnergyLimit)
81 {
82 G4Exception("G4DNARPWBAIonisationModel::InEnergyLimit", "em0102",
83 FatalException, "lowEnergyLimit == highEnergyLimit");
84 }
85 return k >= lowEnergyLimit && k <= highEnergyLimit;
86 }
87
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 void G4DNARPWBAIonisationModel::InitialiseForProton(
90 const G4ParticleDefinition* part)
91 {
92 if(part != fProtonDef)
93 {
94 G4Exception("G4DNARPWBAIonisationModel::InitialiseForProton", "em0002",
95 FatalException, "Model not applicable to particle type.");
96 }
97 // Energy limits
98 G4String fileProton("dna/sigma_ionisation_p_RPWBA");
99 G4double scaleFactor = 1 * cm * cm;
100 const char *path = G4FindDataDir("G4LEDATA");
101 lowEnergyLimit = 100. * MeV;
102 highEnergyLimit = 300. * MeV;
103
104 if(LowEnergyLimit() < lowEnergyLimit || HighEnergyLimit() > highEnergyLimit)
105 {
106 G4ExceptionDescription ed;
107 ed << "Model is applicable from "<<lowEnergyLimit<<" to "<<highEnergyLimit;
108 G4Exception("G4DNARPWBAIonisationModel::InitialiseForProton", "em0004",
109 FatalException, ed);
110 }
111
112 fpTotalCrossSection = make_unique<G4DNACrossSectionDataSet>(
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_cumulated_ionisation_p_RPWBA.dat"
121 : pFullFileName << path << "/dna/sigmadiff_ionisation_p_RPWBA.dat";
122 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
123 if(!pDiffCrossSection)
124 {
125 G4ExceptionDescription exceptionDescription;
126 exceptionDescription << "Missing data file: " + pFullFileName.str();
127 G4Exception("G4DNARPWBAIonisationModel::InitialiseForProton", "em0003",
128 FatalException, exceptionDescription);
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 >> pDiffCrossSectionData[j][tDummy][eDummy];
145
146 if(fasterCode)
147 {
148 pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]] =
149 eDummy;
150 pProbaShellMap[j][tDummy].push_back(
151 pDiffCrossSectionData[j][tDummy][eDummy]);
152 }
153
154 // SI - only if eof is not reached !
155 if(!pDiffCrossSection.eof() && !fasterCode)
156 {
157 pDiffCrossSectionData[j][tDummy][eDummy] *= scaleFactor;
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(const G4ParticleDefinition* particle,
173 const G4DataVector& /*cuts*/)
174 {
175 if(isInitialised)
176 {
177 return;
178 }
179 if(verboseLevel > 3)
180 {
181 G4cout << "Calling G4DNARPWBAIonisationModel::Initialise()"
182 << particle->GetParticleName() << G4endl;
183 }
184
185 InitialiseForProton(particle);
186
187 if(verboseLevel > 0)
188 {
189 G4cout << "RPWBA ionisation model is initialized " << G4endl
190 << "Energy range: " << LowEnergyLimit() / MeV << " MeV - "
191 << HighEnergyLimit() / MeV << " MeV for "
192 << particle->GetParticleName() << G4endl;
193 }
194
195 // Initialize water density pointer
196 if(G4Material::GetMaterial("G4_WATER") != nullptr)
197 {
198 fpMolWaterDensity =
199 G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(
200 G4Material::GetMaterial("G4_WATER"));
201 }
202 else
203 {
204 G4ExceptionDescription exceptionDescription;
205 exceptionDescription << "G4_WATER does not exist :";
206 G4Exception("G4DNARPWBAIonisationModel::Initialise", "em00020",
207 FatalException, exceptionDescription);
208 }
209 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
210 fParticleChangeForGamma = GetParticleChangeForGamma();
211 isInitialised = true;
212 }
213
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
215
216 G4double G4DNARPWBAIonisationModel::CrossSectionPerVolume(
217 const G4Material* material, const G4ParticleDefinition* particleDefinition,
218 G4double ekin, G4double, G4double)
219 {
220 if(particleDefinition != fProtonDef)
221 {
222 G4Exception("G4DNARPWBAIonisationModel::CrossSectionPerVolume", "em0402",
223 FatalException, "Model not applicable to particle type.");
224 }
225 if(verboseLevel > 3)
226 {
227 G4cout << "Calling CrossSectionPerVolume() of G4DNARPWBAIonisationModel"
228 << G4endl;
229 }
230 G4double sigma;
231 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
232
233 if(InEnergyLimit(ekin))
234 {
235 sigma = fpTotalCrossSection->FindValue(ekin);
236 }
237 else
238 {
239 // nput energy is outside this interval the cross section is set to zero
240 // should add a warning or exception ?
241 return 0;
242 }
243
244 if(verboseLevel > 2)
245 {
246 G4cout << "__________________________________" << G4endl;
247 G4cout << "G4DNARPWBAIonisationModel - XS INFO START" << G4endl;
248 G4cout << "Kinetic energy(eV)=" << ekin / eV
249 << " particle : " << fProtonDef->GetParticleName() << G4endl;
250 G4cout << "Cross section per water molecule (cm^2)=" << sigma / cm / cm
251 << G4endl;
252 G4cout << "Cross section per water molecule (cm^-1)="
253 << sigma * waterDensity / (1. / cm) << G4endl;
254 G4cout << "G4DNARPWBAIonisationModel - XS INFO END" << G4endl;
255 }
256
257 return sigma * waterDensity;
258 }
259
260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
261
262 void G4DNARPWBAIonisationModel::SampleSecondaries(
263 std::vector<G4DynamicParticle*>* fvect, const G4MaterialCutsCouple* couple,
264 const G4DynamicParticle* particle, G4double, G4double)
265 {
266 if(verboseLevel > 3)
267 {
268 G4cout << "Calling SampleSecondaries() of G4DNARPWBAIonisationModel"
269 << G4endl;
270 }
271 G4double k = particle->GetKineticEnergy();
272 if(InEnergyLimit(k))
273 {
274 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
275 G4double particleMass = particle->GetDefinition()->GetPDGMass();
276 G4double totalEnergy = k + particleMass;
277 G4double pSquare = k * (totalEnergy + particleMass);
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 == 2 &&
292 particle->GetDefinition() == G4Electron::ElectronDefinition());
293 }
294
295 G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
296
297 // SI: additional protection if tcs interpolation method is modified
298 if(k < bindingEnergy)
299 {
300 return;
301 }
302 //
303 G4double secondaryKinetic;
304 if(!fasterCode)
305 {
306 secondaryKinetic = RandomizeEjectedElectronEnergy(k, ionizationShell);
307 }
308 else
309 {
310 secondaryKinetic =
311 RandomizeEjectedElectronEnergyFromCumulatedDcs(k, ionizationShell);
312 }
313
314 G4int Z = 8; // water Z (6 Oxygen + 2 hydrogen)
315 G4ThreeVector deltaDirection =
316 GetAngularDistribution()->SampleDirectionForShell(
317 particle, secondaryKinetic, Z, ionizationShell, couple->GetMaterial());
318
319 if(secondaryKinetic > 0){
320 auto dp = new G4DynamicParticle(G4Electron::Electron(), deltaDirection,
321 secondaryKinetic);
322 fvect->push_back(dp);
323 }
324
325 if(particle->GetDefinition() == G4Electron::ElectronDefinition()){
326 G4double deltaTotalMomentum = std::sqrt(
327 secondaryKinetic * (secondaryKinetic + 2. * electron_mass_c2));
328
329 G4double finalPx = totalMomentum * primaryDirection.x() -
330 deltaTotalMomentum * deltaDirection.x();
331 G4double finalPy = totalMomentum * primaryDirection.y() -
332 deltaTotalMomentum * deltaDirection.y();
333 G4double finalPz = totalMomentum * primaryDirection.z() -
334 deltaTotalMomentum * deltaDirection.z();
335 G4double finalMomentum =
336 std::sqrt(finalPx * finalPx + finalPy * finalPy + finalPz * finalPz);
337 finalPx /= finalMomentum;
338 finalPy /= finalMomentum;
339 finalPz /= finalMomentum;
340 G4ThreeVector direction;
341 direction.set(finalPx, finalPy, finalPz);
342 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
343 }
344 else
345 {
346 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
347 }
348
349 // AM: sample deexcitation
350 // here we assume that H_{2}O electronic levels are the same as Oxygen.
351 // this can be considered true with a rough 10% error in energy on K-shell,
352
353 size_t secNumberInit; // need to know at a certain point the energy of
354 // secondaries
355 size_t
356 secNumberFinal; // So I'll make the diference and then sum the energies
357
358 G4double scatteredEnergy = k - bindingEnergy - secondaryKinetic;
359
360 // SI: only atomic deexcitation from K shell is considered
361 if((fAtomDeexcitation != nullptr) && ionizationShell == 4)
362 {
363 const G4AtomicShell* shell =
364 fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
365 secNumberInit = fvect->size();
366 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
367 secNumberFinal = fvect->size();
368
369 if(secNumberFinal > secNumberInit){
370 for(size_t i = secNumberInit; i < secNumberFinal; ++i){
371 if(bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
372 {
373 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
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::SampleSecondaries()", "em2050",
386 FatalException, "Negative local energy deposit");
387 }
388
389 // bindingEnergy has been decreased
390 // by the amount of energy taken away by deexc. products
391 if(!statCode){
392 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
393 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
394 }else{
395 fParticleChangeForGamma->SetProposedKineticEnergy(k);
396 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k - scatteredEnergy);
397 }
398 const G4Track* theIncomingTrack =
399 fParticleChangeForGamma->GetCurrentTrack();
400 G4DNAChemistryManager::Instance()->CreateWaterMolecule(
401 eIonizedMolecule, ionizationShell, theIncomingTrack);
402 }
403 }
404
405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
406
407 G4double G4DNARPWBAIonisationModel::RandomizeEjectedElectronEnergy(
408 const G4double& k, const G4int& shell)
409 {
410 G4double maximumKineticEnergyTransfer =
411 4. * (electron_mass_c2 / proton_mass_c2) * k;
412 G4double IonisationEnergyInShell = waterStructure.IonisationEnergy(shell);
413 G4double kIneV = k / eV;
414
415 G4double crossSectionMaximum = 0.;
416 for(G4double value = IonisationEnergyInShell;
417 value <= 4. * IonisationEnergyInShell; value += 0.1 * eV)
418 {
419 G4double differentialCrossSection =
420 DifferentialCrossSection(kIneV, value / eV, shell);
421 if(differentialCrossSection >= crossSectionMaximum)
422 {
423 crossSectionMaximum = differentialCrossSection;
424 }
425 }
426
427 G4double secondaryElectronKineticEnergy;
428 do
429 {
430 secondaryElectronKineticEnergy =
431 G4UniformRand() * maximumKineticEnergyTransfer;
432 } while(G4UniformRand() * crossSectionMaximum >=
433 DifferentialCrossSection(kIneV,
434 (secondaryElectronKineticEnergy +
435 IonisationEnergyInShell) / eV, shell));
436
437 return secondaryElectronKineticEnergy;
438 }
439
440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
441 G4double G4DNARPWBAIonisationModel::DifferentialCrossSection(
442 const G4double& kine, const G4double& energyTransfer,
443 const G4int& ionizationLevelIndex)
444 {
445 G4double k = kine;
446 G4double sigma = 0.;
447 if(energyTransfer >=
448 waterStructure.IonisationEnergy(ionizationLevelIndex) / eV)
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 access - proton case : 100 MeV
463
464 if(k == pTdummyVec.back())
465 {
466 k = k * (1. - 1e-12);
467 }
468 // k should be in eV and energy transfer eV also
469 auto t2 = std::upper_bound(pTdummyVec.begin(), pTdummyVec.end(), k);
470 auto t1 = t2 - 1;
471
472 auto e12 = std::upper_bound(pVecm[(*t1)].begin(), pVecm[(*t1)].end(),
473 energyTransfer);
474 auto e11 = e12 - 1;
475
476 auto e22 = std::upper_bound(pVecm[(*t2)].begin(), pVecm[(*t2)].end(),
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[ionizationLevelIndex][valueT1][valueE11];
488 xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
489 xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
490 xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
491
492 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
493 if(xsProduct != 0.)
494 {
495 sigma =
496 QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12,
497 xs21, xs22, valueT1, valueT2, k, energyTransfer);
498 }
499 }
500 return sigma;
501 }
502
503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
504
505 G4double G4DNARPWBAIonisationModel::Interpolate(const G4double& e1,
506 const G4double& e2,
507 const G4double& e,
508 const G4double& xs1,
509 const G4double& xs2)
510 {
511 G4double value = 0.;
512
513 // Log-log interpolation by default
514
515 if(e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0 &&
516 !fasterCode)
517 {
518 G4double a =
519 (std::log10(xs2) - std::log10(xs1)) / (std::log10(e2) - std::log10(e1));
520 G4double b = std::log10(xs2) - a * std::log10(e2);
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 faster code
535 if((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
536 {
537 G4double d1 = std::log10(xs1);
538 G4double d2 = std::log10(xs2);
539 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
540 }
541
542 // Switch to lin-lin interpolation for faster code
543 // in case one of xs1 or xs2 (=cum proba) value is zero
544
545 if((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
546 {
547 G4double d1 = xs1;
548 G4double d2 = xs2;
549 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
550 }
551 return value;
552 }
553
554 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
555
556 G4double G4DNARPWBAIonisationModel::QuadInterpolator(
557 const G4double& e11, const G4double& e12, const G4double& e21,
558 const G4double& e22, const G4double& xs11, const G4double& xs12,
559 const G4double& xs21, const G4double& xs22, const G4double& t1,
560 const G4double& t2, const G4double& t, const G4double& e)
561 {
562 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
563 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
564 G4double value =
565 Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
566
567 return value;
568 }
569
570 G4double G4DNARPWBAIonisationModel::GetPartialCrossSection(
571 const G4Material* /*material*/, G4int level,
572 const G4ParticleDefinition* particle, G4double kineticEnergy)
573 {
574 if(fpTotalCrossSection != nullptr && particle != fProtonDef)
575 {
576 G4Exception("G4DNARPWBAIonisationModel::GetPartialCrossSection", "em0010",
577 FatalException, "Model not applicable to particle type.");
578 }
579 return fpTotalCrossSection->GetComponent(level)->FindValue(kineticEnergy);
580 }
581
582 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
583
584 G4int G4DNARPWBAIonisationModel::RandomSelect(G4double k)
585 {
586 if(fpTotalCrossSection == nullptr)
587 {
588 G4Exception("G4DNARPWBAIonisationModel::RandomSelect", "em0010",
589 FatalException, "Model not applicable to particle type.");
590 }
591 else
592 {
593 auto valuesBuffer = new G4double[fpTotalCrossSection->NumberOfComponents()];
594 const auto n = (G4int)fpTotalCrossSection->NumberOfComponents();
595 G4int i(n);
596 G4double value = 0.;
597 while(i > 0)
598 {
599 --i;
600 valuesBuffer[i] = fpTotalCrossSection->GetComponent(i)->FindValue(k);
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........oooOO0OOooo........oooOO0OOooo......
622
623 G4double
624 G4DNARPWBAIonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(
625 const G4double& k, const G4int& shell)
626 {
627 G4double random = G4UniformRand();
628 G4double secondaryKineticEnergy =
629 TransferedEnergy(k / eV, shell, random) * eV -
630 waterStructure.IonisationEnergy(shell);
631 if(secondaryKineticEnergy < 0.)
632 {
633 return 0.;
634 }
635 return secondaryKineticEnergy;
636 }
637
638 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
639
640 G4double G4DNARPWBAIonisationModel::TransferedEnergy(G4double k,
641 G4int ionizationLevelIndex,
642 const G4double& random)
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 - proton case : 100 MeV
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(), pTdummyVec.end(), k);
663 auto k1 = k2 - 1;
664
665 // SI : the following condition avoids situations where random > last vector
666 // element,
667 // for eg. when the last element is zero
668 if(random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back() &&
669 random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
670 {
671 auto prob12 = std::upper_bound(
672 pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
673 pProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
674 auto prob11 = prob12 - 1;
675 auto prob22 = std::upper_bound(
676 pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
677 pProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
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[ionizationLevelIndex][valueK1][valuePROB11];
689 nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
690 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
691 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
692 }
693
694 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always
695 // k1<k2)
696
697 if(random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
698 {
699 auto prob22 = std::upper_bound(
700 pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
701 pProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
702 auto prob21 = prob22 - 1;
703
704 valueK1 = *k1;
705 valueK2 = *k2;
706 valuePROB21 = *prob21;
707 valuePROB22 = *prob22;
708 nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
709 nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
710
711 G4double interpolatedvalue2 =
712 Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
713
714 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
715 return value;
716 }
717 G4double nrjTransfProduct =
718 nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
719
720 if(nrjTransfProduct != 0.)
721 {
722 nrj = QuadInterpolator(valuePROB11, valuePROB12, valuePROB21, valuePROB22,
723 nrjTransf11, nrjTransf12, nrjTransf21, nrjTransf22,
724 valueK1, valueK2, k, random);
725 }
726 return nrj;
727 }
728