| 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 // Authors: S. Meylan and C. Villagrasa (IRSN,
27 // Models come from
28 // M. Bug et al, Rad. Phys and Chem. 130, 459-
29 //
30
31 #include "G4DNAPTBIonisationModel.hh"
32
33 #include "G4DNAChemistryManager.hh"
34 #include "G4DNAMaterialManager.hh"
35 #include "G4LossTableManager.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4UAtomicDeexcitation.hh"
39 G4DNAPTBIonisationModel::G4DNAPTBIonisationMod
40
41
42 : G4VDNAModel(nam, applyToMaterial)
43 {
44 if (isAuger) {
45 // create the PTB Auger model
46 fpDNAPTBAugerModel = std::make_unique<G4DN
47 }
48 fpTHF = G4Material::GetMaterial("THF", false
49 fpPY = G4Material::GetMaterial("PY", false);
50 fpPU = G4Material::GetMaterial("PU", false);
51 fpTMP = G4Material::GetMaterial("TMP", false
52 fpG4_WATER = G4Material::GetMaterial("G4_WAT
53 fpBackbone_THF = G4Material::GetMaterial("ba
54 fpCytosine_PY = G4Material::GetMaterial("cyt
55 fpThymine_PY = G4Material::GetMaterial("thym
56 fpAdenine_PU = G4Material::GetMaterial("aden
57 fpBackbone_TMP = G4Material::GetMaterial("ba
58 fpGuanine_PU = G4Material::GetMaterial("guan
59 fpN2 = G4Material::GetMaterial("N2", false);
60 }
61
62 //....oooOO0OOooo........oooOO0OOooo........oo
63
64 void G4DNAPTBIonisationModel::Initialise(const
65 const
66 {
67 if (isInitialised) {
68 return;
69 }
70 if (verboseLevel > 3) {
71 G4cout << "Calling G4DNAPTBIonisationModel
72 }
73
74 G4double scaleFactor = 1e-16 * cm * cm;
75 G4double scaleFactorBorn = (1.e-22 / 3.343)
76
77 G4ParticleDefinition* electronDef = G4Electr
78 G4ParticleDefinition* protonDef = G4Proton::
79
80 //******************************************
81 // Cross section data
82 //******************************************
83 std::size_t index;
84 if (particle == electronDef) {
85 // Raw materials
86 //
87 // MPietrzak
88 if (fpN2 != nullptr) {
89 index = fpN2->GetIndex();
90 AddCrossSectionData(index, particle, "dn
91 "dna/sigmadiff_cumul
92 SetLowELimit(index, particle, 15.5 * eV)
93 SetHighELimit(index, particle, 1.02 * Me
94 }
95
96 // MPietrzak
97
98 if (fpTHF != nullptr) {
99 index = fpTHF->GetIndex();
100 AddCrossSectionData(index, particle, "dn
101 "dna/sigmadiff_cumul
102 SetLowELimit(index, particle, 12. * eV);
103 SetHighELimit(index, particle, 1. * keV)
104 }
105
106 if (fpPY != nullptr) {
107 index = fpPY->GetIndex();
108 AddCrossSectionData(index, particle, "dn
109 "dna/sigmadiff_cumul
110 SetLowELimit(index, particle, 12. * eV);
111 SetHighELimit(index, particle, 1. * keV)
112 }
113
114 if (fpPU != nullptr) {
115 index = fpPU->GetIndex();
116 AddCrossSectionData(index, particle, "dn
117 "dna/sigmadiff_cumul
118 SetLowELimit(index, particle, 12. * eV);
119 SetHighELimit(index, particle, 1. * keV)
120 }
121 if (fpTMP != nullptr) {
122 index = fpTMP->GetIndex();
123 AddCrossSectionData(index, particle, "dn
124 "dna/sigmadiff_cumul
125 SetLowELimit(index, particle, 12. * eV);
126 SetHighELimit(index, particle, 1. * keV)
127 }
128
129 if (fpG4_WATER != nullptr) {
130 index = fpG4_WATER->GetIndex();
131 AddCrossSectionData(index, particle, "dn
132 "dna/sigmadiff_ionis
133 SetLowELimit(index, particle, 12. * eV);
134 SetHighELimit(index, particle, 1. * keV)
135 }
136 // DNA materials
137 //
138 if (fpBackbone_THF != nullptr) {
139 index = fpBackbone_THF->GetIndex();
140 AddCrossSectionData(index, particle, "dn
141 "dna/sigmadiff_cumul
142 SetLowELimit(index, particle, 12. * eV);
143 SetHighELimit(index, particle, 1. * keV)
144 }
145
146 if (fpCytosine_PY != nullptr) {
147 index = fpCytosine_PY->GetIndex();
148 AddCrossSectionData(index, particle, "dn
149 "dna/sigmadiff_cumul
150 SetLowELimit(index, particle, 12. * eV);
151 SetHighELimit(index, particle, 1. * keV)
152 }
153
154 if (fpThymine_PY != nullptr) {
155 index = fpThymine_PY->GetIndex();
156 AddCrossSectionData(index, particle, "dn
157 "dna/sigmadiff_cumul
158 SetLowELimit(index, particle, 12. * eV);
159 SetHighELimit(index, particle, 1. * keV)
160 }
161
162 if (fpAdenine_PU != nullptr) {
163 index = fpAdenine_PU->GetIndex();
164 AddCrossSectionData(index, particle, "dn
165 "dna/sigmadiff_cumul
166 SetLowELimit(index, particle, 12. * eV);
167 SetHighELimit(index, particle, 1. * keV)
168 }
169
170 if (fpGuanine_PU != nullptr) {
171 index = fpGuanine_PU->GetIndex();
172 AddCrossSectionData(index, particle, "dn
173 "dna/sigmadiff_cumul
174 SetLowELimit(index, particle, 12. * eV);
175 SetHighELimit(index, particle, 1. * keV)
176 }
177
178 if (fpBackbone_TMP != nullptr) {
179 index = fpBackbone_TMP->GetIndex();
180 AddCrossSectionData(index, particle, "dn
181 "dna/sigmadiff_cumul
182 SetLowELimit(index, particle, 12. * eV);
183 SetHighELimit(index, particle, 1. * keV)
184 }
185 }
186
187 else if (particle == protonDef) {
188 G4String particleName = particle->GetParti
189
190 // Raw materials
191 //
192 if (fpTHF != nullptr) {
193 index = fpTHF->GetIndex();
194 AddCrossSectionData(index, particle, "dn
195 "dna/sigmadiff_cumul
196 SetLowELimit(index, particle, 70. * keV)
197 SetHighELimit(index, particle, 10. * MeV
198 }
199
200 if (fpPY != nullptr) {
201 index = fpPY->GetIndex();
202 AddCrossSectionData(index, particle, "dn
203 "dna/sigmadiff_cumul
204 SetLowELimit(index, particle, 70. * keV)
205 SetHighELimit(index, particle, 10. * MeV
206 }
207 /*
208 AddCrossSectionData("PU",
209 particleName,
210 "dna/sigma_ion
211 "dna/sigmadiff
212 scaleFactor);
213 SetLowELimit("PU", particleName2,
214 SetHighELimit("PU", particleName2,
215 */
216
217 if (fpTMP != nullptr) {
218 index = fpTMP->GetIndex();
219 AddCrossSectionData(index, particle, "dn
220 "dna/sigmadiff_cumul
221 SetLowELimit(index, particle, 70. * keV)
222 SetHighELimit(index, particle, 10. * MeV
223 }
224 }
225 // *****************************************
226 // deal with composite materials
227 // *****************************************
228 if (!G4DNAMaterialManager::Instance()->IsLoc
229 LoadCrossSectionData(particle);
230 G4DNAMaterialManager::Instance()->SetMaste
231 fpModelData = this;
232 }
233 else {
234 auto dataModel = dynamic_cast<G4DNAPTBIoni
235 G4DNAMaterialManager::Instance()->GetMod
236 if (dataModel == nullptr) {
237 G4cout << "G4DNAPTBIonisationModel::Init
238 G4Exception("G4DNAPTBIonisationModel::In
239 "not good modelData");
240 }
241 else {
242 fpModelData = dataModel;
243 }
244 }
245 // initialise DNAPTBAugerModel
246 if (fpDNAPTBAugerModel) {
247 fpDNAPTBAugerModel->Initialise();
248 }
249 fParticleChangeForGamma = GetParticleChangeF
250 isInitialised = true;
251 }
252
253 //....oooOO0OOooo........oooOO0OOooo........oo
254
255 G4double G4DNAPTBIonisationModel::CrossSection
256
257
258
259 {
260 // initialise the cross section value (outpu
261 G4double sigma(0);
262
263 // Get the current particle name
264 const G4String& particleName = p->GetParticl
265 const std::size_t& matID = material->GetInde
266
267 // Set the low and high energy limits
268 G4double lowLim = fpModelData->GetLowELimit(
269 G4double highLim = fpModelData->GetHighELimi
270
271 // Check that we are in the correct energy r
272 if (ekin >= lowLim && ekin < highLim) {
273 // Get the map with all the model data tab
274 auto tableData = fpModelData->GetData();
275 if ((*tableData)[matID][p] == nullptr) {
276 G4Exception("G4DNAPTBIonisationModel::Cr
277 "No model is registered");
278 }
279 // Retrieve the cross section value for th
280 sigma = (*tableData)[matID][p]->FindValue(
281
282 if (verboseLevel > 2) {
283 G4cout << "_____________________________
284 G4cout << "°°° G4DNAPTBIonisationMode
285 G4cout << "°°° Kinetic energy(eV)=" <
286 G4cout << "°°° Cross section per " <<
287 << G4endl;
288 G4cout << "°°° G4DNAPTBIonisationMode
289 }
290 }
291
292 // Return the cross section value
293 auto MolDensity =
294 (*G4DNAMolecularMaterial::Instance()->GetN
295 return sigma * MolDensity;
296 }
297
298 //....oooOO0OOooo........oooOO0OOooo........oo
299
300 void G4DNAPTBIonisationModel::SampleSecondarie
301
302
303
304 {
305 // Get the current particle energy
306 G4double k = aDynamicParticle->GetKineticEne
307 const std::size_t& materialID = pCouple->Get
308
309 // Get the current particle name
310 const auto& p = aDynamicParticle->GetDefinit
311 const auto& materialName = pCouple->GetMater
312 // Get the energy limits
313 G4double lowLim = fpModelData->GetLowELimit(
314 G4double highLim = fpModelData->GetHighELimi
315
316 // Check if we are in the correct energy ran
317 if (k >= lowLim && k < highLim) {
318 G4ParticleMomentum primaryDirection = aDyn
319 G4double particleMass = aDynamicParticle->
320 G4double totalEnergy = k + particleMass;
321 G4double pSquare = k * (totalEnergy + part
322 G4double totalMomentum = std::sqrt(pSquare
323
324 // Get the ionisation shell from a random
325 G4int ionizationShell = fpModelData->Rando
326
327 // Get the binding energy from the ptbStru
328 G4double bindingEnergy = ptbStructure.Ioni
329
330 // Initialize the secondary kinetic energy
331 G4double secondaryKinetic(-1000 * eV);
332
333 if (fpG4_WATER == nullptr || materialID !=
334 // Get the energy of the secondary parti
335 secondaryKinetic = fpModelData->Randomiz
336 aDynamicParticle->GetDefinition(), k /
337 }
338 else {
339 secondaryKinetic = fpModelData->Randomiz
340 aDynamicParticle->GetDefinition(), k,
341 }
342
343 if (secondaryKinetic <= 0) {
344 G4cout << "Fatal error *****************
345 << G4endl;
346 G4cout << "secondaryKinetic: " << second
347 G4cout << "k: " << k / eV << G4endl;
348 G4cout << "shell: " << ionizationShell <
349 G4cout << "material:" << materialName <<
350 G4Exception("G4DNAPTBIonisationModel::Sa
351 "Fatal error:: scatteredEner
352 }
353
354 G4double cosTheta = 0.;
355 G4double phi = 0.;
356 RandomizeEjectedElectronDirection(aDynamic
357 cosTheta
358
359 G4double sinTheta = std::sqrt(1. - cosThet
360 G4double dirX = sinTheta * std::cos(phi);
361 G4double dirY = sinTheta * std::sin(phi);
362 G4double dirZ = cosTheta;
363 G4ThreeVector deltaDirection(dirX, dirY, d
364 deltaDirection.rotateUz(primaryDirection);
365
366 // The model is written only for electron
367 // incident electron after each ionization
368 // introduced within this model the follow
369 //
370 // Check if the particle is an electron
371 if (aDynamicParticle->GetDefinition() == G
372 // If yes do the following code until ne
373
374 G4double deltaTotalMomentum =
375 std::sqrt(secondaryKinetic * (secondar
376 G4double finalPx =
377 totalMomentum * primaryDirection.x() -
378 G4double finalPy =
379 totalMomentum * primaryDirection.y() -
380 G4double finalPz =
381 totalMomentum * primaryDirection.z() -
382 G4double finalMomentum = std::sqrt(final
383 finalPx /= finalMomentum;
384 finalPy /= finalMomentum;
385 finalPz /= finalMomentum;
386
387 G4ThreeVector direction(finalPx, finalPy
388 if (direction.unit().getX() > 1 || direc
389 {
390 G4cout << "Fatal error ***************
391 G4cout << "direction problem " << dire
392 G4Exception("G4DNAPTBIonisationModel::
393 "Fatal error:: direction p
394 }
395
396 // Give the new direction to the particl
397 fParticleChangeForGamma->ProposeMomentum
398 }
399 // If the particle is not an electron
400 else {
401 fParticleChangeForGamma->ProposeMomentum
402 }
403
404 // note that secondaryKinetic is the energ
405 G4double scatteredEnergy = k - bindingEner
406
407 if (scatteredEnergy <= 0) {
408 G4cout << "Fatal error *****************
409 G4cout << "k: " << k / eV << G4endl;
410 G4cout << "secondaryKinetic: " << second
411 G4cout << "shell: " << ionizationShell <
412 G4cout << "bindingEnergy: " << bindingEn
413 G4cout << "scatteredEnergy: " << scatter
414 G4cout << "material: " << materialName <
415 G4Exception("G4DNAPTBIonisationModel::Sa
416 "Fatal error:: scatteredEner
417 }
418
419 // Set the new energy of the particle
420 fParticleChangeForGamma->SetProposedKineti
421
422 // Set the energy deposited by the ionizat
423 fParticleChangeForGamma->ProposeLocalEnerg
424
425 // Create the new particle with its charac
426 auto dp = new G4DynamicParticle(G4Electron
427 fvect->push_back(dp);
428
429 // Check if the auger model is activated (
430 if (fpDNAPTBAugerModel) {
431 // run the PTB Auger model
432 if (materialName != "G4_WATER") {
433 fpDNAPTBAugerModel->ComputeAugerEffect
434 }
435 }
436 }
437 }
438
439 //....oooOO0OOooo........oooOO0OOooo........oo
440
441 void G4DNAPTBIonisationModel::ReadDiffCSFile(c
442 c
443 c
444 {
445 // To read and save the informations contain
446
447 // get the path of the G4LEDATA data folder
448 const char* path = G4FindDataDir("G4LEDATA")
449 // if it is not found then quit and print er
450 if (path == nullptr) {
451 G4Exception("G4DNAPTBIonisationModel::Read
452 "G4LEDATA environment variable
453 return;
454 }
455
456 // build the fullFileName path of the data f
457 std::ostringstream fullFileName;
458 fullFileName << path << "/" << file << ".dat
459
460 // open the data file
461 std::ifstream diffCrossSection(fullFileName.
462 // error if file is not there
463 std::stringstream endPath;
464 if (!diffCrossSection) {
465 endPath << "Missing data file: " << file;
466 G4Exception("G4DNAPTBIonisationModel::Init
467 endPath.str().c_str());
468 }
469
470 // load data from the file
471 fTMapWithVec[materialID][p].push_back(0.);
472
473 G4String line;
474
475 // read the file until we reach the end of f
476 // fill fTMapWithVec, diffCrossSectionData,
477 // fEMapWithVector
478 while (std::getline(diffCrossSection, line))
479 // check if the line is comment or empty
480 //
481 std::istringstream testIss(line);
482 G4String test;
483 testIss >> test;
484 // check first caracter to determine if fo
485 if (test == "#") {
486 // skip the line by beginning a new whil
487 continue;
488 }
489 // check if line is empty
490 if (line.empty()) {
491 // skip the line by beginning a new whil
492 continue;
493 }
494 //
495 // end of the check
496
497 // transform the line into a iss
498 std::istringstream iss(line);
499
500 // Initialise the variables to be filled
501 G4double T;
502 G4double E;
503
504 // Filled T and E with the first two numbe
505 iss >> T >> E;
506
507 // Fill the fTMapWithVec container with al
508 // Duplicate must be avoided and this is t
509 if (T != fTMapWithVec[materialID][p].back(
510
511 // iterate on each shell of the correspond
512 for (int shell = 0, eshell = ptbStructure.
513 // map[material][particle][shell][T][E]=
514 // Fill the map with the informations of
515 iss >> diffCrossSectionData[materialID][
516
517 if (fpG4_WATER != nullptr && fpG4_WATER-
518 // map[material][particle][shell][T][C
519 // Fill the map
520 fEnergySecondaryData[materialID][p][sh
521 [diffCrossSectionD
522
523 // map[material][particle][shell][T]=C
524 // Fill the vector within the map
525 fProbaShellMap[materialID][p][shell][T
526 diffCrossSectionData[materialID][p][
527 }
528 else {
529 diffCrossSectionData[materialID][p][sh
530 fEMapWithVector[materialID][p][T].push
531 }
532 }
533 }
534 }
535
536 //....oooOO0OOooo........oooOO0OOooo........oo
537
538 G4double G4DNAPTBIonisationModel::RandomizeEje
539 const G4ParticleDefinition* particleDefiniti
540 {
541 if (particleDefinition == G4Electron::Electr
542 // G4double Tcut=25.0E-6;
543 G4double maximumEnergyTransfer;
544 ((k + ptbStructure.IonisationEnergy(shell,
545 ? maximumEnergyTransfer = k
546 : maximumEnergyTransfer = (k + ptbStruct
547
548 // SI : original method
549 /*
550 G4double crossSectionMaximum = 0.;
551 for(G4double value=waterStructure.Ionisati
552 value+=0.1*eV)
553 {
554 G4double differentialCrossSection = Diff
555 value/eV, shell); if(differentialCrossSect
556 differentialCrossSection;
557 }
558 */
559
560 // SI : alternative method
561
562 // if (k > Tcut)
563 //{
564 G4double crossSectionMaximum = 0.;
565
566 G4double minEnergy = ptbStructure.Ionisati
567 G4double maxEnergy = maximumEnergyTransfer
568 G4int nEnergySteps = 50;
569 G4double value(minEnergy);
570 G4double stpEnergy(std::pow(maxEnergy / va
571 G4int step(nEnergySteps);
572 while (step > 0) {
573 step--;
574 G4double differentialCrossSection =
575 DifferentialCrossSection(particleDefin
576 if (differentialCrossSection >= crossSec
577 crossSectionMaximum = differentialCros
578 value *= stpEnergy;
579 }
580 //
581
582 G4double secondaryElectronKineticEnergy =
583
584 do {
585 secondaryElectronKineticEnergy =
586 G4UniformRand()
587 * (maximumEnergyTransfer - ptbStructur
588
589 } while (
590 G4UniformRand() * crossSectionMaximum >
591 particleDefinition, k / eV,
592 (secondaryElectronKineticEnergy + ptbS
593 shell, materialID));
594
595 return secondaryElectronKineticEnergy;
596 }
597
598 if (particleDefinition == G4Proton::ProtonDe
599 G4double maximumKineticEnergyTransfer = 4.
600
601 G4double crossSectionMaximum = 0.;
602 for (G4double value = ptbStructure.Ionisat
603 value <= 4. * ptbStructure.Ionisation
604 {
605 G4double differentialCrossSection =
606 DifferentialCrossSection(particleDefin
607 if (differentialCrossSection >= crossSec
608 crossSectionMaximum = differentialCros
609 }
610
611 G4double secondaryElectronKineticEnergy =
612 do {
613 secondaryElectronKineticEnergy = G4Unifo
614 } while (
615 G4UniformRand() * crossSectionMaximum >=
616 particleDefinition, k / eV,
617 (secondaryElectronKineticEnergy + ptbS
618 shell, materialID));
619
620 return secondaryElectronKineticEnergy;
621 }
622
623 return 0;
624 }
625
626 //....oooOO0OOooo........oooOO0OOooo........oo
627
628 void G4DNAPTBIonisationModel::RandomizeEjected
629
630
631 {
632 if (p == G4Electron::ElectronDefinition()) {
633 phi = twopi * G4UniformRand();
634 if (secKinetic < 50. * eV)
635 cosTheta = (2. * G4UniformRand()) - 1.;
636 else if (secKinetic <= 200. * eV) {
637 if (G4UniformRand() <= 0.1)
638 cosTheta = (2. * G4UniformRand()) - 1.
639 else
640 cosTheta = G4UniformRand() * (std::sqr
641 }
642 else {
643 G4double sin2O = (1. - secKinetic / k) /
644 cosTheta = std::sqrt(1. - sin2O);
645 }
646 }
647 else if (p == G4Proton::ProtonDefinition())
648 G4double maxSecKinetic = 4. * (electron_ma
649 phi = twopi * G4UniformRand();
650
651 // cosTheta = std::sqrt(secKinetic / maxSe
652
653 // Restriction below 100 eV from Emfietzog
654
655 (secKinetic > 100 * eV) ? cosTheta = std::
656 : cosTheta = (2. *
657 }
658 }
659
660 //....oooOO0OOooo........oooOO0OOooo........oo
661
662 double G4DNAPTBIonisationModel::DifferentialCr
663
664
665
666 {
667 G4double sigma = 0.;
668 G4double shellEnergy(ptbStructure.Ionisation
669 G4double kSE(energyTransfer - shellEnergy);
670
671 if (energyTransfer >= shellEnergy) {
672 G4double valueT1 = 0;
673 G4double valueT2 = 0;
674 G4double valueE21 = 0;
675 G4double valueE22 = 0;
676 G4double valueE12 = 0;
677 G4double valueE11 = 0;
678
679 G4double xs11 = 0;
680 G4double xs12 = 0;
681 G4double xs21 = 0;
682 G4double xs22 = 0;
683
684 if (p == G4Electron::ElectronDefinition())
685 // k should be in eV and energy transfer
686 auto t2 =
687 std::upper_bound(fTMapWithVec[material
688 auto t1 = t2 - 1;
689
690 // SI : the following condition avoids s
691 if (kSE <= fEMapWithVector[materialID][p
692 && kSE <= fEMapWithVector[materialID
693 {
694 auto e12 = std::upper_bound(fEMapWithV
695 fEMapWithV
696 auto e11 = e12 - 1;
697
698 auto e22 = std::upper_bound(fEMapWithV
699 fEMapWithV
700 auto e21 = e22 - 1;
701
702 valueT1 = *t1;
703 valueT2 = *t2;
704 valueE21 = *e21;
705 valueE22 = *e22;
706 valueE12 = *e12;
707 valueE11 = *e11;
708
709 xs11 = diffCrossSectionData[materialID
710 xs12 = diffCrossSectionData[materialID
711 xs21 = diffCrossSectionData[materialID
712 xs22 = diffCrossSectionData[materialID
713 }
714 }
715
716 if (p == G4Proton::ProtonDefinition()) {
717 // k should be in eV and energy transfer
718 auto t2 =
719 std::upper_bound(fTMapWithVec[material
720 auto t1 = t2 - 1;
721
722 auto e12 = std::upper_bound(fEMapWithVec
723 fEMapWithVec
724 auto e11 = e12 - 1;
725
726 auto e22 = std::upper_bound(fEMapWithVec
727 fEMapWithVec
728 auto e21 = e22 - 1;
729
730 valueT1 = *t1;
731 valueT2 = *t2;
732 valueE21 = *e21;
733 valueE22 = *e22;
734 valueE12 = *e12;
735 valueE11 = *e11;
736
737 xs11 = diffCrossSectionData[materialID][
738 xs12 = diffCrossSectionData[materialID][
739 xs21 = diffCrossSectionData[materialID][
740 xs22 = diffCrossSectionData[materialID][
741 }
742
743 G4double xsProduct = xs11 * xs12 * xs21 *
744
745 if (xsProduct != 0.) {
746 sigma = QuadInterpolator(valueE11, value
747 valueT1, valueT
748 }
749 }
750
751 return sigma;
752 }
753
754 //....oooOO0OOooo........oooOO0OOooo........oo
755
756 G4double G4DNAPTBIonisationModel::RandomizeEje
757 const G4ParticleDefinition* p, G4double k, G
758 {
759 // k should be in eV
760
761 // Schematic explanation.
762 // We will do an interpolation to get a fina
763 // 1/ We choose a random number between 0 an
764 // 2/ We look for T_lower and T_upper.
765 // 3/ We look for the cumulated correspondin
766 //
767 // T_low | CS_low_1 -> E_low_1
768 // | CS_low_2 -> E_low_2
769 // T_up | CS_up_1 -> E_up_1
770 // | CS_up_2 -> E_up_2
771 //
772 // 4/ We interpolate to get our E value.
773 //
774 // T_low | CS_low_1 -> E_low_1 -----
775 // | |----> E
776 // | CS_low_2 -> E_low_2 -----
777 //
778 // T_up | CS_up_1 -> E_up_1 -------
779 // | |----> E
780 // | CS_up_2 -> E_up_2 -------
781
782 // Initialize some values
783 //
784 G4double ejectedElectronEnergy = 0.;
785 G4double valueK1 = 0;
786 G4double valueK2 = 0;
787 G4double valueCumulCS21 = 0;
788 G4double valueCumulCS22 = 0;
789 G4double valueCumulCS12 = 0;
790 G4double valueCumulCS11 = 0;
791 G4double secElecE11 = 0;
792 G4double secElecE12 = 0;
793 G4double secElecE21 = 0;
794 G4double secElecE22 = 0;
795 G4String particleName = p->GetParticleName()
796
797 // *****************************************
798 // Get a random number between 0 and 1 to co
799 // *****************************************
800 //
801 // It will allow us to choose an ejected ele
802 G4double random = G4UniformRand();
803
804 // *****************************************
805 // Take the input from the data tables
806 // *****************************************
807
808 // Cumulated tables are like this: T E cumul
809 // We have two sets of loaded data: fTMapWit
810 // energy) and fProbaShellMap which contains
811 // specific T energy value which could not b
812 // values.
813
814 // First, we select the upper and lower T da
815 auto k2 =
816 std::upper_bound(fTMapWithVec[materialID][
817 auto k1 = k2 - 1;
818
819 // Check if we have found a k2 value (0 if w
820 // A missing k2 value can be caused by a ene
821 // Ex : table done for 12*eV -> 1000*eV and
822 // then k2 = 0 and k1 = max of the table.
823 // To detect this, we check that k1 is not s
824 if (*k1 > *k2) {
825 // Error
826 G4cerr << "**************** Fatal error **
827 G4cerr << "G4DNAPTBIonisationModel::Random
828 G4cerr << "You have *k1 > *k2 with k1 " <<
829 G4cerr
830 << "This may be because the energy of th
831 << G4endl;
832 G4cerr << "Particle energy (eV): " << k <<
833 exit(EXIT_FAILURE);
834 }
835
836 // We have a random number and we select the
837 // random number. But we need to do that for
838 //
839 // First one.
840 auto cumulCS12 =
841 std::upper_bound(fProbaShellMap[materialID
842 fProbaShellMap[materialID
843 auto cumulCS11 = cumulCS12 - 1;
844 // Second one.
845 auto cumulCS22 =
846 std::upper_bound(fProbaShellMap[materialID
847 fProbaShellMap[materialID
848 auto cumulCS21 = cumulCS22 - 1;
849
850 // Now that we have the "values" through poi
851 valueK1 = *k1;
852 valueK2 = *k2;
853 valueCumulCS11 = *cumulCS11;
854 valueCumulCS12 = *cumulCS12;
855 valueCumulCS21 = *cumulCS21;
856 valueCumulCS22 = *cumulCS22;
857
858 // *****************************************
859 // Do the interpolation to get the ejected e
860 // *****************************************
861
862 // Here we will get four E values correspond
863 // previously selected. But we need to take
864 // by using the ionisation cross section tab
865 // differential cross sections (or cumulated
866 // When looking for the cumulated cross sect
867 // (for the lower T), the upper_bound method
868 // method will return the last E value prese
869 // being the highest, we will later perform
870 // T) and a small E value (for the upper T).
871 // equal to zero for the lower T then it mea
872 // secondary electron. But, in our situation
873 // interpolate T value between Tupper Tlower
874 // between 0 and the E value (upper T).
875 //
876 if (cumulCS12 == fProbaShellMap[materialID][
877 // Here we are in the special case and we
878 // interpolation.
879 secElecE11 = 0;
880 secElecE12 = 0;
881 secElecE21 = fEnergySecondaryData[material
882 secElecE22 = fEnergySecondaryData[material
883
884 valueCumulCS11 = 0;
885 valueCumulCS12 = 0;
886 }
887 else {
888 // No special case, interpolation will hap
889 secElecE11 = fEnergySecondaryData[material
890 secElecE12 = fEnergySecondaryData[material
891 secElecE21 = fEnergySecondaryData[material
892 secElecE22 = fEnergySecondaryData[material
893 }
894
895 ejectedElectronEnergy =
896 QuadInterpolator(valueCumulCS11, valueCumu
897 secElecE12, secElecE21, s
898
899 // *****************************************
900 // Some tests for debugging
901 // *****************************************
902
903 G4double bindingEnergy(ptbStructure.Ionisati
904 if (k - ejectedElectronEnergy - bindingEnerg
905 G4cout << "k " << k << G4endl;
906 G4cout << "material ID : " << materialID <
907 G4cout << "secondaryKin " << ejectedElectr
908 G4cout << "shell " << ionizationLevelIndex
909 G4cout << "bindingEnergy " << bindingEnerg
910 G4cout << "scatteredEnergy " << k - ejecte
911 G4cout << "rand " << random << G4endl;
912 G4cout << "surrounding k values: valueK1 v
913 G4cout << "surrounding E values: secElecE1
914 << secElecE11 << " " << secElecE12
915 << G4endl;
916 G4cout
917 << "surrounding cumulCS values: valueCum
918 << valueCumulCS11 << " " << valueCumulCS
919 << " " << G4endl;
920 G4ExceptionDescription errmsg;
921 errmsg << "*****************************"
922 errmsg << "Fatal error, EXIT." << G4endl;
923 G4Exception("G4DNAPTBIonisationModel::Rand
924 FatalException, errmsg);
925 exit(EXIT_FAILURE);
926 }
927
928 return ejectedElectronEnergy * eV;
929 }
930
931 //....oooOO0OOooo........oooOO0OOooo........oo
932
933 G4double G4DNAPTBIonisationModel::LogLogInterp
934
935 {
936 G4double value(0);
937
938 // Switch to log-lin interpolation for faste
939
940 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0)
941 G4double d1 = std::log10(xs1);
942 G4double d2 = std::log10(xs2);
943 value = std::pow(10., (d1 + (d2 - d1) * (e
944 }
945
946 // Switch to lin-lin interpolation for faste
947 // in case one of xs1 or xs2 (=cum proba) va
948
949 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0)
950 G4double d1 = xs1;
951 G4double d2 = xs2;
952 value = (d1 + (d2 - d1) * (e - e1) / (e2 -
953 }
954
955 return value;
956 }
957
958 //....oooOO0OOooo........oooOO0OOooo........oo
959
960 G4double G4DNAPTBIonisationModel::QuadInterpol
961
962
963
964 {
965 G4double interpolatedvalue1, interpolatedval
966 (xs11 != xs12) ? interpolatedvalue1 = LogLog
967 : interpolatedvalue1 = xs11;
968
969 (xs21 != xs22) ? interpolatedvalue2 = LogLog
970 : interpolatedvalue2 = xs21;
971
972 (interpolatedvalue1 != interpolatedvalue2)
973 ? value = LogLogInterpolate(t1, t2, t, int
974 : value = interpolatedvalue1;
975 return value;
976 }
977