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