| 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 "G4DNABornIonisationModel2.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 "G4DNABornAngle.hh"
36 #include "G4DeltaAngle.hh"
37 #include "G4Exp.hh"
38
39 //....oooOO0OOooo........oooOO0OOooo........oo
40
41 using namespace std;
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44
45 G4DNABornIonisationModel2::G4DNABornIonisation
46
47 G4VEmModel(nam)
48 {
49 verboseLevel = 0;
50 // Verbosity scale:
51 // 0 = nothing
52 // 1 = warning for energy non-conservation
53 // 2 = details of energy budget
54 // 3 = calculation of cross sections, file o
55 // 4 = entering in methods
56
57 if (verboseLevel > 0)
58 {
59 G4cout << "Born ionisation model is constr
60 }
61
62 // Mark this model as "applicable" for atomi
63
64 SetDeexcitationFlag(true);
65 fAtomDeexcitation = nullptr;
66 fParticleChangeForGamma = nullptr;
67 fpMolWaterDensity = nullptr;
68 fTableData = nullptr;
69 fLowEnergyLimit = 0;
70 fHighEnergyLimit = 0;
71 fParticleDef = nullptr;
72
73 // Define default angular generator
74
75 SetAngularDistribution(new G4DNABornAngle())
76
77 // Selection of computation method
78
79 fasterCode = false;
80
81 // Selection of stationary mode
82
83 statCode = false;
84
85 // Selection of SP scaling
86
87 spScaling = true;
88 }
89
90 //....oooOO0OOooo........oooOO0OOooo........oo
91
92 G4DNABornIonisationModel2::~G4DNABornIonisatio
93 {
94 // Cross section
95
96
97 delete fTableData;
98
99 // Final state
100
101 fVecm.clear();
102 }
103
104 //....oooOO0OOooo........oooOO0OOooo........oo
105
106 void G4DNABornIonisationModel2::Initialise(con
107 con
108 {
109
110 if (verboseLevel > 3)
111 {
112 G4cout << "Calling G4DNABornIonisationMode
113 }
114
115 if(fParticleDef != nullptr && particle != fP
116 {
117 G4ExceptionDescription description;
118 description << "You are trying to initiali
119 "for particle "
120 << particle->GetParticleName()
121 << G4endl;
122 description << "G4DNABornIonisationModel2
123 "for particle:" << fParticleDef->GetPartic
124 G4Exception("G4DNABornIonisationModel2::In
125 FatalException,description);
126 }
127
128 fParticleDef = particle;
129
130 // Energy limits
131 const char* path = G4FindDataDir("G4LEDATA")
132
133 // ***
134
135 G4String particleName = particle->GetParticl
136 std::ostringstream fullFileName;
137 fullFileName << path;
138
139 if(particleName == "e-")
140 {
141 fTableFile = "dna/sigma_ionisation_e_born"
142 fLowEnergyLimit = 11.*eV;
143 fHighEnergyLimit = 1.*MeV;
144
145 if (fasterCode)
146 {
147 fullFileName << "/dna/sigmadiff_cumulate
148 }
149 else
150 {
151 fullFileName << "/dna/sigmadiff_ionisati
152 }
153 }
154 else if(particleName == "proton")
155 {
156 fTableFile = "dna/sigma_ionisation_p_born"
157 fLowEnergyLimit = 500. * keV;
158 fHighEnergyLimit = 100. * MeV;
159
160 if (fasterCode)
161 {
162 fullFileName << "/dna/sigmadiff_cumulate
163 }
164 else
165 {
166 fullFileName << "/dna/sigmadiff_ionisati
167 }
168 }
169
170 // Cross section
171
172 G4double scaleFactor = (1.e-22 / 3.343) * m*
173 fTableData = new G4DNACrossSectionDataSet(ne
174 fTableData->LoadData(fTableFile);
175
176 // Final state
177
178 std::ifstream diffCrossSection(fullFileName.
179
180 if (!diffCrossSection)
181 {
182 G4ExceptionDescription description;
183 description << "Missing data file:" << G4e
184 G4Exception("G4DNABornIonisationModel2::In
185 FatalException,description);
186 }
187
188 // Clear the arrays for re-initialization ca
189 // March 25th, 2014 - Vaclav Stepan, Sebasti
190
191 fTdummyVec.clear();
192 fVecm.clear();
193
194 for (int j=0; j<5; j++)
195 {
196 fProbaShellMap[j].clear();
197 fDiffCrossSectionData[j].clear();
198 fNrjTransfData[j].clear();
199 }
200
201 //
202
203 fTdummyVec.push_back(0.);
204 while(!diffCrossSection.eof())
205 {
206 G4double tDummy;
207 G4double eDummy;
208 diffCrossSection>>tDummy>>eDummy;
209 if (tDummy != fTdummyVec.back()) fTdummyVe
210
211 G4double tmp;
212 for (int j=0; j<5; j++)
213 {
214 diffCrossSection>> tmp;
215
216 fDiffCrossSectionData[j][tDummy][eDummy]
217
218 if (fasterCode)
219 {
220 fNrjTransfData[j][tDummy][fDiffCrossSe
221 fProbaShellMap[j][tDummy].push_back(fD
222 }
223
224 // SI - only if eof is not reached
225 if (!diffCrossSection.eof() && !fasterCo
226
227 if (!fasterCode) fVecm[tDummy].push_back
228
229 }
230 }
231
232 //
233 SetLowEnergyLimit(fLowEnergyLimit);
234 SetHighEnergyLimit(fHighEnergyLimit);
235
236 if( verboseLevel>0 )
237 {
238 G4cout << "Born ionisation model is initia
239 << "Energy range: "
240 << LowEnergyLimit() / eV << " eV - "
241 << HighEnergyLimit() / keV << " keV for "
242 << particle->GetParticleName()
243 << G4endl;
244 }
245
246 // Initialize water density pointer
247
248 fpMolWaterDensity = G4DNAMolecularMaterial::
249 GetNumMolPerVolTableFor(G4Material::GetMater
250
251 // AD
252
253 fAtomDeexcitation = G4LossTableManager::Inst
254
255 if (isInitialised)
256 { return;}
257 fParticleChangeForGamma = GetParticleChangeF
258 isInitialised = true;
259 }
260
261 //....oooOO0OOooo........oooOO0OOooo........oo
262
263 G4double G4DNABornIonisationModel2::CrossSecti
264
265
266
267
268 {
269 if (verboseLevel > 3)
270 {
271 G4cout << "Calling CrossSectionPerVolume()
272 << G4endl;
273 }
274
275 if (particleDefinition != fParticleDef) retu
276
277 // Calculate total cross section for model
278
279 G4double sigma=0;
280
281 G4double waterDensity = (*fpMolWaterDensity)
282
283 if (ekin >= fLowEnergyLimit && ekin <= fHigh
284 {
285 sigma = fTableData->FindValue(ekin);
286
287 // ICRU49 electronic SP scaling - ZF, SI
288
289 if (particleDefinition == G4Proton::Proton
290 {
291 G4double A = 1.39241700556072800000E-009
292 G4double B = -8.52610412942622630000E-00
293 sigma = sigma * G4Exp(A*(ekin/eV)+B);
294 }
295 //
296 }
297
298 if (verboseLevel > 2)
299 {
300 G4cout << "_______________________________
301 G4cout << "G4DNABornIonisationModel2 - XS
302 G4cout << "Kinetic energy(eV)=" << ekin/eV
303 G4cout << "Cross section per water molecul
304 G4cout << "Cross section per water molecul
305 G4cout << "G4DNABornIonisationModel2 - XS
306 }
307
308 return sigma*waterDensity;
309 }
310
311 //....oooOO0OOooo........oooOO0OOooo........oo
312
313 void G4DNABornIonisationModel2::SampleSecondar
314
315
316
317
318 {
319
320 if (verboseLevel > 3)
321 {
322 G4cout << "Calling SampleSecondaries() of
323 << G4endl;
324 }
325
326 G4double k = particle->GetKineticEnergy();
327
328 if (k >= fLowEnergyLimit && k <= fHighEnergy
329 {
330 G4ParticleMomentum primaryDirection = part
331 G4double particleMass = particle->GetDefin
332 G4double totalEnergy = k + particleMass;
333 G4double pSquare = k * (totalEnergy + part
334 G4double totalMomentum = std::sqrt(pSquare
335
336 G4int ionizationShell = 0;
337
338 if (!fasterCode) ionizationShell = RandomS
339
340 // SI: The following protection is necessa
341 // sigmadiff_ionisation_e_born.dat has no
342 // sigmadiff_cumulated_ionisation_e_born.
343 // this is due to the fact that the max a
344 // strictly above this value have non zer
345
346 if (fasterCode)
347 do
348 {
349 ionizationShell = RandomSelect(k);
350 } while (k<19*eV && ionizationShell==2 &&
351
352 G4double secondaryKinetic=-1000*eV;
353
354 if (!fasterCode)
355 {
356 secondaryKinetic = RandomizeEjectedElect
357 }
358 else
359 {
360 secondaryKinetic = RandomizeEjectedElect
361 }
362
363 G4int Z = 8;
364
365 G4ThreeVector deltaDirection =
366 GetAngularDistribution()->SampleDirectionF
367 Z, ionizationShell,
368 couple->GetMaterial());
369
370 if (secondaryKinetic>0)
371 {
372 auto dp = new G4DynamicParticle (G4Elec
373 fvect->push_back(dp);
374 }
375
376 if (particle->GetDefinition() == G4Electro
377 {
378 G4double deltaTotalMomentum = std::sqrt(
379
380 G4double finalPx = totalMomentum*primary
381 G4double finalPy = totalMomentum*primary
382 G4double finalPz = totalMomentum*primary
383 G4double finalMomentum = std::sqrt(final
384 finalPx /= finalMomentum;
385 finalPy /= finalMomentum;
386 finalPz /= finalMomentum;
387
388 G4ThreeVector direction;
389 direction.set(finalPx,finalPy,finalPz);
390
391 fParticleChangeForGamma->ProposeMomentum
392 }
393
394 else fParticleChangeForGamma->ProposeMomen
395
396 // AM: sample deexcitation
397 // here we assume that H_{2}O electronic l
398 // this can be considered true with a roug
399
400 std::size_t secNumberInit = 0;
401 std::size_t secNumberFinal = 0;
402
403 G4double bindingEnergy = 0;
404 bindingEnergy = waterStructure.IonisationE
405
406 // SI: additional protection if tcs interp
407 if (k<bindingEnergy) return;
408 //
409
410 G4double scatteredEnergy = k-bindingEnergy
411
412 // SI: only atomic deexcitation from K she
413 if((fAtomDeexcitation != nullptr) && ioniz
414 {
415 const G4AtomicShell* shell =
416 fAtomDeexcitation->GetAtomicShell(Z, G
417 secNumberInit = fvect->size();
418 fAtomDeexcitation->GenerateParticles(fve
419 secNumberFinal = fvect->size();
420
421 if(secNumberFinal > secNumberInit)
422 {
423 for (std::size_t i=secNumberInit; i<secNumbe
424 {
425 //Check if there is enough residual
426 if (bindingEnergy >= ((*fvect)[i])->
427 {
428 //Ok, this is a valid secondary:
429 bindingEnergy -= ((*fvect)[i])->GetKine
430 }
431 else
432 {
433 //Invalid secondary: not enough energy
434 //Keep its energy in the local deposit
435 delete (*fvect)[i];
436 (*fvect)[i]=nullptr;
437 }
438 }
439 }
440
441 }
442
443 //This should never happen
444 if(bindingEnergy < 0.0)
445 G4Exception("G4DNAEmfietzoglouIonisatioMo
446 "em2050",FatalException,"Nega
447
448 //bindingEnergy has been decreased
449 //by the amount of energy taken away by de
450 if (!statCode)
451 {
452 fParticleChangeForGamma->SetProposedKine
453 fParticleChangeForGamma->ProposeLocalEne
454 }
455 else
456 {
457 fParticleChangeForGamma->SetProposedKine
458 fParticleChangeForGamma->ProposeLocalEne
459 }
460
461 // TEST //////////////////////////
462 // if (secondaryKinetic<0) abort();
463 // if (scatteredEnergy<0) abort();
464 // if (k-scatteredEnergy-secondaryKinetic-
465 // if (k-scatteredEnergy<0) abort();
466 /////////////////////////////////
467
468 const G4Track * theIncomingTrack = fPartic
469 G4DNAChemistryManager::Instance()->CreateW
470 ionizationShell,
471 theIncomingTrack);
472 }
473 }
474
475 //....oooOO0OOooo........oooOO0OOooo........oo
476
477 G4double G4DNABornIonisationModel2::RandomizeE
478
479
480 {
481 // G4cout << "*** SLOW computation for " <<
482
483 if (particleDefinition == G4Electron::Electr
484 {
485 G4double maximumEnergyTransfer = 0.;
486 if ((k + waterStructure.IonisationEnergy(s
487 maximumEnergyTransfer = k;
488 else
489 maximumEnergyTransfer = (k + waterStruct
490
491 // SI : original method
492 /*
493 G4double crossSectionMaximum = 0.;
494 for(G4double value=waterStructure.Ionisat
495 {
496 G4double differentialCrossSection = Diffe
497 if(differentialCrossSection >= crossSecti
498 }
499 */
500
501 // SI : alternative method
502 G4double crossSectionMaximum = 0.;
503
504 G4double minEnergy = waterStructure.Ionisa
505 G4double maxEnergy = maximumEnergyTransfer
506 G4int nEnergySteps = 50;
507
508 G4double value(minEnergy);
509 G4double stpEnergy(std::pow(maxEnergy / va
510 1. / static_ca
511 G4int step(nEnergySteps);
512 while (step > 0)
513 {
514 step--;
515 G4double differentialCrossSection =
516 DifferentialCrossSection(particleDef
517 k / eV,
518 value / eV,
519 shell);
520 if (differentialCrossSection >= crossSec
521 crossSectionMaximum = differentialCros
522 value *= stpEnergy;
523 }
524 //
525
526 G4double secondaryElectronKineticEnergy =
527 do
528 {
529 secondaryElectronKineticEnergy = G4Unifo
530 } while(G4UniformRand()*crossSectionMaximu
531 DifferentialCrossSection(particleDefin
532 (secondaryElectronKineticEnergy+wa
533
534 return secondaryElectronKineticEnergy;
535
536 }
537
538 if (particleDefinition == G4Proton::ProtonDe
539 {
540 G4double maximumKineticEnergyTransfer = 4.
541 * (electron_mass_c2 / proton_mass_c2)
542
543 G4double crossSectionMaximum = 0.;
544 for (G4double value = waterStructure.Ionis
545 value <= 4. * waterStructure.Ionisatio
546 {
547 G4double differentialCrossSection =
548 DifferentialCrossSection(particleDef
549 k / eV,
550 value / eV,
551 shell);
552 if (differentialCrossSection >= crossSec
553 crossSectionMaximum = differentialCros
554 }
555
556 G4double secondaryElectronKineticEnergy =
557 do
558 {
559 secondaryElectronKineticEnergy = G4Unifo
560 } while(G4UniformRand()*crossSectionMaximu
561 DifferentialCrossSection(particleDefin
562 (secondaryElectronKineticEnergy+wa
563
564 return secondaryElectronKineticEnergy;
565 }
566
567 return 0;
568 }
569
570 //....oooOO0OOooo........oooOO0OOooo........oo
571
572 // The following section is not used anymore b
573 // GetAngularDistribution()->SampleDirectionFo
574
575 /*
576 void G4DNABornIonisationModel2::RandomizeEjec
577 G4double k,
578 G4double secKinetic,
579 G4double & cosTheta,
580 G4double & phi )
581 {
582 if (particleDefinition == G4Electron::Electro
583 {
584 phi = twopi * G4UniformRand();
585 if (secKinetic < 50.*eV) cosTheta = (2.*G4Uni
586 else if (secKinetic <= 200.*eV)
587 {
588 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4
589 else cosTheta = G4UniformRand()*(std::sqrt(2.
590 }
591 else
592 {
593 G4double sin2O = (1.-secKinetic/k) / (1.+secK
594 cosTheta = std::sqrt(1.-sin2O);
595 }
596 }
597
598 else if (particleDefinition == G4Proton::Prot
599 {
600 G4double maxSecKinetic = 4.* (electron_mass_c
601 phi = twopi * G4UniformRand();
602
603 // cosTheta = std::sqrt(secKinetic / maxSecKi
604
605 // Restriction below 100 eV from Emfietzoglou
606
607 if (secKinetic>100*eV) cosTheta = std::sqrt(s
608 else cosTheta = (2.*G4UniformRand())-1.;
609
610 }
611 }
612 */
613
614 //....oooOO0OOooo........oooOO0OOooo........oo
615 G4double G4DNABornIonisationModel2::Differenti
616
617
618
619 {
620 G4double sigma = 0.;
621
622 if (energyTransfer >= waterStructure.Ionisat
623 {
624 G4double valueT1 = 0;
625 G4double valueT2 = 0;
626 G4double valueE21 = 0;
627 G4double valueE22 = 0;
628 G4double valueE12 = 0;
629 G4double valueE11 = 0;
630
631 G4double xs11 = 0;
632 G4double xs12 = 0;
633 G4double xs21 = 0;
634 G4double xs22 = 0;
635
636 // Protection against out of boundary acce
637 if (k==fTdummyVec.back()) k=k*(1.-1e-12);
638 //
639
640 // k should be in eV and energy transfer e
641
642 auto t2 = std::upper_bound(fTdummyVec.begi
643
644
645
646 auto t1 = t2 - 1;
647
648 // SI : the following condition avoids sit
649
650 if (energyTransfer <= fVecm[(*t1)].back()
651 && energyTransfer <= fVecm[(*t2)].back
652 {
653 auto e12 = std::upper_bound(fVecm[(*t1)]
654
655
656 auto e11 = e12 - 1;
657
658 auto e22 = std::upper_bound(fVecm[(*t2)]
659
660
661 auto e21 = e22 - 1;
662
663 valueT1 = *t1;
664 valueT2 = *t2;
665 valueE21 = *e21;
666 valueE22 = *e22;
667 valueE12 = *e12;
668 valueE11 = *e11;
669
670 xs11 = fDiffCrossSectionData[ionizationL
671 xs12 = fDiffCrossSectionData[ionizationL
672 xs21 = fDiffCrossSectionData[ionizationL
673 xs22 = fDiffCrossSectionData[ionizationL
674
675 }
676
677 G4double xsProduct = xs11 * xs12 * xs21 *
678 if (xsProduct != 0.)
679 {
680 sigma = QuadInterpolator(valueE11,
681 valueE12,
682 valueE21,
683 valueE22,
684 xs11,
685 xs12,
686 xs21,
687 xs22,
688 valueT1,
689 valueT2,
690 k,
691 energyTransfer)
692 }
693 }
694
695 return sigma;
696 }
697
698 //....oooOO0OOooo........oooOO0OOooo........oo
699
700 G4double G4DNABornIonisationModel2::Interpolat
701
702
703
704
705 {
706 G4double value = 0.;
707
708 // Log-log interpolation by default
709
710 if (e1 != 0 && e2 != 0 && (std::log10(e2) -
711 && !fasterCode)
712 {
713 G4double a = (std::log10(xs2) - std::log10
714 / (std::log10(e2) - std::log10(e1));
715 G4double b = std::log10(xs2) - a * std::lo
716 G4double sigma = a * std::log10(e) + b;
717 value = (std::pow(10., sigma));
718 }
719
720 // Switch to lin-lin interpolation
721 /*
722 if ((e2-e1)!=0)
723 {
724 G4double d1 = xs1;
725 G4double d2 = xs2;
726 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1)
727 }
728 */
729
730 // Switch to log-lin interpolation for faste
731 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 &
732 {
733 G4double d1 = std::log10(xs1);
734 G4double d2 = std::log10(xs2);
735 value = std::pow(10., (d1 + (d2 - d1) * (e
736 }
737
738 // Switch to lin-lin interpolation for faste
739 // in case one of xs1 or xs2 (=cum proba) va
740
741 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0)
742 {
743 G4double d1 = xs1;
744 G4double d2 = xs2;
745 value = (d1 + (d2 - d1) * (e - e1) / (e2 -
746 }
747
748 /*
749 G4cout
750 << e1 << " "
751 << e2 << " "
752 << e << " "
753 << xs1 << " "
754 << xs2 << " "
755 << value
756 << G4endl;
757 */
758
759 return value;
760 }
761
762 //....oooOO0OOooo........oooOO0OOooo........oo
763
764 G4double G4DNABornIonisationModel2::QuadInterp
765
766
767
768
769
770
771
772
773
774
775
776 {
777 G4double interpolatedvalue1 = Interpolate(e1
778 G4double interpolatedvalue2 = Interpolate(e2
779 G4double value = Interpolate(t1,
780 t2,
781 t,
782 interpolatedval
783 interpolatedval
784
785 return value;
786 }
787
788 //....oooOO0OOooo........oooOO0OOooo........oo
789
790 G4double G4DNABornIonisationModel2::GetPartial
791
792
793
794 {
795 return fTableData->GetComponent(level)->Find
796 }
797
798 //....oooOO0OOooo........oooOO0OOooo........oo
799
800 G4int G4DNABornIonisationModel2::RandomSelect(
801 {
802 G4int level = 0;
803
804 auto valuesBuffer = new G4double[fTableData
805 const auto n = (G4int)fTableData->NumberOfC
806 G4int i(n);
807 G4double value = 0.;
808
809 while (i > 0)
810 {
811 i--;
812 valuesBuffer[i] = fTableData->GetComponent
813 value += valuesBuffer[i];
814 }
815
816 value *= G4UniformRand();
817
818 i = n;
819
820 while (i > 0)
821 {
822 i--;
823
824 if (valuesBuffer[i] > value)
825 {
826 delete[] valuesBuffer;
827 return i;
828 }
829 value -= valuesBuffer[i];
830 }
831
832
833 delete[] valuesBuffer;
834
835 return level;
836 }
837
838 //....oooOO0OOooo........oooOO0OOooo........oo
839
840 G4double G4DNABornIonisationModel2::RandomizeE
841
842
843 {
844 // G4cout << "*** FAST computation for " <<
845
846 G4double secondaryElectronKineticEnergy = 0.
847
848 G4double random = G4UniformRand();
849
850 secondaryElectronKineticEnergy = TransferedE
851
852
853
854 - waterStructure.IonisationEnergy(shell)
855
856 // G4cout << TransferedEnergy(particleDefini
857 if (secondaryElectronKineticEnergy < 0.)
858 return 0.;
859 //
860
861 return secondaryElectronKineticEnergy;
862 }
863
864 //....oooOO0OOooo........oooOO0OOooo........oo
865
866 G4double G4DNABornIonisationModel2::Transfered
867
868
869
870 {
871
872 G4double nrj = 0.;
873
874 G4double valueK1 = 0;
875 G4double valueK2 = 0;
876 G4double valuePROB21 = 0;
877 G4double valuePROB22 = 0;
878 G4double valuePROB12 = 0;
879 G4double valuePROB11 = 0;
880
881 G4double nrjTransf11 = 0;
882 G4double nrjTransf12 = 0;
883 G4double nrjTransf21 = 0;
884 G4double nrjTransf22 = 0;
885
886 // Protection against out of boundary access
887 if (k==fTdummyVec.back()) k=k*(1.-1e-12);
888 //
889
890 // k should be in eV
891 auto k2 = std::upper_bound(fTdummyVec.begin(
892
893
894 auto k1 = k2 - 1;
895
896 /*
897 G4cout << "----> k=" << k
898 << " " << *k1
899 << " " << *k2
900 << " " << random
901 << " " << ionizationLevelIndex
902 << " " << eProbaShellMap[ionizationLevelInd
903 << " " << eProbaShellMap[ionizationLevelInd
904 << G4endl;
905 */
906
907 // SI : the following condition avoids situa
908 if (random <= fProbaShellMap[ionizationLevel
909 && random <= fProbaShellMap[ionizationLe
910 {
911 auto prob12 =
912 std::upper_bound(fProbaShellMap[ioniza
913 fProbaShellMap[ioniza
914 random);
915
916 auto prob11 = prob12 - 1;
917
918 auto prob22 =
919 std::upper_bound(fProbaShellMap[ioniza
920 fProbaShellMap[ioniza
921 random);
922
923 auto prob21 = prob22 - 1;
924
925 valueK1 = *k1;
926 valueK2 = *k2;
927 valuePROB21 = *prob21;
928 valuePROB22 = *prob22;
929 valuePROB12 = *prob12;
930 valuePROB11 = *prob11;
931
932 /*
933 G4cout << " " << random << " " <<
934 << valuePROB12 << " " << valuePROB21 << "
935 */
936
937 nrjTransf11 = fNrjTransfData[ionizationLev
938 nrjTransf12 = fNrjTransfData[ionizationLev
939 nrjTransf21 = fNrjTransfData[ionizationLev
940 nrjTransf22 = fNrjTransfData[ionizationLev
941
942 /*
943 G4cout << " " << ionizationLevelIn
944 << random << " " <<valueK1 << " " << valu
945
946 G4cout << " " << random << " " <<
947 << nrjTransf12 << " " << nrjTransf21 << "
948 */
949
950 }
951 // Avoids cases where cum xs is zero for k1
952 if (random > fProbaShellMap[ionizationLevelI
953 {
954 auto prob22 =
955 std::upper_bound(fProbaShellMap[ioniza
956 fProbaShellMap[ioniza
957 random);
958
959 auto prob21 = prob22 - 1;
960
961 valueK1 = *k1;
962 valueK2 = *k2;
963 valuePROB21 = *prob21;
964 valuePROB22 = *prob22;
965
966 // G4cout << " " << random << " " <
967
968 nrjTransf21 = fNrjTransfData[ionizationLev
969 nrjTransf22 = fNrjTransfData[ionizationLev
970
971 G4double interpolatedvalue2 = Interpolate(
972
973
974
975
976
977 // zeros are explicitly set
978
979 G4double value = Interpolate(valueK1, valu
980
981 /*
982 G4cout << " " << ionizationLevelIn
983 << random << " " <<valueK1 << " " << valu
984
985 G4cout << " " << random << " " <<
986 << nrjTransf12 << " " << nrjTransf21 << "
987
988 G4cout << "ici" << " " << value << G4endl
989 */
990
991 return value;
992 }
993
994 // End electron and proton cases
995
996 G4double nrjTransfProduct = nrjTransf11 * nr
997 * nrjTransf22;
998 //G4cout << "nrjTransfProduct=" << nrjTransf
999
1000 if (nrjTransfProduct != 0.)
1001 {
1002 nrj = QuadInterpolator(valuePROB11,
1003 valuePROB12,
1004 valuePROB21,
1005 valuePROB22,
1006 nrjTransf11,
1007 nrjTransf12,
1008 nrjTransf21,
1009 nrjTransf22,
1010 valueK1,
1011 valueK2,
1012 k,
1013 random);
1014 }
1015 // G4cout << nrj << endl;
1016
1017 return nrj;
1018 }
1019