| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // CPA100 ionisation model class for electrons
27 //
28 // Based on the work of M. Terrissol and M. C. Bordage
29 //
30 // Users are requested to cite the following papers:
31 // - M. Terrissol, A. Baudre, Radiat. Prot. Dosim. 31 (1990) 175-177
32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terrissol, X. Franceries,
33 // M. Bardies, N. Lampe, S. Incerti, Phys. Med. 32 (2016) 1833-1840
34 //
35 // Authors of this class:
36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bordes, S. Incerti
37 //
38 // 15.01.2014: creation
39 //
40 // Based on the study by S. Zein et. al. Nucl. Inst. Meth. B 488 (2021) 70-82
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........oooOO0OOooo........oooOO0OOooo....
57
58 using namespace std;
59
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61
62 G4DNACPA100IonisationModel::G4DNACPA100IonisationModel(const G4ParticleDefinition*,
63 const G4String& nam)
64 : G4VDNAModel(nam, "all")
65 {
66 fpGuanine = G4Material::GetMaterial("G4_GUANINE", false);
67 fpG4_WATER = G4Material::GetMaterial("G4_WATER", false);
68 fpDeoxyribose = G4Material::GetMaterial("G4_DEOXYRIBOSE", false);
69 fpCytosine = G4Material::GetMaterial("G4_CYTOSINE", false);
70 fpThymine = G4Material::GetMaterial("G4_THYMINE", false);
71 fpAdenine = G4Material::GetMaterial("G4_ADENINE", false);
72 fpPhosphate = G4Material::GetMaterial("G4_PHOSPHORIC_ACID", false);
73 fpParticle = G4Electron::ElectronDefinition();
74 }
75
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78 void G4DNACPA100IonisationModel::Initialise(const G4ParticleDefinition* p,
79 const G4DataVector& /*cuts*/)
80 {
81 if (isInitialised) {
82 return;
83 }
84 if (verboseLevel > 3) {
85 G4cout << "Calling G4DNACPA100IonisationModel::Initialise()" << G4endl;
86 }
87
88 if (!G4DNAMaterialManager::Instance()->IsLocked()) {
89 if (p != fpParticle) {
90 std::ostringstream oss;
91 oss << " Model is not applied for this particle " << p->GetParticleName();
92 G4Exception("G4DNACPA100IonisationModel::G4DNACPA100IonisationModel", "CPA001",
93 FatalException, oss.str().c_str());
94 }
95
96 const char* path = G4FindDataDir("G4LEDATA");
97
98 if (path == nullptr) {
99 G4Exception("G4DNACPA100IonisationModel::Initialise", "em0006", FatalException,
100 "G4LEDATA environment variable not set.");
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/sigmadiff_cumulated_ionisation_e_cpa100_rel"
109 : eFullFileName = "/dna/sigmadiff_ionisation_e_cpa100_rel";
110 AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_form_rel", eFullFileName,
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/sigmadiff_cumulated_elastic_e_cpa100_guanine"
120 : eFullFileName = "/dna/sigmadiff_ionisation_e_cpa100_guanine";
121 }
122 AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_guanine", eFullFileName,
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_cumulated_ionisation_e_cpa100_deoxyribose";
132 }
133 AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_deoxyribose", eFullFileName,
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/sigmadiff_cumulated_ionisation_e_cpa100_cytosine"
143 : eFullFileName = "/dna/sigmadiff_ionisation_e_cpa100_cytosine";
144 }
145 AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_cytosine", eFullFileName,
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/sigmadiff_cumulated_ionisation_e_cpa100_thymine"
155 : eFullFileName = "/dna/sigmadiff_ionisation_e_cpa100_thymine";
156 }
157 AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_thymine", eFullFileName,
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/sigmadiff_cumulated_ionisation_e_cpa100_adenine"
167 : eFullFileName = "/dna/sigmadiff_ionisation_e_cpa100_adenine";
168 }
169 AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_adenine", eFullFileName,
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_cumulated_ionisation_e_cpa100_phosphoric_acid";
179 }
180 AddCrossSectionData(index, p, "dna/sigma_ionisation_e_cpa100_phosphoric_acid",eFullFileName,
181 1. * cm * cm);
182 SetLowELimit(index, p, 11 * eV);
183 SetHighELimit(index, p, 1 * MeV);
184 }
185 LoadCrossSectionData(p);
186 G4DNAMaterialManager::Instance()->SetMasterDataModel(DNAModelType::fDNAIonisation, this);
187 fpModelData = this;
188 }
189 else {
190 auto dataModel = dynamic_cast<G4DNACPA100IonisationModel*>(
191 G4DNAMaterialManager::Instance()->GetModel(DNAModelType::fDNAIonisation));
192 if (dataModel == nullptr) {
193 G4cout << "G4DNACPA100IonisationModel::CrossSectionPerVolume:: not good modelData" << G4endl;
194 throw;
195 }
196 fpModelData = dataModel;
197 }
198
199 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
200
201 fParticleChangeForGamma = GetParticleChangeForGamma();
202 isInitialised = true;
203 }
204
205 G4double G4DNACPA100IonisationModel::CrossSectionPerVolume(const G4Material* material,
206 const G4ParticleDefinition* p,
207 G4double ekin, G4double, G4double)
208 {
209 // initialise the cross section value (output value)
210 G4double sigma(0);
211
212 // Get the current particle name
213 const G4String& particleName = p->GetParticleName();
214
215 if (p != fpParticle) {
216 G4Exception("G4DNACPA100IonisationModel::CrossSectionPerVolume", "em00223", FatalException,
217 "No model is registered for this particle");
218 }
219
220 auto matID = material->GetIndex();
221
222 // Set the low and high energy limits
223 G4double lowLim = fpModelData->GetLowELimit(matID, p);
224 G4double highLim = fpModelData->GetHighELimit(matID, p);
225
226 // Check that we are in the correct energy range
227 if (ekin >= lowLim && ekin < highLim) {
228 // Get the map with all the model data tables
229 auto tableData = fpModelData->GetData();
230
231 if ((*tableData)[matID][p] == nullptr) {
232 G4Exception("G4DNACPA100IonisationModel::CrossSectionPerVolume", "em00236", FatalException,
233 "No model is registered");
234 }
235 else {
236 sigma = (*tableData)[matID][p]->FindValue(ekin);
237 }
238
239 if (verboseLevel > 2) {
240 auto MolDensity =
241 (*G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(material))[matID];
242 G4cout << "__________________________________" << G4endl;
243 G4cout << "°°° G4DNACPA100IonisationModel - XS INFO START" << G4endl;
244 G4cout << "°°° Kinetic energy(eV)=" << ekin / eV << " particle : " << particleName << G4endl;
245 G4cout << "°°° lowLim (eV) = " << lowLim / eV << " highLim (eV) : " << highLim / eV << G4endl;
246 G4cout << "°°° Materials = " << (*G4Material::GetMaterialTable())[matID]->GetName() << G4endl;
247 G4cout << "°°° Cross section per " << matID << " index molecule (cm^2)=" << sigma / cm / cm
248 << G4endl;
249 G4cout << "°°° Cross section per Phosphate molecule (cm^-1)="
250 << sigma * MolDensity / (1. / cm) << G4endl;
251 G4cout << "°°° G4DNACPA100IonisationModel - XS INFO END" << G4endl;
252 }
253 }
254
255 auto MolDensity = (*G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(material))[matID];
256 return sigma * MolDensity;
257 }
258
259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260
261 void G4DNACPA100IonisationModel::SampleSecondaries(
262 std::vector<G4DynamicParticle*>* fvect,
263 const G4MaterialCutsCouple* couple, // must be set!
264 const G4DynamicParticle* particle, G4double, G4double)
265 {
266 if (verboseLevel > 3) {
267 G4cout << "Calling SampleSecondaries() of G4DNACPA100IonisationModel" << G4endl;
268 }
269 auto k = particle->GetKineticEnergy();
270
271 const G4Material* material = couple->GetMaterial();
272
273 auto MatID = material->GetIndex();
274
275 auto p = particle->GetDefinition();
276
277 auto lowLim = fpModelData->GetLowELimit(MatID, p);
278 auto highLim = fpModelData->GetHighELimit(MatID, p);
279
280 // Check if we are in the correct energy range
281 if (k >= lowLim && k < highLim) {
282 const auto& primaryDirection = particle->GetMomentumDirection();
283 auto particleMass = particle->GetDefinition()->GetPDGMass();
284 auto totalEnergy = k + particleMass;
285 auto pSquare = k * (totalEnergy + particleMass);
286 auto totalMomentum = std::sqrt(pSquare);
287 G4int shell = -1;
288 G4double bindingEnergy, secondaryKinetic;
289 shell = fpModelData->RandomSelectShell(k, p, MatID);
290 bindingEnergy = iStructure.IonisationEnergy(shell, MatID);
291
292 if (k < bindingEnergy) {
293 return;
294 }
295
296 auto info = std::make_tuple(MatID, k, shell);
297
298 secondaryKinetic = -1000 * eV;
299 if (fpG4_WATER->GetIndex() != MatID) {//for DNA material useDcs = false
300 secondaryKinetic = fpModelData->RandomizeEjectedElectronEnergyFromanalytical(info);
301 }else if(fasterCode){
302 secondaryKinetic = fpModelData->RandomizeEjectedElectronEnergyFromCumulatedDcs(info);
303 }else{
304 secondaryKinetic = fpModelData->RandomizeEjectedElectronEnergy(info);
305 }
306
307 G4double cosTheta = 0.;
308 G4double phi = 0.;
309 RandomizeEjectedElectronDirection(particle->GetDefinition(), k, secondaryKinetic, cosTheta,
310 phi);
311
312 G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
313 G4double dirX = sinTheta * std::cos(phi);
314 G4double dirY = sinTheta * std::sin(phi);
315 G4double dirZ = cosTheta;
316 G4ThreeVector deltaDirection(dirX, dirY, dirZ);
317 deltaDirection.rotateUz(primaryDirection);
318
319 // SI - For atom. deexc. tagging - 23/05/2017
320 if (secondaryKinetic > 0) {
321 auto dp = new G4DynamicParticle(G4Electron::Electron(), deltaDirection, secondaryKinetic);
322 fvect->push_back(dp);
323 }
324
325 if (particle->GetDefinition() != fpParticle) {
326 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
327 }
328 else {
329 G4double deltaTotalMomentum =
330 std::sqrt(secondaryKinetic * (secondaryKinetic + 2. * electron_mass_c2));
331 G4double finalPx =
332 totalMomentum * primaryDirection.x() - deltaTotalMomentum * deltaDirection.x();
333 G4double finalPy =
334 totalMomentum * primaryDirection.y() - deltaTotalMomentum * deltaDirection.y();
335 G4double finalPz =
336 totalMomentum * primaryDirection.z() - deltaTotalMomentum * deltaDirection.z();
337 G4double finalMomentum = std::sqrt(finalPx * finalPx + finalPy * finalPy + finalPz * finalPz);
338 finalPx /= finalMomentum;
339 finalPy /= finalMomentum;
340 finalPz /= finalMomentum;
341
342 G4ThreeVector direction;
343 direction.set(finalPx, finalPy, finalPz);
344
345 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
346 }
347
348 // SI - For atom. deexc. tagging - 23/05/2017
349
350 // AM: sample deexcitation
351 // here we assume that H_{2}O electronic levels are the same of Oxigen.
352 // this can be considered true with a rough 10% error in energy on K-shell,
353
354 G4double scatteredEnergy = k - bindingEnergy - secondaryKinetic;
355
356 // SI: only atomic deexcitation from K shell is considered
357 // Hoang: only for water
358 if (material == G4Material::GetMaterial("G4_WATER")) {
359 std::size_t secNumberInit = 0; // need to know at a certain point the energy of secondaries
360 std::size_t secNumberFinal = 0; // So I'll make the diference and then sum the energies
361 if ((fAtomDeexcitation != nullptr) && shell == 4) {
362 G4int Z = 8;
363 auto Kshell = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
364 secNumberInit = fvect->size();
365 fAtomDeexcitation->GenerateParticles(fvect, Kshell, Z, 0, 0);
366 secNumberFinal = fvect->size();
367 if (secNumberFinal > secNumberInit) {
368 for (std::size_t i = secNumberInit; i < secNumberFinal; ++i) {
369 // Check if there is enough residual energy
370 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy()) {
371 // Ok, this is a valid secondary: keep it
372 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
373 }
374 else {
375 // Invalid secondary: not enough energy to create it!
376 // Keep its energy in the local deposit
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::SampleSecondaries()", "em2050", FatalException,
388 "Negative local energy deposit");
389 }
390 if (!statCode) {
391 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
392 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
393 }
394 else {
395 fParticleChangeForGamma->SetProposedKineticEnergy(k);
396 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k - scatteredEnergy);
397 }
398
399 // only water for chemistry
400 if (fpG4_WATER != nullptr && material == G4Material::GetMaterial("G4_WATER")) {
401 const G4Track* theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
402 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule, shell,
403 theIncomingTrack);
404 }
405 }
406 }
407
408 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
409
410 G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergy(PartKineticInMat info)
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(shell, MatID);
417 (k + IonLevel) / 2. > k ? maximumEnergyTransfer = k : maximumEnergyTransfer = (k + IonLevel) / 2.;
418
419 G4double crossSectionMaximum = 0.;
420
421 G4double minEnergy = IonLevel;
422 G4double maxEnergy = maximumEnergyTransfer;
423
424 // nEnergySteps can be optimized - 100 by default
425 G4int nEnergySteps = 50;
426
427 G4double value(minEnergy);
428 G4double stpEnergy(std::pow(maxEnergy / value, 1. / static_cast<G4double>(nEnergySteps - 1)));
429 G4int step(nEnergySteps);
430 G4double differentialCrossSection = 0.;
431 while (step > 0) {
432 step--;
433 differentialCrossSection = DifferentialCrossSection(info, value / eV);
434
435 if (differentialCrossSection > 0) {
436 crossSectionMaximum = differentialCrossSection;
437 break;
438 }
439 value *= stpEnergy;
440 }
441
442 G4double secondaryElectronKineticEnergy = 0.;
443 do {
444 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer - IonLevel);
445 } while (G4UniformRand() * crossSectionMaximum
446 > DifferentialCrossSection(info, (secondaryElectronKineticEnergy + IonLevel) / eV));
447
448 return secondaryElectronKineticEnergy;
449 }
450
451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
452
453 void G4DNACPA100IonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition*,
454 G4double k, G4double secKinetic,
455 G4double& cosTheta,
456 G4double& phi)
457 {
458 phi = twopi * G4UniformRand();
459 G4double sin2O = (1. - secKinetic / k) / (1. + secKinetic / (2. * electron_mass_c2));
460 cosTheta = std::sqrt(1. - sin2O);
461 }
462
463 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
464
465 G4double G4DNACPA100IonisationModel::DifferentialCrossSection(PartKineticInMat info,
466 const G4double& energyTransfer)
467 {
468 auto MatID = std::get<0>(info);
469 auto k = std::get<1>(info) / eV; // in eV unit
470 auto shell = std::get<2>(info);
471 G4double sigma = 0.;
472 G4double shellEnergy = iStructure.IonisationEnergy(shell, MatID);
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[MatID][fpParticle].begin(),
489 fTMapWithVec[MatID][fpParticle].end(), k);
490 auto t1 = t2 - 1;
491
492 if (kSE <= fEMapWithVector[MatID][fpParticle][(*t1)].back()
493 && kSE <= fEMapWithVector[MatID][fpParticle][(*t2)].back())
494 {
495 auto e12 = std::upper_bound(fEMapWithVector[MatID][fpParticle][(*t1)].begin(),
496 fEMapWithVector[MatID][fpParticle][(*t1)].end(), kSE);
497 auto e11 = e12 - 1;
498
499 auto e22 = std::upper_bound(fEMapWithVector[MatID][fpParticle][(*t2)].begin(),
500 fEMapWithVector[MatID][fpParticle][(*t2)].end(), kSE);
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][fpParticle][shell][valueT1][valueE11];
511 xs12 = diffCrossSectionData[MatID][fpParticle][shell][valueT1][valueE12];
512 xs21 = diffCrossSectionData[MatID][fpParticle][shell][valueT2][valueE21];
513 xs22 = diffCrossSectionData[MatID][fpParticle][shell][valueT2][valueE22];
514 }
515
516 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
517
518 if (xsProduct != 0.) {
519 sigma = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12, xs21, xs22,
520 valueT1, valueT2, k, kSE);
521 }
522 }
523
524 return sigma;
525 }
526
527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
528
529 G4double G4DNACPA100IonisationModel::Interpolate(G4double e1, G4double e2, G4double e, G4double xs1,
530 G4double xs2)
531 {
532 G4double value = 0.;
533
534 // Log-log interpolation by default
535
536 if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0 && !fasterCode) {
537 G4double a = (std::log10(xs2) - std::log10(xs1)) / (std::log10(e2) - std::log10(e1));
538 G4double b = std::log10(xs2) - a * std::log10(e2);
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 faster code
554
555 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode) {
556 G4double d1 = std::log10(xs1);
557 G4double d2 = std::log10(xs2);
558 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
559 }
560
561 // Switch to lin-lin interpolation for faster code
562 // in case one of xs1 or xs2 (=cum proba) value is zero
563
564 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode) {
565 G4double d1 = xs1;
566 G4double d2 = xs2;
567 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
568 }
569 return value;
570 }
571
572 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
573
574 G4double G4DNACPA100IonisationModel::QuadInterpolator(G4double e11, G4double e12, G4double e21,
575 G4double e22, G4double xs11, G4double xs12,
576 G4double xs21, G4double xs22, G4double t1,
577 G4double t2, G4double t, G4double e)
578 {
579 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
580 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
581 G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
582
583 return value;
584 }
585
586 G4double
587 G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(PartKineticInMat info)
588 {
589 auto MatID = std::get<0>(info);
590 auto shell = std::get<2>(info);
591 G4double secondaryElectronKineticEnergy =
592 RandomTransferedEnergy(info) * eV - iStructure.IonisationEnergy(shell, MatID);
593 if (secondaryElectronKineticEnergy < 0.) {
594 return 0.;
595 }
596 return secondaryElectronKineticEnergy;
597 }
598
599 G4double G4DNACPA100IonisationModel::RandomTransferedEnergy(PartKineticInMat info)
600 {
601 auto materialID = std::get<0>(info);
602 auto k = std::get<1>(info) / eV; // data table in eV
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].back()) {
617 k = k * (1. - 1e-12);
618 }
619
620 G4double random = G4UniformRand();
621 auto k2 = std::upper_bound(fTMapWithVec[materialID][fpParticle].begin(),
622 fTMapWithVec[materialID][fpParticle].end(), k);
623 auto k1 = k2 - 1;
624
625 if (random <= fProbaShellMap[materialID][fpParticle][shell][(*k1)].back()
626 && random <= fProbaShellMap[materialID][fpParticle][shell][(*k2)].back())
627 {
628 auto cumulCS12 =
629 std::upper_bound(fProbaShellMap[materialID][fpParticle][shell][(*k1)].begin(),
630 fProbaShellMap[materialID][fpParticle][shell][(*k1)].end(), random);
631 auto cumulCS11 = cumulCS12 - 1;
632 // Second one.
633 auto cumulCS22 =
634 std::upper_bound(fProbaShellMap[materialID][fpParticle][shell][(*k2)].begin(),
635 fProbaShellMap[materialID][fpParticle][shell][(*k2)].end(), random);
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[materialID][fpParticle][shell][valueK1][valueCumulCS11];
646 secElecE12 = fEnergySecondaryData[materialID][fpParticle][shell][valueK1][valueCumulCS12];
647 secElecE21 = fEnergySecondaryData[materialID][fpParticle][shell][valueK2][valueCumulCS21];
648 secElecE22 = fEnergySecondaryData[materialID][fpParticle][shell][valueK2][valueCumulCS22];
649
650 if (valueCumulCS11 == 0. && valueCumulCS12 == 1.) {
651 auto interpolatedvalue2 =
652 Interpolate(valueCumulCS21, valueCumulCS22, random, secElecE21, secElecE22);
653 G4double valueNrjTransf = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
654 return valueNrjTransf;
655 }
656 }
657
658 if (random > fProbaShellMap[materialID][fpParticle][shell][(*k1)].back()) {
659 auto cumulCS22 =
660 std::upper_bound(fProbaShellMap[materialID][fpParticle][shell][(*k2)].begin(),
661 fProbaShellMap[materialID][fpParticle][shell][(*k2)].end(), random);
662 auto cumulCS21 = cumulCS22 - 1;
663 valueK1 = *k1;
664 valueK2 = *k2;
665 valueCumulCS21 = *cumulCS21;
666 valueCumulCS22 = *cumulCS22;
667
668 secElecE21 = fEnergySecondaryData[materialID][fpParticle][shell][valueK2][valueCumulCS21];
669 secElecE22 = fEnergySecondaryData[materialID][fpParticle][shell][valueK2][valueCumulCS22];
670
671 G4double interpolatedvalue2 =
672 Interpolate(valueCumulCS21, valueCumulCS22, random, secElecE21, secElecE22);
673
674 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
675 return value;
676 }
677 G4double nrjTransfProduct = secElecE11 * secElecE12 * secElecE21 * secElecE22;
678
679 if (nrjTransfProduct != 0.) {
680 ejectedElectronEnergy =
681 QuadInterpolator(valueCumulCS11, valueCumulCS12, valueCumulCS21, valueCumulCS22, secElecE11,
682 secElecE12, secElecE21, secElecE22, valueK1, valueK2, k, random);
683 }
684 return ejectedElectronEnergy;
685 }
686
687 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
688
689 G4double
690 G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergyFromanalytical(PartKineticInMat info)
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 ***** (optimized)
696 // Composition sampling method based on eq 7 in (Guerra et al. 2015) the RBEBV
697
698 //// Defining constants
699 G4double alfa = 1. / 137; // fine structure constant
700 G4double e_charge = 1.6e-19; // electron charge
701 G4double e_mass = 9.1e-31; // electron mass in kg
702 G4double c = 3e8; // speed of light in vacuum constant c (m/s)
703 G4double mc2 = e_mass * c * c / e_charge; //
704
705 G4double BB = iStructure.IonisationEnergy(shell, MatID); // binding energy of the shell (eV)
706
707 if (tt <= BB) return 0.;
708
709 G4double b_prime = BB / mc2; // binding energy divided by mc2
710 G4double beta_b2 = 1. - 1. / ((1 + b_prime) * (1 + b_prime)); // binding energy Beta
711
712 //// Indicent energy
713 //// tt is the incident electron energy
714
715 G4double t_prime = tt / mc2; // incident energy divided by mc2
716 G4double t = tt / BB; // reduced incident energy by binding energy
717
718 G4double D = (1 + 2 * t_prime) / ((1 + t_prime / 2) * (1 + t_prime / 2));
719 G4double F = b_prime * b_prime / ((1 + t_prime / 2) * (1 + t_prime / 2));
720
721 G4double beta_t2 = 1 - 1 / ((1 + t_prime) * (1 + t_prime)); // incident energy Beta
722
723 G4double PHI_R = std::cos(std::sqrt(alfa * alfa / (beta_t2 + beta_b2))
724 * std::log(beta_t2 / beta_b2)); // relativistic Vriens function phi
725 G4double G_R = std::log(beta_t2 / (1 - beta_t2)) - beta_t2 - std::log(2 * b_prime);
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 * t + 1) / (2 * t * tplus1));
731 G4double ZH2max = 1 - PHI_R * D / 4;
732
733 G4double A1_p = ZH1max * tminus1 / tplus1; // A1'
734 G4double A2_p = ZH2max * tminus1 / (t * tplus1); // A2'
735 G4double A3_p = ((tplus12 - 4) / tplus12) * G_R; // A3'
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 * tplus1));
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........oooOO0OOooo........oooOO0OOooo....
801
802 void G4DNACPA100IonisationModel::ReadDiffCSFile(const std::size_t& materialID,
803 const G4ParticleDefinition* p, const G4String& file,
804 const G4double& scaleFactor)
805 {
806 const char* path = G4FindDataDir("G4LEDATA");
807 if (path == nullptr) {
808 G4Exception("G4DNACPA100IonisationModel::ReadAllDiffCSFiles", "em0006", FatalException,
809 "G4LEDATA environment variable was not set.");
810 return;
811 }
812
813 std::ostringstream fullFileName;
814 fullFileName << path << "/" << file << ".dat";
815
816 std::ifstream diffCrossSection(fullFileName.str().c_str());
817 std::stringstream endPath;
818 if (!diffCrossSection) {
819 endPath << "Missing data file: " << file;
820 G4Exception("G4DNACPA100IonisationModel::Initialise", "em0003", FatalException,
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 transferred
838 if (T != fTMapWithVec[materialID][p].back()) {
839 fTMapWithVec[materialID][p].push_back(T);
840 }
841
842 auto eshell = (G4int)iStructure.NumberOfLevels(materialID);
843 for (G4int shell = 0; shell < eshell; ++shell) {
844 diffCrossSection >> diffCrossSectionData[materialID][p][shell][T][E];
845 if (fasterCode) {
846 fEnergySecondaryData[materialID][p][shell][T]
847 [diffCrossSectionData[materialID][p][shell][T][E]] = E;
848
849 fProbaShellMap[materialID][p][shell][T].push_back(
850 diffCrossSectionData[materialID][p][shell][T][E]);
851 }
852 else {
853 diffCrossSectionData[materialID][p][shell][T][E] *= scaleFactor;
854 fEMapWithVector[materialID][p][T].push_back(E);
855 }
856 }
857 }
858 }
859