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