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