| 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 // CPA100 ionisation model class for electrons
27 //
28 // Based on the work of M. Terrissol and M. C.
29 //
30 // Users are requested to cite the following p
31 // - M. Terrissol, A. Baudre, Radiat. Prot. Do
32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terr
33 // M. Bardies, N. Lampe, S. Incerti, Phys. M
34 //
35 // Authors of this class:
36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bor
37 //
38 // 15.01.2014: creation
39 //
40 // Based on the study by S. Zein et. al. Nucl.
41 // 1/2/2023 : Hoang added modification
42
43 #include "G4DNACPA100IonisationModel.hh"
44
45 #include "G4DNAChemistryManager.hh"
46 #include "G4DNAMaterialManager.hh"
47 #include "G4DNAMolecularMaterial.hh"
48 #include "G4EnvironmentUtils.hh"
49 #include "G4LossTableManager.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4UAtomicDeexcitation.hh"
53
54 #include <fstream>
55
56 //....oooOO0OOooo........oooOO0OOooo........oo
57
58 using namespace std;
59
60 //....oooOO0OOooo........oooOO0OOooo........oo
61
62 G4DNACPA100IonisationModel::G4DNACPA100Ionisat
63
64 : G4VDNAModel(nam, "all")
65 {
66 fpGuanine = G4Material::GetMaterial("G4_GUAN
67 fpG4_WATER = G4Material::GetMaterial("G4_WAT
68 fpDeoxyribose = G4Material::GetMaterial("G4_
69 fpCytosine = G4Material::GetMaterial("G4_CYT
70 fpThymine = G4Material::GetMaterial("G4_THYM
71 fpAdenine = G4Material::GetMaterial("G4_ADEN
72 fpPhosphate = G4Material::GetMaterial("G4_PH
73 fpParticle = G4Electron::ElectronDefinition(
74 }
75
76 //....oooOO0OOooo........oooOO0OOooo........oo
77
78 void G4DNACPA100IonisationModel::Initialise(co
79 co
80 {
81 if (isInitialised) {
82 return;
83 }
84 if (verboseLevel > 3) {
85 G4cout << "Calling G4DNACPA100IonisationMo
86 }
87
88 if (!G4DNAMaterialManager::Instance()->IsLoc
89 if (p != fpParticle) {
90 std::ostringstream oss;
91 oss << " Model is not applied for this p
92 G4Exception("G4DNACPA100IonisationModel:
93 FatalException, oss.str().c_
94 }
95
96 const char* path = G4FindDataDir("G4LEDATA
97
98 if (path == nullptr) {
99 G4Exception("G4DNACPA100IonisationModel:
100 "G4LEDATA environment variab
101 return;
102 }
103
104 std::size_t index;
105 if (fpG4_WATER != nullptr) {
106 index = fpG4_WATER->GetIndex();
107 G4String eFullFileName = "";
108 fasterCode ? eFullFileName = "/dna/sigma
109 : eFullFileName = "/dna/sigma
110 AddCrossSectionData(index, p, "dna/sigma
111 1.e-20 * m * m);
112 SetLowELimit(index, p, 11 * eV);
113 SetHighELimit(index, p, 255955 * eV);
114 }
115 if (fpGuanine != nullptr) {
116 index = fpGuanine->GetIndex();
117 G4String eFullFileName = "";
118 if(useDcs) {
119 fasterCode ? eFullFileName = "/dna/sig
120 : eFullFileName = "/dna/sig
121 }
122 AddCrossSectionData(index, p, "dna/sigma
123 1. * cm * cm);
124 SetLowELimit(index, p, 11 * eV);
125 SetHighELimit(index, p, 1 * MeV);
126 }
127 if (fpDeoxyribose != nullptr) {
128 index = fpDeoxyribose->GetIndex();
129 G4String eFullFileName = "";
130 if(useDcs) {
131 eFullFileName = "/dna/sigmadiff_cumula
132 }
133 AddCrossSectionData(index, p, "dna/sigma
134 1. * cm * cm);
135 SetLowELimit(index, p, 11 * eV);
136 SetHighELimit(index, p, 1 * MeV);
137 }
138 if (fpCytosine != nullptr) {
139 index = fpCytosine->GetIndex();
140 G4String eFullFileName = "";
141 if(useDcs) {
142 fasterCode ? eFullFileName = "/dna/sig
143 : eFullFileName = "/dna/sig
144 }
145 AddCrossSectionData(index, p, "dna/sigma
146 1. * cm * cm);
147 SetLowELimit(index, p, 11 * eV);
148 SetHighELimit(index, p, 1 * MeV);
149 }
150 if (fpThymine != nullptr) {
151 index = fpThymine->GetIndex();
152 G4String eFullFileName = "";
153 if(useDcs) {
154 fasterCode ? eFullFileName = "/dna/sig
155 : eFullFileName = "/dna/sig
156 }
157 AddCrossSectionData(index, p, "dna/sigma
158 1. * cm * cm);
159 SetLowELimit(index, p, 11 * eV);
160 SetHighELimit(index, p, 1 * MeV);
161 }
162 if (fpAdenine != nullptr) {
163 index = fpAdenine->GetIndex();
164 G4String eFullFileName = "";
165 if(useDcs) {
166 fasterCode ? eFullFileName = "/dna/sig
167 : eFullFileName = "/dna/sig
168 }
169 AddCrossSectionData(index, p, "dna/sigma
170 1. * cm * cm);
171 SetLowELimit(index, p, 11 * eV);
172 SetHighELimit(index, p, 1 * MeV);
173 }
174 if (fpPhosphate != nullptr) {
175 index = fpPhosphate->GetIndex();
176 G4String eFullFileName = "";
177 if(useDcs) {
178 eFullFileName = "dna/sigmadiff_cumulat
179 }
180 AddCrossSectionData(index, p, "dna/sigma
181 1. * cm * cm);
182 SetLowELimit(index, p, 11 * eV);
183 SetHighELimit(index, p, 1 * MeV);
184 }
185 LoadCrossSectionData(p);
186 G4DNAMaterialManager::Instance()->SetMaste
187 fpModelData = this;
188 }
189 else {
190 auto dataModel = dynamic_cast<G4DNACPA100I
191 G4DNAMaterialManager::Instance()->GetMod
192 if (dataModel == nullptr) {
193 G4cout << "G4DNACPA100IonisationModel::C
194 throw;
195 }
196 fpModelData = dataModel;
197 }
198
199 fAtomDeexcitation = G4LossTableManager::Inst
200
201 fParticleChangeForGamma = GetParticleChangeF
202 isInitialised = true;
203 }
204
205 G4double G4DNACPA100IonisationModel::CrossSect
206
207
208 {
209 // initialise the cross section value (outpu
210 G4double sigma(0);
211
212 // Get the current particle name
213 const G4String& particleName = p->GetParticl
214
215 if (p != fpParticle) {
216 G4Exception("G4DNACPA100IonisationModel::C
217 "No model is registered for th
218 }
219
220 auto matID = material->GetIndex();
221
222 // Set the low and high energy limits
223 G4double lowLim = fpModelData->GetLowELimit(
224 G4double highLim = fpModelData->GetHighELimi
225
226 // Check that we are in the correct energy r
227 if (ekin >= lowLim && ekin < highLim) {
228 // Get the map with all the model data tab
229 auto tableData = fpModelData->GetData();
230
231 if ((*tableData)[matID][p] == nullptr) {
232 G4Exception("G4DNACPA100IonisationModel:
233 "No model is registered");
234 }
235 else {
236 sigma = (*tableData)[matID][p]->FindValu
237 }
238
239 if (verboseLevel > 2) {
240 auto MolDensity =
241 (*G4DNAMolecularMaterial::Instance()->
242 G4cout << "_____________________________
243 G4cout << "°°° G4DNACPA100IonisationM
244 G4cout << "°°° Kinetic energy(eV)=" <
245 G4cout << "°°° lowLim (eV) = " << low
246 G4cout << "°°° Materials = " << (*G4M
247 G4cout << "°°° Cross section per " <<
248 << G4endl;
249 G4cout << "°°° Cross section per Phos
250 << sigma * MolDensity / (1. / cm)
251 G4cout << "°°° G4DNACPA100IonisationM
252 }
253 }
254
255 auto MolDensity = (*G4DNAMolecularMaterial::
256 return sigma * MolDensity;
257 }
258
259 //....oooOO0OOooo........oooOO0OOooo........oo
260
261 void G4DNACPA100IonisationModel::SampleSeconda
262 std::vector<G4DynamicParticle*>* fvect,
263 const G4MaterialCutsCouple* couple, // must
264 const G4DynamicParticle* particle, G4double,
265 {
266 if (verboseLevel > 3) {
267 G4cout << "Calling SampleSecondaries() of
268 }
269 auto k = particle->GetKineticEnergy();
270
271 const G4Material* material = couple->GetMate
272
273 auto MatID = material->GetIndex();
274
275 auto p = particle->GetDefinition();
276
277 auto lowLim = fpModelData->GetLowELimit(MatI
278 auto highLim = fpModelData->GetHighELimit(Ma
279
280 // Check if we are in the correct energy ran
281 if (k >= lowLim && k < highLim) {
282 const auto& primaryDirection = particle->G
283 auto particleMass = particle->GetDefinitio
284 auto totalEnergy = k + particleMass;
285 auto pSquare = k * (totalEnergy + particle
286 auto totalMomentum = std::sqrt(pSquare);
287 G4int shell = -1;
288 G4double bindingEnergy, secondaryKinetic;
289 shell = fpModelData->RandomSelectShell(k,
290 bindingEnergy = iStructure.IonisationEnerg
291
292 if (k < bindingEnergy) {
293 return;
294 }
295
296 auto info = std::make_tuple(MatID, k, shel
297
298 secondaryKinetic = -1000 * eV;
299 if (fpG4_WATER->GetIndex() != MatID) {//fo
300 secondaryKinetic = fpModelData->Randomiz
301 }else if(fasterCode){
302 secondaryKinetic = fpModelData->Random
303 }else{
304 secondaryKinetic = fpModelData->Random
305 }
306
307 G4double cosTheta = 0.;
308 G4double phi = 0.;
309 RandomizeEjectedElectronDirection(particle
310 phi);
311
312 G4double sinTheta = std::sqrt(1. - cosThet
313 G4double dirX = sinTheta * std::cos(phi);
314 G4double dirY = sinTheta * std::sin(phi);
315 G4double dirZ = cosTheta;
316 G4ThreeVector deltaDirection(dirX, dirY, d
317 deltaDirection.rotateUz(primaryDirection);
318
319 // SI - For atom. deexc. tagging - 23/05/2
320 if (secondaryKinetic > 0) {
321 auto dp = new G4DynamicParticle(G4Electr
322 fvect->push_back(dp);
323 }
324
325 if (particle->GetDefinition() != fpParticl
326 fParticleChangeForGamma->ProposeMomentum
327 }
328 else {
329 G4double deltaTotalMomentum =
330 std::sqrt(secondaryKinetic * (secondar
331 G4double finalPx =
332 totalMomentum * primaryDirection.x() -
333 G4double finalPy =
334 totalMomentum * primaryDirection.y() -
335 G4double finalPz =
336 totalMomentum * primaryDirection.z() -
337 G4double finalMomentum = std::sqrt(final
338 finalPx /= finalMomentum;
339 finalPy /= finalMomentum;
340 finalPz /= finalMomentum;
341
342 G4ThreeVector direction;
343 direction.set(finalPx, finalPy, finalPz)
344
345 fParticleChangeForGamma->ProposeMomentum
346 }
347
348 // SI - For atom. deexc. tagging - 23/05/2
349
350 // AM: sample deexcitation
351 // here we assume that H_{2}O electronic l
352 // this can be considered true with a roug
353
354 G4double scatteredEnergy = k - bindingEner
355
356 // SI: only atomic deexcitation from K she
357 // Hoang: only for water
358 if (material == G4Material::GetMaterial("G
359 std::size_t secNumberInit = 0; // need
360 std::size_t secNumberFinal = 0; // So I
361 if ((fAtomDeexcitation != nullptr) && sh
362 G4int Z = 8;
363 auto Kshell = fAtomDeexcitation->GetAt
364 secNumberInit = fvect->size();
365 fAtomDeexcitation->GenerateParticles(f
366 secNumberFinal = fvect->size();
367 if (secNumberFinal > secNumberInit) {
368 for (std::size_t i = secNumberInit;
369 // Check if there is enough residu
370 if (bindingEnergy >= ((*fvect)[i])
371 // Ok, this is a valid secondary
372 bindingEnergy -= ((*fvect)[i])->
373 }
374 else {
375 // Invalid secondary: not enough
376 // Keep its energy in the local
377 delete (*fvect)[i];
378 (*fvect)[i] = nullptr;
379 }
380 }
381 }
382 }
383 }
384
385 // This should never happen
386 if (bindingEnergy < 0.0) {
387 G4Exception("G4DNACPA100IonisatioModel1:
388 "Negative local energy depos
389 }
390 if (!statCode) {
391 fParticleChangeForGamma->SetProposedKine
392 fParticleChangeForGamma->ProposeLocalEne
393 }
394 else {
395 fParticleChangeForGamma->SetProposedKine
396 fParticleChangeForGamma->ProposeLocalEne
397 }
398
399 // only water for chemistry
400 if (fpG4_WATER != nullptr && material == G
401 const G4Track* theIncomingTrack = fParti
402 G4DNAChemistryManager::Instance()->Creat
403
404 }
405 }
406 }
407
408 //....oooOO0OOooo........oooOO0OOooo........oo
409
410 G4double G4DNACPA100IonisationModel::Randomize
411 {
412 auto MatID = std::get<0>(info);
413 auto k = std::get<1>(info);
414 auto shell = std::get<2>(info);
415 G4double maximumEnergyTransfer = 0.;
416 auto IonLevel = iStructure.IonisationEnergy(
417 (k + IonLevel) / 2. > k ? maximumEnergyTrans
418
419 G4double crossSectionMaximum = 0.;
420
421 G4double minEnergy = IonLevel;
422 G4double maxEnergy = maximumEnergyTransfer;
423
424 // nEnergySteps can be optimized - 100 by de
425 G4int nEnergySteps = 50;
426
427 G4double value(minEnergy);
428 G4double stpEnergy(std::pow(maxEnergy / valu
429 G4int step(nEnergySteps);
430 G4double differentialCrossSection = 0.;
431 while (step > 0) {
432 step--;
433 differentialCrossSection = DifferentialCro
434
435 if (differentialCrossSection > 0) {
436 crossSectionMaximum = differentialCrossS
437 break;
438 }
439 value *= stpEnergy;
440 }
441
442 G4double secondaryElectronKineticEnergy = 0.
443 do {
444 secondaryElectronKineticEnergy = G4Uniform
445 } while (G4UniformRand() * crossSectionMaxim
446 > DifferentialCrossSection(info, (s
447
448 return secondaryElectronKineticEnergy;
449 }
450
451 //....oooOO0OOooo........oooOO0OOooo........oo
452
453 void G4DNACPA100IonisationModel::RandomizeEjec
454
455
456
457 {
458 phi = twopi * G4UniformRand();
459 G4double sin2O = (1. - secKinetic / k) / (1.
460 cosTheta = std::sqrt(1. - sin2O);
461 }
462
463 //....oooOO0OOooo........oooOO0OOooo........oo
464
465 G4double G4DNACPA100IonisationModel::Different
466
467 {
468 auto MatID = std::get<0>(info);
469 auto k = std::get<1>(info) / eV; // in eV u
470 auto shell = std::get<2>(info);
471 G4double sigma = 0.;
472 G4double shellEnergy = iStructure.Ionisation
473 G4double kSE(energyTransfer - shellEnergy);
474
475 if (energyTransfer >= shellEnergy) {
476 G4double valueT1 = 0;
477 G4double valueT2 = 0;
478 G4double valueE21 = 0;
479 G4double valueE22 = 0;
480 G4double valueE12 = 0;
481 G4double valueE11 = 0;
482
483 G4double xs11 = 0;
484 G4double xs12 = 0;
485 G4double xs21 = 0;
486 G4double xs22 = 0;
487
488 auto t2 = std::upper_bound(fTMapWithVec[Ma
489 fTMapWithVec[Ma
490 auto t1 = t2 - 1;
491
492 if (kSE <= fEMapWithVector[MatID][fpPartic
493 && kSE <= fEMapWithVector[MatID][fpPar
494 {
495 auto e12 = std::upper_bound(fEMapWithVec
496 fEMapWithVec
497 auto e11 = e12 - 1;
498
499 auto e22 = std::upper_bound(fEMapWithVec
500 fEMapWithVec
501 auto e21 = e22 - 1;
502
503 valueT1 = *t1;
504 valueT2 = *t2;
505 valueE21 = *e21;
506 valueE22 = *e22;
507 valueE12 = *e12;
508 valueE11 = *e11;
509
510 xs11 = diffCrossSectionData[MatID][fpPar
511 xs12 = diffCrossSectionData[MatID][fpPar
512 xs21 = diffCrossSectionData[MatID][fpPar
513 xs22 = diffCrossSectionData[MatID][fpPar
514 }
515
516 G4double xsProduct = xs11 * xs12 * xs21 *
517
518 if (xsProduct != 0.) {
519 sigma = QuadInterpolator(valueE11, value
520 valueT1, valueT
521 }
522 }
523
524 return sigma;
525 }
526
527 //....oooOO0OOooo........oooOO0OOooo........oo
528
529 G4double G4DNACPA100IonisationModel::Interpola
530
531 {
532 G4double value = 0.;
533
534 // Log-log interpolation by default
535
536 if (e1 != 0 && e2 != 0 && (std::log10(e2) -
537 G4double a = (std::log10(xs2) - std::log10
538 G4double b = std::log10(xs2) - a * std::lo
539 G4double sigma = a * std::log10(e) + b;
540 value = (std::pow(10., sigma));
541 }
542
543 // Switch to lin-lin interpolation
544 /*
545 if ((e2-e1)!=0)
546 {
547 G4double d1 = xs1;
548 G4double d2 = xs2;
549 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1
550 }
551 */
552
553 // Switch to log-lin interpolation for faste
554
555 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 &
556 G4double d1 = std::log10(xs1);
557 G4double d2 = std::log10(xs2);
558 value = std::pow(10., (d1 + (d2 - d1) * (e
559 }
560
561 // Switch to lin-lin interpolation for faste
562 // in case one of xs1 or xs2 (=cum proba) va
563
564 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0)
565 G4double d1 = xs1;
566 G4double d2 = xs2;
567 value = (d1 + (d2 - d1) * (e - e1) / (e2 -
568 }
569 return value;
570 }
571
572 //....oooOO0OOooo........oooOO0OOooo........oo
573
574 G4double G4DNACPA100IonisationModel::QuadInter
575
576
577
578 {
579 G4double interpolatedvalue1 = Interpolate(e1
580 G4double interpolatedvalue2 = Interpolate(e2
581 G4double value = Interpolate(t1, t2, t, inte
582
583 return value;
584 }
585
586 G4double
587 G4DNACPA100IonisationModel::RandomizeEjectedEl
588 {
589 auto MatID = std::get<0>(info);
590 auto shell = std::get<2>(info);
591 G4double secondaryElectronKineticEnergy =
592 RandomTransferedEnergy(info) * eV - iStruc
593 if (secondaryElectronKineticEnergy < 0.) {
594 return 0.;
595 }
596 return secondaryElectronKineticEnergy;
597 }
598
599 G4double G4DNACPA100IonisationModel::RandomTra
600 {
601 auto materialID = std::get<0>(info);
602 auto k = std::get<1>(info) / eV; // data ta
603 auto shell = std::get<2>(info);
604 G4double ejectedElectronEnergy = 0.;
605 G4double valueK1 = 0;
606 G4double valueK2 = 0;
607 G4double valueCumulCS21 = 0;
608 G4double valueCumulCS22 = 0;
609 G4double valueCumulCS12 = 0;
610 G4double valueCumulCS11 = 0;
611 G4double secElecE11 = 0;
612 G4double secElecE12 = 0;
613 G4double secElecE21 = 0;
614 G4double secElecE22 = 0;
615
616 if (k == fTMapWithVec[materialID][fpParticle
617 k = k * (1. - 1e-12);
618 }
619
620 G4double random = G4UniformRand();
621 auto k2 = std::upper_bound(fTMapWithVec[mate
622 fTMapWithVec[mate
623 auto k1 = k2 - 1;
624
625 if (random <= fProbaShellMap[materialID][fpP
626 && random <= fProbaShellMap[materialID][
627 {
628 auto cumulCS12 =
629 std::upper_bound(fProbaShellMap[material
630 fProbaShellMap[material
631 auto cumulCS11 = cumulCS12 - 1;
632 // Second one.
633 auto cumulCS22 =
634 std::upper_bound(fProbaShellMap[material
635 fProbaShellMap[material
636 auto cumulCS21 = cumulCS22 - 1;
637
638 valueK1 = *k1;
639 valueK2 = *k2;
640 valueCumulCS11 = *cumulCS11;
641 valueCumulCS12 = *cumulCS12;
642 valueCumulCS21 = *cumulCS21;
643 valueCumulCS22 = *cumulCS22;
644
645 secElecE11 = fEnergySecondaryData[material
646 secElecE12 = fEnergySecondaryData[material
647 secElecE21 = fEnergySecondaryData[material
648 secElecE22 = fEnergySecondaryData[material
649
650 if (valueCumulCS11 == 0. && valueCumulCS12
651 auto interpolatedvalue2 =
652 Interpolate(valueCumulCS21, valueCumul
653 G4double valueNrjTransf = Interpolate(va
654 return valueNrjTransf;
655 }
656 }
657
658 if (random > fProbaShellMap[materialID][fpPa
659 auto cumulCS22 =
660 std::upper_bound(fProbaShellMap[material
661 fProbaShellMap[material
662 auto cumulCS21 = cumulCS22 - 1;
663 valueK1 = *k1;
664 valueK2 = *k2;
665 valueCumulCS21 = *cumulCS21;
666 valueCumulCS22 = *cumulCS22;
667
668 secElecE21 = fEnergySecondaryData[material
669 secElecE22 = fEnergySecondaryData[material
670
671 G4double interpolatedvalue2 =
672 Interpolate(valueCumulCS21, valueCumulCS
673
674 G4double value = Interpolate(valueK1, valu
675 return value;
676 }
677 G4double nrjTransfProduct = secElecE11 * sec
678
679 if (nrjTransfProduct != 0.) {
680 ejectedElectronEnergy =
681 QuadInterpolator(valueCumulCS11, valueCu
682 secElecE12, secElecE21,
683 }
684 return ejectedElectronEnergy;
685 }
686
687 //....oooOO0OOooo........oooOO0OOooo........oo
688
689 G4double
690 G4DNACPA100IonisationModel::RandomizeEjectedEl
691 {
692 auto MatID = std::get<0>(info);
693 auto tt = std::get<1>(info);
694 auto shell = std::get<2>(info);
695 // ***** METHOD by M. C. Bordage ***** (opti
696 // Composition sampling method based on eq
697
698 //// Defining constants
699 G4double alfa = 1. / 137; // fine structure
700 G4double e_charge = 1.6e-19; // electron ch
701 G4double e_mass = 9.1e-31; // electron mass
702 G4double c = 3e8; // speed of light in vacu
703 G4double mc2 = e_mass * c * c / e_charge; /
704
705 G4double BB = iStructure.IonisationEnergy(sh
706
707 if (tt <= BB) return 0.;
708
709 G4double b_prime = BB / mc2; // binding ene
710 G4double beta_b2 = 1. - 1. / ((1 + b_prime)
711
712 //// Indicent energy
713 //// tt is the incident electron energy
714
715 G4double t_prime = tt / mc2; // incident en
716 G4double t = tt / BB; // reduced incident e
717
718 G4double D = (1 + 2 * t_prime) / ((1 + t_pri
719 G4double F = b_prime * b_prime / ((1 + t_pri
720
721 G4double beta_t2 = 1 - 1 / ((1 + t_prime) *
722
723 G4double PHI_R = std::cos(std::sqrt(alfa * a
724 * std::log(beta_t2
725 G4double G_R = std::log(beta_t2 / (1 - beta_
726
727 G4double tplus1 = t + 1;
728 G4double tminus1 = t - 1;
729 G4double tplus12 = tplus1 * tplus1;
730 G4double ZH1max = 1 + F - (PHI_R * D * (2 *
731 G4double ZH2max = 1 - PHI_R * D / 4;
732
733 G4double A1_p = ZH1max * tminus1 / tplus1;
734 G4double A2_p = ZH2max * tminus1 / (t * tplu
735 G4double A3_p = ((tplus12 - 4) / tplus12) *
736
737 G4double AAA = A1_p + A2_p + A3_p;
738
739 G4double AA1_R = A1_p / AAA;
740 G4double AA2_R = (A1_p + A2_p) / AAA;
741
742 G4int FF = 0;
743 G4double fx = 0;
744 G4double gx = 0;
745 G4double gg = 0;
746 G4double wx = 0;
747
748 G4double r1 = 0;
749 G4double r2 = 0;
750 G4double r3 = 0;
751
752 //
753
754 do {
755 r1 = G4UniformRand();
756 r2 = G4UniformRand();
757 r3 = G4UniformRand();
758
759 if (r1 > AA2_R)
760 FF = 3;
761 else if ((r1 > AA1_R) && (r1 < AA2_R))
762 FF = 2;
763 else
764 FF = 1;
765
766 switch (FF) {
767 case 1: {
768 fx = r2 * tminus1 / tplus1;
769 wx = 1 / (1 - fx) - 1;
770 gg = PHI_R * D * (wx + 1) / tplus1;
771 gx = 1 - gg;
772 gx = gx - gg * (wx + 1) / (2 * (t - wx
773 gx = gx + F * (wx + 1) * (wx + 1);
774 gx = gx / ZH1max;
775 break;
776 }
777
778 case 2: {
779 fx = tplus1 + r2 * tminus1;
780 wx = t * tminus1 * r2 / fx;
781 gx = 1 - (PHI_R * D * (t - wx) / (2 *
782 gx = gx / ZH2max;
783 break;
784 }
785
786 case 3: {
787 fx = 1 - r2 * (tplus12 - 4) / tplus12;
788 wx = std::sqrt(1 / fx) - 1;
789 gg = (wx + 1) / (t - wx);
790 gx = (1 + gg * gg * gg) / 2;
791 break;
792 }
793 } // switch
794
795 } while (r3 > gx);
796
797 return wx * BB;
798 }
799
800 //....oooOO0OOooo........oooOO0OOooo........oo
801
802 void G4DNACPA100IonisationModel::ReadDiffCSFil
803
804
805 {
806 const char* path = G4FindDataDir("G4LEDATA")
807 if (path == nullptr) {
808 G4Exception("G4DNACPA100IonisationModel::R
809 "G4LEDATA environment variable
810 return;
811 }
812
813 std::ostringstream fullFileName;
814 fullFileName << path << "/" << file << ".dat
815
816 std::ifstream diffCrossSection(fullFileName.
817 std::stringstream endPath;
818 if (!diffCrossSection) {
819 endPath << "Missing data file: " << file;
820 G4Exception("G4DNACPA100IonisationModel::I
821 endPath.str().c_str());
822 }
823
824 // load data from the file
825 fTMapWithVec[materialID][p].push_back(0.);
826
827 G4String line;
828
829 while (!diffCrossSection.eof()) {
830 G4double T, E;
831 diffCrossSection >> T >> E;
832
833 if (T != fTMapWithVec[materialID][p].back(
834 fTMapWithVec[materialID][p].push_back(T)
835 }
836
837 // T is incident energy, E is the energy t
838 if (T != fTMapWithVec[materialID][p].back(
839 fTMapWithVec[materialID][p].push_back(T)
840 }
841
842 auto eshell = (G4int)iStructure.NumberOfLe
843 for (G4int shell = 0; shell < eshell; ++sh
844 diffCrossSection >> diffCrossSectionData
845 if (fasterCode) {
846 fEnergySecondaryData[materialID][p][sh
847 [diffCrossSectionD
848
849 fProbaShellMap[materialID][p][shell][T
850 diffCrossSectionData[materialID][p][
851 }
852 else {
853 diffCrossSectionData[materialID][p][sh
854 fEMapWithVector[materialID][p][T].push
855 }
856 }
857 }
858 }
859