| 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 // Authors: S. Meylan and C. Villagrasa (IRSN, France)
27 // Models come from
28 // M. Bug et al, Rad. Phys and Chem. 130, 459-479 (2017)
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::G4DNAPTBIonisationModel(const G4String& applyToMaterial,
40 const G4ParticleDefinition*, const G4String& nam,
41 const G4bool isAuger)
42 : G4VDNAModel(nam, applyToMaterial)
43 {
44 if (isAuger) {
45 // create the PTB Auger model
46 fpDNAPTBAugerModel = std::make_unique<G4DNAPTBAugerModel>("e-_G4DNAPTBAugerModel");
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_WATER", false);
53 fpBackbone_THF = G4Material::GetMaterial("backbone_THF", false);
54 fpCytosine_PY = G4Material::GetMaterial("cytosine_PY", false);
55 fpThymine_PY = G4Material::GetMaterial("thymine_PY", false);
56 fpAdenine_PU = G4Material::GetMaterial("adenine_PU", false);
57 fpBackbone_TMP = G4Material::GetMaterial("backbone_TMP", false);
58 fpGuanine_PU = G4Material::GetMaterial("guanine_PU", false);
59 fpN2 = G4Material::GetMaterial("N2", false);
60 }
61
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63
64 void G4DNAPTBIonisationModel::Initialise(const G4ParticleDefinition* particle,
65 const G4DataVector& /*cuts*/)
66 {
67 if (isInitialised) {
68 return;
69 }
70 if (verboseLevel > 3) {
71 G4cout << "Calling G4DNAPTBIonisationModel::Initialise()" << G4endl;
72 }
73
74 G4double scaleFactor = 1e-16 * cm * cm;
75 G4double scaleFactorBorn = (1.e-22 / 3.343) * m * m;
76
77 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
78 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
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, "dna/sigma_ionisation_e-_PTB_N2",
91 "dna/sigmadiff_cumulated_ionisation_e-_PTB_N2", scaleFactor);
92 SetLowELimit(index, particle, 15.5 * eV);
93 SetHighELimit(index, particle, 1.02 * MeV);
94 }
95
96 // MPietrzak
97
98 if (fpTHF != nullptr) {
99 index = fpTHF->GetIndex();
100 AddCrossSectionData(index, particle, "dna/sigma_ionisation_e-_PTB_THF",
101 "dna/sigmadiff_cumulated_ionisation_e-_PTB_THF", scaleFactor);
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, "dna/sigma_ionisation_e-_PTB_PY",
109 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PY", scaleFactor);
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, "dna/sigma_ionisation_e-_PTB_PU",
117 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU", scaleFactor);
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, "dna/sigma_ionisation_e-_PTB_TMP",
124 "dna/sigmadiff_cumulated_ionisation_e-_PTB_TMP", scaleFactor);
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, "dna/sigma_ionisation_e_born",
132 "dna/sigmadiff_ionisation_e_born", scaleFactorBorn);
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, "dna/sigma_ionisation_e-_PTB_THF",
141 "dna/sigmadiff_cumulated_ionisation_e-_PTB_THF", scaleFactor * 33. / 30);
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, "dna/sigma_ionisation_e-_PTB_PY",
149 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PY", scaleFactor * 42. / 30);
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, "dna/sigma_ionisation_e-_PTB_PY",
157 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PY", scaleFactor * 48. / 30);
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, "dna/sigma_ionisation_e-_PTB_PU",
165 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU", scaleFactor * 50. / 44);
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, "dna/sigma_ionisation_e-_PTB_PU",
173 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU", scaleFactor * 56. / 44);
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, "dna/sigma_ionisation_e-_PTB_TMP",
181 "dna/sigmadiff_cumulated_ionisation_e-_PTB_TMP", scaleFactor * 33. / 50);
182 SetLowELimit(index, particle, 12. * eV);
183 SetHighELimit(index, particle, 1. * keV);
184 }
185 }
186
187 else if (particle == protonDef) {
188 G4String particleName = particle->GetParticleName();
189
190 // Raw materials
191 //
192 if (fpTHF != nullptr) {
193 index = fpTHF->GetIndex();
194 AddCrossSectionData(index, particle, "dna/sigma_ionisation_p_HKS_THF",
195 "dna/sigmadiff_cumulated_ionisation_p_PTB_THF", scaleFactor);
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, "dna/sigma_ionisation_p_HKS_PY",
203 "dna/sigmadiff_cumulated_ionisation_p_PTB_PY", scaleFactor);
204 SetLowELimit(index, particle, 70. * keV);
205 SetHighELimit(index, particle, 10. * MeV);
206 }
207 /*
208 AddCrossSectionData("PU",
209 particleName,
210 "dna/sigma_ionisation_e-_PTB_PU",
211 "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU",
212 scaleFactor);
213 SetLowELimit("PU", particleName2, 70.*keV);
214 SetHighELimit("PU", particleName2, 10.*keV);
215 */
216
217 if (fpTMP != nullptr) {
218 index = fpTMP->GetIndex();
219 AddCrossSectionData(index, particle, "dna/sigma_ionisation_p_HKS_TMP",
220 "dna/sigmadiff_cumulated_ionisation_p_PTB_TMP", scaleFactor);
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()->IsLocked()) {
229 LoadCrossSectionData(particle);
230 G4DNAMaterialManager::Instance()->SetMasterDataModel(DNAModelType::fDNAIonisation, this);
231 fpModelData = this;
232 }
233 else {
234 auto dataModel = dynamic_cast<G4DNAPTBIonisationModel*>(
235 G4DNAMaterialManager::Instance()->GetModel(DNAModelType::fDNAIonisation));
236 if (dataModel == nullptr) {
237 G4cout << "G4DNAPTBIonisationModel::Initialise:: not good modelData" << G4endl;
238 G4Exception("G4DNAPTBIonisationModel::Initialise", "PTB0004", FatalException,
239 "not good modelData");
240 }
241 else {
242 fpModelData = dataModel;
243 }
244 }
245 // initialise DNAPTBAugerModel
246 if (fpDNAPTBAugerModel) {
247 fpDNAPTBAugerModel->Initialise();
248 }
249 fParticleChangeForGamma = GetParticleChangeForGamma();
250 isInitialised = true;
251 }
252
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254
255 G4double G4DNAPTBIonisationModel::CrossSectionPerVolume(const G4Material* material,
256 const G4ParticleDefinition* p,
257 G4double ekin, G4double /*emin*/,
258 G4double /*emax*/)
259 {
260 // initialise the cross section value (output value)
261 G4double sigma(0);
262
263 // Get the current particle name
264 const G4String& particleName = p->GetParticleName();
265 const std::size_t& matID = material->GetIndex();
266
267 // Set the low and high energy limits
268 G4double lowLim = fpModelData->GetLowELimit(matID, p);
269 G4double highLim = fpModelData->GetHighELimit(matID, p);
270
271 // Check that we are in the correct energy range
272 if (ekin >= lowLim && ekin < highLim) {
273 // Get the map with all the model data tables
274 auto tableData = fpModelData->GetData();
275 if ((*tableData)[matID][p] == nullptr) {
276 G4Exception("G4DNAPTBIonisationModel::CrossSectionPerVolume", "em00236", FatalException,
277 "No model is registered");
278 }
279 // Retrieve the cross section value for the current material, particle and energy values
280 sigma = (*tableData)[matID][p]->FindValue(ekin);
281
282 if (verboseLevel > 2) {
283 G4cout << "__________________________________" << G4endl;
284 G4cout << "°°° G4DNAPTBIonisationModel - XS INFO START" << G4endl;
285 G4cout << "°°° Kinetic energy(eV)=" << ekin / eV << " particle : " << particleName << G4endl;
286 G4cout << "°°° Cross section per " << matID << " index molecule (cm^2)=" << sigma / cm / cm
287 << G4endl;
288 G4cout << "°°° G4DNAPTBIonisationModel - XS INFO END" << G4endl;
289 }
290 }
291
292 // Return the cross section value
293 auto MolDensity =
294 (*G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(material))[matID];
295 return sigma * MolDensity;
296 }
297
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299
300 void G4DNAPTBIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
301 const G4MaterialCutsCouple* pCouple,
302 const G4DynamicParticle* aDynamicParticle,
303 G4double /*tmin*/, G4double /*tmax*/)
304 {
305 // Get the current particle energy
306 G4double k = aDynamicParticle->GetKineticEnergy();
307 const std::size_t& materialID = pCouple->GetMaterial()->GetIndex();
308
309 // Get the current particle name
310 const auto& p = aDynamicParticle->GetDefinition();
311 const auto& materialName = pCouple->GetMaterial()->GetName();
312 // Get the energy limits
313 G4double lowLim = fpModelData->GetLowELimit(materialID, p);
314 G4double highLim = fpModelData->GetHighELimit(materialID, p);
315
316 // Check if we are in the correct energy range
317 if (k >= lowLim && k < highLim) {
318 G4ParticleMomentum primaryDirection = aDynamicParticle->GetMomentumDirection();
319 G4double particleMass = aDynamicParticle->GetDefinition()->GetPDGMass();
320 G4double totalEnergy = k + particleMass;
321 G4double pSquare = k * (totalEnergy + particleMass);
322 G4double totalMomentum = std::sqrt(pSquare);
323
324 // Get the ionisation shell from a random sampling
325 G4int ionizationShell = fpModelData->RandomSelectShell(k, p, materialID);
326
327 // Get the binding energy from the ptbStructure class
328 G4double bindingEnergy = ptbStructure.IonisationEnergy(ionizationShell, materialID);
329
330 // Initialize the secondary kinetic energy to a negative value.
331 G4double secondaryKinetic(-1000 * eV);
332
333 if (fpG4_WATER == nullptr || materialID != fpG4_WATER->GetIndex()) {
334 // Get the energy of the secondary particle
335 secondaryKinetic = fpModelData->RandomizeEjectedElectronEnergyFromCumulated(
336 aDynamicParticle->GetDefinition(), k / eV, ionizationShell, materialID);
337 }
338 else {
339 secondaryKinetic = fpModelData->RandomizeEjectedElectronEnergy(
340 aDynamicParticle->GetDefinition(), k, ionizationShell, materialID);
341 }
342
343 if (secondaryKinetic <= 0) {
344 G4cout << "Fatal error *************************************** " << secondaryKinetic / eV
345 << G4endl;
346 G4cout << "secondaryKinetic: " << secondaryKinetic / eV << G4endl;
347 G4cout << "k: " << k / eV << G4endl;
348 G4cout << "shell: " << ionizationShell << G4endl;
349 G4cout << "material:" << materialName << G4endl;
350 G4Exception("G4DNAPTBIonisationModel::SampleSecondaries", "em0026", FatalException,
351 "Fatal error:: scatteredEnergy <= 0");
352 }
353
354 G4double cosTheta = 0.;
355 G4double phi = 0.;
356 RandomizeEjectedElectronDirection(aDynamicParticle->GetDefinition(), k, secondaryKinetic,
357 cosTheta, phi);
358
359 G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
360 G4double dirX = sinTheta * std::cos(phi);
361 G4double dirY = sinTheta * std::sin(phi);
362 G4double dirZ = cosTheta;
363 G4ThreeVector deltaDirection(dirX, dirY, dirZ);
364 deltaDirection.rotateUz(primaryDirection);
365
366 // The model is written only for electron and thus we want the change the direction of the
367 // incident electron after each ionization. However, if other particle are going to be
368 // introduced within this model the following should be added:
369 //
370 // Check if the particle is an electron
371 if (aDynamicParticle->GetDefinition() == G4Electron::ElectronDefinition()) {
372 // If yes do the following code until next commented "else" statement
373
374 G4double deltaTotalMomentum =
375 std::sqrt(secondaryKinetic * (secondaryKinetic + 2. * electron_mass_c2));
376 G4double finalPx =
377 totalMomentum * primaryDirection.x() - deltaTotalMomentum * deltaDirection.x();
378 G4double finalPy =
379 totalMomentum * primaryDirection.y() - deltaTotalMomentum * deltaDirection.y();
380 G4double finalPz =
381 totalMomentum * primaryDirection.z() - deltaTotalMomentum * deltaDirection.z();
382 G4double finalMomentum = std::sqrt(finalPx * finalPx + finalPy * finalPy + finalPz * finalPz);
383 finalPx /= finalMomentum;
384 finalPy /= finalMomentum;
385 finalPz /= finalMomentum;
386
387 G4ThreeVector direction(finalPx, finalPy, finalPz);
388 if (direction.unit().getX() > 1 || direction.unit().getY() > 1 || direction.unit().getZ() > 1)
389 {
390 G4cout << "Fatal error ****************************" << G4endl;
391 G4cout << "direction problem " << direction.unit() << G4endl;
392 G4Exception("G4DNAPTBIonisationModel::SampleSecondaries", "em0017", FatalException,
393 "Fatal error:: direction problem");
394 }
395
396 // Give the new direction to the particle
397 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
398 }
399 // If the particle is not an electron
400 else {
401 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
402 }
403
404 // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
405 G4double scatteredEnergy = k - bindingEnergy - secondaryKinetic;
406
407 if (scatteredEnergy <= 0) {
408 G4cout << "Fatal error ****************************" << G4endl;
409 G4cout << "k: " << k / eV << G4endl;
410 G4cout << "secondaryKinetic: " << secondaryKinetic / eV << G4endl;
411 G4cout << "shell: " << ionizationShell << G4endl;
412 G4cout << "bindingEnergy: " << bindingEnergy / eV << G4endl;
413 G4cout << "scatteredEnergy: " << scatteredEnergy / eV << G4endl;
414 G4cout << "material: " << materialName << G4endl;
415 G4Exception("G4DNAPTBIonisationModel::SampleSecondaries", "em0016", FatalException,
416 "Fatal error:: scatteredEnergy <= 0");
417 }
418
419 // Set the new energy of the particle
420 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
421
422 // Set the energy deposited by the ionization
423 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k - scatteredEnergy - secondaryKinetic);
424
425 // Create the new particle with its characteristics
426 auto dp = new G4DynamicParticle(G4Electron::Electron(), deltaDirection, secondaryKinetic);
427 fvect->push_back(dp);
428
429 // Check if the auger model is activated (ie instanciated)
430 if (fpDNAPTBAugerModel) {
431 // run the PTB Auger model
432 if (materialName != "G4_WATER") {
433 fpDNAPTBAugerModel->ComputeAugerEffect(fvect, materialName, bindingEnergy);
434 }
435 }
436 }
437 }
438
439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
440
441 void G4DNAPTBIonisationModel::ReadDiffCSFile(const std::size_t& materialID,
442 const G4ParticleDefinition* p, const G4String& file,
443 const G4double& scaleFactor)
444 {
445 // To read and save the informations contained within the differential cross section files
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 error message
450 if (path == nullptr) {
451 G4Exception("G4DNAPTBIonisationModel::ReadAllDiffCSFiles", "em0006", FatalException,
452 "G4LEDATA environment variable was not set.");
453 return;
454 }
455
456 // build the fullFileName path of the data file
457 std::ostringstream fullFileName;
458 fullFileName << path << "/" << file << ".dat";
459
460 // open the data file
461 std::ifstream diffCrossSection(fullFileName.str().c_str());
462 // error if file is not there
463 std::stringstream endPath;
464 if (!diffCrossSection) {
465 endPath << "Missing data file: " << file;
466 G4Exception("G4DNAPTBIonisationModel::Initialise", "em0003", FatalException,
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 file point
476 // fill fTMapWithVec, diffCrossSectionData, fEnergyTransferData, fProbaShellMap and
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 following information is data or comments
485 if (test == "#") {
486 // skip the line by beginning a new while loop.
487 continue;
488 }
489 // check if line is empty
490 if (line.empty()) {
491 // skip the line by beginning a new while loop.
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 numbers of each file line
505 iss >> T >> E;
506
507 // Fill the fTMapWithVec container with all the different T values contained within the file.
508 // Duplicate must be avoided and this is the purpose of the if statement
509 if (T != fTMapWithVec[materialID][p].back()) fTMapWithVec[materialID][p].push_back(T);
510
511 // iterate on each shell of the corresponding material
512 for (int shell = 0, eshell = ptbStructure.NumberOfLevels(materialID); shell < eshell; ++shell) {
513 // map[material][particle][shell][T][E]=diffCrossSectionValue
514 // Fill the map with the informations of the input file
515 iss >> diffCrossSectionData[materialID][p][shell][T][E];
516
517 if (fpG4_WATER != nullptr && fpG4_WATER->GetIndex() != materialID) {
518 // map[material][particle][shell][T][CS]=E
519 // Fill the map
520 fEnergySecondaryData[materialID][p][shell][T]
521 [diffCrossSectionData[materialID][p][shell][T][E]] = E;
522
523 // map[material][particle][shell][T]=CS_vector
524 // Fill the vector within the map
525 fProbaShellMap[materialID][p][shell][T].push_back(
526 diffCrossSectionData[materialID][p][shell][T][E]);
527 }
528 else {
529 diffCrossSectionData[materialID][p][shell][T][E] *= scaleFactor;
530 fEMapWithVector[materialID][p][T].push_back(E);
531 }
532 }
533 }
534 }
535
536 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
537
538 G4double G4DNAPTBIonisationModel::RandomizeEjectedElectronEnergy(
539 const G4ParticleDefinition* particleDefinition, G4double k, G4int shell, const std::size_t& materialID)
540 {
541 if (particleDefinition == G4Electron::ElectronDefinition()) {
542 // G4double Tcut=25.0E-6;
543 G4double maximumEnergyTransfer;
544 ((k + ptbStructure.IonisationEnergy(shell, materialID)) / 2. > k)
545 ? maximumEnergyTransfer = k
546 : maximumEnergyTransfer = (k + ptbStructure.IonisationEnergy(shell, materialID)) / 2.;
547
548 // SI : original method
549 /*
550 G4double crossSectionMaximum = 0.;
551 for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer;
552 value+=0.1*eV)
553 {
554 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV,
555 value/eV, shell); if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum =
556 differentialCrossSection;
557 }
558 */
559
560 // SI : alternative method
561
562 // if (k > Tcut)
563 //{
564 G4double crossSectionMaximum = 0.;
565
566 G4double minEnergy = ptbStructure.IonisationEnergy(shell, materialID);
567 G4double maxEnergy = maximumEnergyTransfer;
568 G4int nEnergySteps = 50;
569 G4double value(minEnergy);
570 G4double stpEnergy(std::pow(maxEnergy / value, 1. / static_cast<G4double>(nEnergySteps - 1)));
571 G4int step(nEnergySteps);
572 while (step > 0) {
573 step--;
574 G4double differentialCrossSection =
575 DifferentialCrossSection(particleDefinition, k / eV, value / eV, shell, materialID);
576 if (differentialCrossSection >= crossSectionMaximum)
577 crossSectionMaximum = differentialCrossSection;
578 value *= stpEnergy;
579 }
580 //
581
582 G4double secondaryElectronKineticEnergy = 0.;
583
584 do {
585 secondaryElectronKineticEnergy =
586 G4UniformRand()
587 * (maximumEnergyTransfer - ptbStructure.IonisationEnergy(shell, materialID));
588
589 } while (
590 G4UniformRand() * crossSectionMaximum > DifferentialCrossSection(
591 particleDefinition, k / eV,
592 (secondaryElectronKineticEnergy + ptbStructure.IonisationEnergy(shell, materialID)) / eV,
593 shell, materialID));
594
595 return secondaryElectronKineticEnergy;
596 }
597
598 if (particleDefinition == G4Proton::ProtonDefinition()) {
599 G4double maximumKineticEnergyTransfer = 4. * (electron_mass_c2 / proton_mass_c2) * k;
600
601 G4double crossSectionMaximum = 0.;
602 for (G4double value = ptbStructure.IonisationEnergy(shell, materialID);
603 value <= 4. * ptbStructure.IonisationEnergy(shell, materialID); value += 0.1 * eV)
604 {
605 G4double differentialCrossSection =
606 DifferentialCrossSection(particleDefinition, k / eV, value / eV, shell, materialID);
607 if (differentialCrossSection >= crossSectionMaximum)
608 crossSectionMaximum = differentialCrossSection;
609 }
610
611 G4double secondaryElectronKineticEnergy = 0.;
612 do {
613 secondaryElectronKineticEnergy = G4UniformRand() * maximumKineticEnergyTransfer;
614 } while (
615 G4UniformRand() * crossSectionMaximum >= DifferentialCrossSection(
616 particleDefinition, k / eV,
617 (secondaryElectronKineticEnergy + ptbStructure.IonisationEnergy(shell, materialID)) / eV,
618 shell, materialID));
619
620 return secondaryElectronKineticEnergy;
621 }
622
623 return 0;
624 }
625
626 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
627
628 void G4DNAPTBIonisationModel::RandomizeEjectedElectronDirection(const G4ParticleDefinition* p,
629 G4double k, G4double secKinetic,
630 G4double& cosTheta, G4double& phi)
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::sqrt(2.) / 2);
641 }
642 else {
643 G4double sin2O = (1. - secKinetic / k) / (1. + secKinetic / (2. * electron_mass_c2));
644 cosTheta = std::sqrt(1. - sin2O);
645 }
646 }
647 else if (p == G4Proton::ProtonDefinition()) {
648 G4double maxSecKinetic = 4. * (electron_mass_c2 / proton_mass_c2) * k;
649 phi = twopi * G4UniformRand();
650
651 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
652
653 // Restriction below 100 eV from Emfietzoglou (2000)
654
655 (secKinetic > 100 * eV) ? cosTheta = std::sqrt(secKinetic / maxSecKinetic)
656 : cosTheta = (2. * G4UniformRand()) - 1.;
657 }
658 }
659
660 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
661
662 double G4DNAPTBIonisationModel::DifferentialCrossSection(const G4ParticleDefinition* p, G4double k,
663 G4double energyTransfer,
664 G4int ionizationLevelIndex,
665 const std::size_t& materialID)
666 {
667 G4double sigma = 0.;
668 G4double shellEnergy(ptbStructure.IonisationEnergy(ionizationLevelIndex, materialID));
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 eV also
686 auto t2 =
687 std::upper_bound(fTMapWithVec[materialID][p].begin(), fTMapWithVec[materialID][p].end(), k);
688 auto t1 = t2 - 1;
689
690 // SI : the following condition avoids situations where energyTransfer >last vector element
691 if (kSE <= fEMapWithVector[materialID][p][(*t1)].back()
692 && kSE <= fEMapWithVector[materialID][p][(*t2)].back())
693 {
694 auto e12 = std::upper_bound(fEMapWithVector[materialID][p][(*t1)].begin(),
695 fEMapWithVector[materialID][p][(*t1)].end(), kSE);
696 auto e11 = e12 - 1;
697
698 auto e22 = std::upper_bound(fEMapWithVector[materialID][p][(*t2)].begin(),
699 fEMapWithVector[materialID][p][(*t2)].end(), kSE);
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][p][ionizationLevelIndex][valueT1][valueE11];
710 xs12 = diffCrossSectionData[materialID][p][ionizationLevelIndex][valueT1][valueE12];
711 xs21 = diffCrossSectionData[materialID][p][ionizationLevelIndex][valueT2][valueE21];
712 xs22 = diffCrossSectionData[materialID][p][ionizationLevelIndex][valueT2][valueE22];
713 }
714 }
715
716 if (p == G4Proton::ProtonDefinition()) {
717 // k should be in eV and energy transfer eV also
718 auto t2 =
719 std::upper_bound(fTMapWithVec[materialID][p].begin(), fTMapWithVec[materialID][p].end(), k);
720 auto t1 = t2 - 1;
721
722 auto e12 = std::upper_bound(fEMapWithVector[materialID][p][(*t1)].begin(),
723 fEMapWithVector[materialID][p][(*t1)].end(), kSE);
724 auto e11 = e12 - 1;
725
726 auto e22 = std::upper_bound(fEMapWithVector[materialID][p][(*t2)].begin(),
727 fEMapWithVector[materialID][p][(*t2)].end(), kSE);
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][p][ionizationLevelIndex][valueT1][valueE11];
738 xs12 = diffCrossSectionData[materialID][p][ionizationLevelIndex][valueT1][valueE12];
739 xs21 = diffCrossSectionData[materialID][p][ionizationLevelIndex][valueT2][valueE21];
740 xs22 = diffCrossSectionData[materialID][p][ionizationLevelIndex][valueT2][valueE22];
741 }
742
743 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
744
745 if (xsProduct != 0.) {
746 sigma = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12, xs21, xs22,
747 valueT1, valueT2, k, kSE);
748 }
749 }
750
751 return sigma;
752 }
753
754 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
755
756 G4double G4DNAPTBIonisationModel::RandomizeEjectedElectronEnergyFromCumulated(
757 const G4ParticleDefinition* p, G4double k, G4int ionizationLevelIndex, const std::size_t& materialID)
758 {
759 // k should be in eV
760
761 // Schematic explanation.
762 // We will do an interpolation to get a final E value (ejected electron energy).
763 // 1/ We choose a random number between 0 and 1 (ie we select a cumulated cross section).
764 // 2/ We look for T_lower and T_upper.
765 // 3/ We look for the cumulated corresponding cross sections and their associated E values.
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_low --
776 // | CS_low_2 -> E_low_2 ----- |
777 // | ---> E_final
778 // T_up | CS_up_1 -> E_up_1 ------- |
779 // | |----> E_up ---
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 compare with the cumulated CS
799 // ***************************************************************************
800 //
801 // It will allow us to choose an ejected electron energy with respect to the CS.
802 G4double random = G4UniformRand();
803
804 // **********************************************
805 // Take the input from the data tables
806 // **********************************************
807
808 // Cumulated tables are like this: T E cumulatedCS1 cumulatedCS2 cumulatedCS3
809 // We have two sets of loaded data: fTMapWithVec which contains data about T (incident particle
810 // energy) and fProbaShellMap which contains cumulated cross section data. Since we already have a
811 // specific T energy value which could not be explicitly in the table, we must interpolate all the
812 // values.
813
814 // First, we select the upper and lower T data values surrounding our T value (ie "k").
815 auto k2 =
816 std::upper_bound(fTMapWithVec[materialID][p].begin(), fTMapWithVec[materialID][p].end(), k);
817 auto k1 = k2 - 1;
818
819 // Check if we have found a k2 value (0 if we did not found it).
820 // A missing k2 value can be caused by a energy to high for the data table,
821 // Ex : table done for 12*eV -> 1000*eV and k=2000*eV
822 // then k2 = 0 and k1 = max of the table.
823 // To detect this, we check that k1 is not superior to k2.
824 if (*k1 > *k2) {
825 // Error
826 G4cerr << "**************** Fatal error ******************" << G4endl;
827 G4cerr << "G4DNAPTBIonisationModel::RandomizeEjectedElectronEnergyFromCumulated" << G4endl;
828 G4cerr << "You have *k1 > *k2 with k1 " << *k1 << " and k2 " << *k2 << G4endl;
829 G4cerr
830 << "This may be because the energy of the incident particle is to high for the data table."
831 << G4endl;
832 G4cerr << "Particle energy (eV): " << k << G4endl;
833 exit(EXIT_FAILURE);
834 }
835
836 // We have a random number and we select the cumulated cross section data values surrounding our
837 // random number. But we need to do that for each T value (ie two T values) previously selected.
838 //
839 // First one.
840 auto cumulCS12 =
841 std::upper_bound(fProbaShellMap[materialID][p][ionizationLevelIndex][(*k1)].begin(),
842 fProbaShellMap[materialID][p][ionizationLevelIndex][(*k1)].end(), random);
843 auto cumulCS11 = cumulCS12 - 1;
844 // Second one.
845 auto cumulCS22 =
846 std::upper_bound(fProbaShellMap[materialID][p][ionizationLevelIndex][(*k2)].begin(),
847 fProbaShellMap[materialID][p][ionizationLevelIndex][(*k2)].end(), random);
848 auto cumulCS21 = cumulCS22 - 1;
849
850 // Now that we have the "values" through pointers, we access them.
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 electron energy
860 // *************************************************************
861
862 // Here we will get four E values corresponding to our four cumulated cross section values
863 // previously selected. But we need to take into account a specific case: we have selected a shell
864 // by using the ionisation cross section table and, since we get two T values, we could have
865 // differential cross sections (or cumulated) equal to 0 for the lower T and not for the upper T.
866 // When looking for the cumulated cross section values which surround the selected random number
867 // (for the lower T), the upper_bound method will only found 0 values. Thus, the upper_bound
868 // method will return the last E value present in the table for the selected T. The last E value
869 // being the highest, we will later perform an interpolation between a high E value (for the lower
870 // T) and a small E value (for the upper T). This is inconsistent because if the cross section are
871 // equal to zero for the lower T then it means it is not possible to ionize and, thus, to have a
872 // secondary electron. But, in our situation, it is possible to ionize for the upper T AND for an
873 // interpolate T value between Tupper Tlower. That's why the final E value should be interpolate
874 // between 0 and the E value (upper T).
875 //
876 if (cumulCS12 == fProbaShellMap[materialID][p][ionizationLevelIndex][(*k1)].end()) {
877 // Here we are in the special case and we force Elower1 and Elower2 to be equal at 0 for the
878 // interpolation.
879 secElecE11 = 0;
880 secElecE12 = 0;
881 secElecE21 = fEnergySecondaryData[materialID][p][ionizationLevelIndex][valueK2][valueCumulCS21];
882 secElecE22 = fEnergySecondaryData[materialID][p][ionizationLevelIndex][valueK2][valueCumulCS22];
883
884 valueCumulCS11 = 0;
885 valueCumulCS12 = 0;
886 }
887 else {
888 // No special case, interpolation will happen as usual.
889 secElecE11 = fEnergySecondaryData[materialID][p][ionizationLevelIndex][valueK1][valueCumulCS11];
890 secElecE12 = fEnergySecondaryData[materialID][p][ionizationLevelIndex][valueK1][valueCumulCS12];
891 secElecE21 = fEnergySecondaryData[materialID][p][ionizationLevelIndex][valueK2][valueCumulCS21];
892 secElecE22 = fEnergySecondaryData[materialID][p][ionizationLevelIndex][valueK2][valueCumulCS22];
893 }
894
895 ejectedElectronEnergy =
896 QuadInterpolator(valueCumulCS11, valueCumulCS12, valueCumulCS21, valueCumulCS22, secElecE11,
897 secElecE12, secElecE21, secElecE22, valueK1, valueK2, k, random);
898
899 // **********************************************
900 // Some tests for debugging
901 // **********************************************
902
903 G4double bindingEnergy(ptbStructure.IonisationEnergy(ionizationLevelIndex, materialID) / eV);
904 if (k - ejectedElectronEnergy - bindingEnergy <= 0 || ejectedElectronEnergy <= 0) {
905 G4cout << "k " << k << G4endl;
906 G4cout << "material ID : " << materialID << G4endl;
907 G4cout << "secondaryKin " << ejectedElectronEnergy << G4endl;
908 G4cout << "shell " << ionizationLevelIndex << G4endl;
909 G4cout << "bindingEnergy " << bindingEnergy << G4endl;
910 G4cout << "scatteredEnergy " << k - ejectedElectronEnergy - bindingEnergy << G4endl;
911 G4cout << "rand " << random << G4endl;
912 G4cout << "surrounding k values: valueK1 valueK2\n" << valueK1 << " " << valueK2 << G4endl;
913 G4cout << "surrounding E values: secElecE11 secElecE12 secElecE21 secElecE22\n"
914 << secElecE11 << " " << secElecE12 << " " << secElecE21 << " " << secElecE22 << " "
915 << G4endl;
916 G4cout
917 << "surrounding cumulCS values: valueCumulCS11 valueCumulCS12 valueCumulCS21 valueCumulCS22\n"
918 << valueCumulCS11 << " " << valueCumulCS12 << " " << valueCumulCS21 << " " << valueCumulCS22
919 << " " << G4endl;
920 G4ExceptionDescription errmsg;
921 errmsg << "*****************************" << G4endl;
922 errmsg << "Fatal error, EXIT." << G4endl;
923 G4Exception("G4DNAPTBIonisationModel::RandomizeEjectedElectronEnergyFromCumulated", "",
924 FatalException, errmsg);
925 exit(EXIT_FAILURE);
926 }
927
928 return ejectedElectronEnergy * eV;
929 }
930
931 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
932
933 G4double G4DNAPTBIonisationModel::LogLogInterpolate(G4double e1, G4double e2, G4double e,
934 G4double xs1, G4double xs2)
935 {
936 G4double value(0);
937
938 // Switch to log-lin interpolation for faster code
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 - e1) / (e2 - e1)));
944 }
945
946 // Switch to lin-lin interpolation for faster code
947 // in case one of xs1 or xs2 (=cum proba) value is zero
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 - e1));
953 }
954
955 return value;
956 }
957
958 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
959
960 G4double G4DNAPTBIonisationModel::QuadInterpolator(G4double e11, G4double e12, G4double e21,
961 G4double e22, G4double xs11, G4double xs12,
962 G4double xs21, G4double xs22, G4double t1,
963 G4double t2, G4double t, G4double e)
964 {
965 G4double interpolatedvalue1, interpolatedvalue2, value;
966 (xs11 != xs12) ? interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12)
967 : interpolatedvalue1 = xs11;
968
969 (xs21 != xs22) ? interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22)
970 : interpolatedvalue2 = xs21;
971
972 (interpolatedvalue1 != interpolatedvalue2)
973 ? value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2)
974 : value = interpolatedvalue1;
975 return value;
976 }
977