| 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 // 26 //
27 27
28 #include "G4DNARuddIonisationModel.hh" 28 #include "G4DNARuddIonisationModel.hh"
29 #include "G4PhysicalConstants.hh" 29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh"
31 #include "G4UAtomicDeexcitation.hh" 31 #include "G4UAtomicDeexcitation.hh"
32 #include "G4LossTableManager.hh" 32 #include "G4LossTableManager.hh"
33 #include "G4DNAChemistryManager.hh" 33 #include "G4DNAChemistryManager.hh"
34 #include "G4DNAMolecularMaterial.hh" 34 #include "G4DNAMolecularMaterial.hh"
35 #include "G4DNARuddAngle.hh" 35 #include "G4DNARuddAngle.hh"
36 #include "G4DeltaAngle.hh" 36 #include "G4DeltaAngle.hh"
37 #include "G4Exp.hh" 37 #include "G4Exp.hh"
38 #include "G4Pow.hh" 38 #include "G4Pow.hh"
39 #include "G4Log.hh" 39 #include "G4Log.hh"
40 #include "G4Alpha.hh" 40 #include "G4Alpha.hh"
41 41
42 static G4Pow * gpow = G4Pow::GetInstance(); 42 static G4Pow * gpow = G4Pow::GetInstance();
43 //....oooOO0OOooo........oooOO0OOooo........oo 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 44
45 using namespace std; 45 using namespace std;
46 46
47 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 48
49 G4DNARuddIonisationModel::G4DNARuddIonisationM 49 G4DNARuddIonisationModel::G4DNARuddIonisationModel(const G4ParticleDefinition*,
50 50 const G4String& nam) :
51 G4VEmModel(nam) 51 G4VEmModel(nam)
52 { 52 {
53 slaterEffectiveCharge[0] = 0.; 53 slaterEffectiveCharge[0] = 0.;
54 slaterEffectiveCharge[1] = 0.; 54 slaterEffectiveCharge[1] = 0.;
55 slaterEffectiveCharge[2] = 0.; 55 slaterEffectiveCharge[2] = 0.;
56 sCoefficient[0] = 0.; 56 sCoefficient[0] = 0.;
57 sCoefficient[1] = 0.; 57 sCoefficient[1] = 0.;
58 sCoefficient[2] = 0.; 58 sCoefficient[2] = 0.;
59 59
60 lowEnergyLimitForZ1 = 0 * eV; 60 lowEnergyLimitForZ1 = 0 * eV;
61 lowEnergyLimitForZ2 = 0 * eV; 61 lowEnergyLimitForZ2 = 0 * eV;
62 lowEnergyLimitOfModelForZ1 = 100 * eV; 62 lowEnergyLimitOfModelForZ1 = 100 * eV;
63 lowEnergyLimitOfModelForZ2 = 1 * keV; 63 lowEnergyLimitOfModelForZ2 = 1 * keV;
64 killBelowEnergyForZ1 = lowEnergyLimitOfModel 64 killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
65 killBelowEnergyForZ2 = lowEnergyLimitOfModel 65 killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
66 66
67 verboseLevel = 0; 67 verboseLevel = 0;
68 // Verbosity scale: 68 // Verbosity scale:
69 // 0 = nothing 69 // 0 = nothing
70 // 1 = warning for energy non-conservation 70 // 1 = warning for energy non-conservation
71 // 2 = details of energy budget 71 // 2 = details of energy budget
72 // 3 = calculation of cross sections, file o 72 // 3 = calculation of cross sections, file openings, sampling of atoms
73 // 4 = entering in methods 73 // 4 = entering in methods
74 74
75 if (verboseLevel > 0) 75 if (verboseLevel > 0)
76 { 76 {
77 G4cout << "Rudd ionisation model is constr 77 G4cout << "Rudd ionisation model is constructed " << G4endl;
78 } 78 }
79 79
80 // Define default angular generator 80 // Define default angular generator
81 SetAngularDistribution(new G4DNARuddAngle()) 81 SetAngularDistribution(new G4DNARuddAngle());
82 82
83 // Mark this model as "applicable" for atomi 83 // Mark this model as "applicable" for atomic deexcitation
84 SetDeexcitationFlag(true); 84 SetDeexcitationFlag(true);
85 85
86 // Selection of stationary mode 86 // Selection of stationary mode
87 87
88 statCode = false; 88 statCode = false;
89 } 89 }
90 90
91 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 92
93 G4DNARuddIonisationModel::~G4DNARuddIonisation 93 G4DNARuddIonisationModel::~G4DNARuddIonisationModel()
94 { 94 {
95 // Cross section 95 // Cross section
96 96
97 std::map<G4String, G4DNACrossSectionDataSet* 97 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
98 for (pos = tableData.begin(); pos != tableDa 98 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
99 { 99 {
100 G4DNACrossSectionDataSet* table = pos->sec 100 G4DNACrossSectionDataSet* table = pos->second;
101 delete table; 101 delete table;
102 } 102 }
103 103
104 // The following removal is forbidden since 104 // The following removal is forbidden since G4VEnergyLossmodel takes care of deletion
105 // Coverity however will signal this as an e 105 // Coverity however will signal this as an error
106 // if (fAtomDeexcitation) {delete fAtomDeex 106 // if (fAtomDeexcitation) {delete fAtomDeexcitation;}
107 107
108 } 108 }
109 109
110 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 111
112 void G4DNARuddIonisationModel::Initialise(cons 112 void G4DNARuddIonisationModel::Initialise(const G4ParticleDefinition* particle,
113 cons 113 const G4DataVector& /*cuts*/)
114 { 114 {
115 115
116 if (verboseLevel > 3) 116 if (verboseLevel > 3)
117 { 117 {
118 G4cout << "Calling G4DNARuddIonisationMode 118 G4cout << "Calling G4DNARuddIonisationModel::Initialise()" << G4endl;
119 } 119 }
120 120
121 // Energy limits 121 // Energy limits
122 122
123 G4String fileProton("dna/sigma_ionisation_p_ 123 G4String fileProton("dna/sigma_ionisation_p_rudd");
124 G4String fileHydrogen("dna/sigma_ionisation_ 124 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
125 G4String fileAlphaPlusPlus("dna/sigma_ionisa 125 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
126 G4String fileAlphaPlus("dna/sigma_ionisation 126 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
127 G4String fileHelium("dna/sigma_ionisation_he 127 G4String fileHelium("dna/sigma_ionisation_he_rudd");
128 128
129 G4DNAGenericIonsManager *instance; 129 G4DNAGenericIonsManager *instance;
130 instance = G4DNAGenericIonsManager::Instance 130 instance = G4DNAGenericIonsManager::Instance();
131 protonDef = G4Proton::ProtonDefinition(); 131 protonDef = G4Proton::ProtonDefinition();
132 hydrogenDef = instance->GetIon("hydrogen"); 132 hydrogenDef = instance->GetIon("hydrogen");
133 alphaPlusPlusDef = G4Alpha::Alpha(); 133 alphaPlusPlusDef = G4Alpha::Alpha();
134 alphaPlusDef = instance->GetIon("alpha+"); 134 alphaPlusDef = instance->GetIon("alpha+");
135 heliumDef = instance->GetIon("helium"); 135 heliumDef = instance->GetIon("helium");
136 136
137 G4String proton; 137 G4String proton;
138 G4String hydrogen; 138 G4String hydrogen;
139 G4String alphaPlusPlus; 139 G4String alphaPlusPlus;
140 G4String alphaPlus; 140 G4String alphaPlus;
141 G4String helium; 141 G4String helium;
142 142
143 G4double scaleFactor = 1 * m*m; 143 G4double scaleFactor = 1 * m*m;
144 144
145 // LIMITS AND DATA 145 // LIMITS AND DATA
146 146
147 // ***************************************** 147 // ********************************************************
148 148
149 proton = protonDef->GetParticleName(); 149 proton = protonDef->GetParticleName();
150 tableFile[proton] = fileProton; 150 tableFile[proton] = fileProton;
151 151
152 lowEnergyLimit[proton] = lowEnergyLimitForZ1 152 lowEnergyLimit[proton] = lowEnergyLimitForZ1;
153 highEnergyLimit[proton] = 500. * keV; 153 highEnergyLimit[proton] = 500. * keV;
154 154
155 // Cross section 155 // Cross section
156 156
157 auto tableProton = new G4DNACrossSectionDat 157 auto tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
158 eV, 158 eV,
159 scaleFactor ); 159 scaleFactor );
160 tableProton->LoadData(fileProton); 160 tableProton->LoadData(fileProton);
161 tableData[proton] = tableProton; 161 tableData[proton] = tableProton;
162 162
163 // ***************************************** 163 // ********************************************************
164 164
165 hydrogen = hydrogenDef->GetParticleName(); 165 hydrogen = hydrogenDef->GetParticleName();
166 tableFile[hydrogen] = fileHydrogen; 166 tableFile[hydrogen] = fileHydrogen;
167 167
168 lowEnergyLimit[hydrogen] = lowEnergyLimitFor 168 lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
169 highEnergyLimit[hydrogen] = 100. * MeV; 169 highEnergyLimit[hydrogen] = 100. * MeV;
170 170
171 // Cross section 171 // Cross section
172 172
173 auto tableHydrogen = new G4DNACrossSectionD 173 auto tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
174 eV, 174 eV,
175 scaleFactor ); 175 scaleFactor );
176 tableHydrogen->LoadData(fileHydrogen); 176 tableHydrogen->LoadData(fileHydrogen);
177 177
178 tableData[hydrogen] = tableHydrogen; 178 tableData[hydrogen] = tableHydrogen;
179 179
180 // ***************************************** 180 // ********************************************************
181 181
182 alphaPlusPlus = alphaPlusPlusDef->GetParticl 182 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
183 tableFile[alphaPlusPlus] = fileAlphaPlusPlus 183 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
184 184
185 lowEnergyLimit[alphaPlusPlus] = lowEnergyLim 185 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
186 highEnergyLimit[alphaPlusPlus] = 400. * MeV; 186 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
187 187
188 // Cross section 188 // Cross section
189 189
190 auto tableAlphaPlusPlus = new G4DNACrossSec 190 auto tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
191 eV, 191 eV,
192 scaleFactor ); 192 scaleFactor );
193 tableAlphaPlusPlus->LoadData(fileAlphaPlusPl 193 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
194 194
195 tableData[alphaPlusPlus] = tableAlphaPlusPlu 195 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
196 196
197 // ***************************************** 197 // ********************************************************
198 198
199 alphaPlus = alphaPlusDef->GetParticleName(); 199 alphaPlus = alphaPlusDef->GetParticleName();
200 tableFile[alphaPlus] = fileAlphaPlus; 200 tableFile[alphaPlus] = fileAlphaPlus;
201 201
202 lowEnergyLimit[alphaPlus] = lowEnergyLimitFo 202 lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
203 highEnergyLimit[alphaPlus] = 400. * MeV; 203 highEnergyLimit[alphaPlus] = 400. * MeV;
204 204
205 // Cross section 205 // Cross section
206 206
207 auto tableAlphaPlus = new G4DNACrossSection 207 auto tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
208 eV, 208 eV,
209 scaleFactor ); 209 scaleFactor );
210 tableAlphaPlus->LoadData(fileAlphaPlus); 210 tableAlphaPlus->LoadData(fileAlphaPlus);
211 tableData[alphaPlus] = tableAlphaPlus; 211 tableData[alphaPlus] = tableAlphaPlus;
212 212
213 // ***************************************** 213 // ********************************************************
214 214
215 helium = heliumDef->GetParticleName(); 215 helium = heliumDef->GetParticleName();
216 tableFile[helium] = fileHelium; 216 tableFile[helium] = fileHelium;
217 217
218 lowEnergyLimit[helium] = lowEnergyLimitForZ2 218 lowEnergyLimit[helium] = lowEnergyLimitForZ2;
219 highEnergyLimit[helium] = 400. * MeV; 219 highEnergyLimit[helium] = 400. * MeV;
220 220
221 // Cross section 221 // Cross section
222 222
223 auto tableHelium = new G4DNACrossSectionDat 223 auto tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
224 eV, 224 eV,
225 scaleFactor ); 225 scaleFactor );
226 tableHelium->LoadData(fileHelium); 226 tableHelium->LoadData(fileHelium);
227 tableData[helium] = tableHelium; 227 tableData[helium] = tableHelium;
228 228
229 // 229 //
230 230
231 if (particle==protonDef) 231 if (particle==protonDef)
232 { 232 {
233 SetLowEnergyLimit(lowEnergyLimit[proton]); 233 SetLowEnergyLimit(lowEnergyLimit[proton]);
234 SetHighEnergyLimit(highEnergyLimit[proton] 234 SetHighEnergyLimit(highEnergyLimit[proton]);
235 } 235 }
236 236
237 if (particle==hydrogenDef) 237 if (particle==hydrogenDef)
238 { 238 {
239 SetLowEnergyLimit(lowEnergyLimit[hydrogen] 239 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
240 SetHighEnergyLimit(highEnergyLimit[hydroge 240 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
241 } 241 }
242 242
243 if (particle==heliumDef) 243 if (particle==heliumDef)
244 { 244 {
245 SetLowEnergyLimit(lowEnergyLimit[helium]); 245 SetLowEnergyLimit(lowEnergyLimit[helium]);
246 SetHighEnergyLimit(highEnergyLimit[helium] 246 SetHighEnergyLimit(highEnergyLimit[helium]);
247 } 247 }
248 248
249 if (particle==alphaPlusDef) 249 if (particle==alphaPlusDef)
250 { 250 {
251 SetLowEnergyLimit(lowEnergyLimit[alphaPlus 251 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
252 SetHighEnergyLimit(highEnergyLimit[alphaPl 252 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
253 } 253 }
254 254
255 if (particle==alphaPlusPlusDef) 255 if (particle==alphaPlusPlusDef)
256 { 256 {
257 SetLowEnergyLimit(lowEnergyLimit[alphaPlus 257 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
258 SetHighEnergyLimit(highEnergyLimit[alphaPl 258 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
259 } 259 }
260 260
261 if( verboseLevel>0 ) 261 if( verboseLevel>0 )
262 { 262 {
263 G4cout << "Rudd ionisation model is initia 263 G4cout << "Rudd ionisation model is initialized " << G4endl
264 << "Energy range: " 264 << "Energy range: "
265 << LowEnergyLimit() / eV << " eV - " 265 << LowEnergyLimit() / eV << " eV - "
266 << HighEnergyLimit() / keV << " keV for " 266 << HighEnergyLimit() / keV << " keV for "
267 << particle->GetParticleName() 267 << particle->GetParticleName()
268 << G4endl; 268 << G4endl;
269 } 269 }
270 270
271 // Initialize water density pointer 271 // Initialize water density pointer
272 fpWaterDensity = G4DNAMolecularMaterial::Ins 272 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
273 273
274 // 274 //
275 275
276 fAtomDeexcitation = G4LossTableManager::Inst 276 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
277 277
278 if (isInitialised) 278 if (isInitialised)
279 { return;} 279 { return;}
280 fParticleChangeForGamma = GetParticleChangeF 280 fParticleChangeForGamma = GetParticleChangeForGamma();
281 isInitialised = true; 281 isInitialised = true;
282 282
283 } 283 }
284 284
285 //....oooOO0OOooo........oooOO0OOooo........oo 285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286 286
287 G4double G4DNARuddIonisationModel::CrossSectio 287 G4double G4DNARuddIonisationModel::CrossSectionPerVolume(const G4Material* material,
288 288 const G4ParticleDefinition* particleDefinition,
289 289 G4double k,
290 290 G4double,
291 291 G4double)
292 { 292 {
293 if (verboseLevel > 3) 293 if (verboseLevel > 3)
294 { 294 {
295 G4cout << "Calling CrossSectionPerVolume() 295 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationModel"
296 << G4endl; 296 << G4endl;
297 } 297 }
298 298
299 // Calculate total cross section for model 299 // Calculate total cross section for model
300 300
301 if ( 301 if (
302 particleDefinition != protonDef 302 particleDefinition != protonDef
303 && 303 &&
304 particleDefinition != hydrogenDef 304 particleDefinition != hydrogenDef
305 && 305 &&
306 particleDefinition != alphaPlusPlusDef 306 particleDefinition != alphaPlusPlusDef
307 && 307 &&
308 particleDefinition != alphaPlusDef 308 particleDefinition != alphaPlusDef
309 && 309 &&
310 particleDefinition != heliumDef 310 particleDefinition != heliumDef
311 ) 311 )
312 312
313 return 0; 313 return 0;
314 314
315 G4double lowLim = 0; 315 G4double lowLim = 0;
316 316
317 if ( particleDefinition == protonDef 317 if ( particleDefinition == protonDef
318 || particleDefinition == hydrogenDef 318 || particleDefinition == hydrogenDef
319 ) 319 )
320 320
321 lowLim = lowEnergyLimitOfModelForZ1; 321 lowLim = lowEnergyLimitOfModelForZ1;
322 322
323 if ( particleDefinition == alphaPlusPlusDef 323 if ( particleDefinition == alphaPlusPlusDef
324 || particleDefinition == alphaPlusDef 324 || particleDefinition == alphaPlusDef
325 || particleDefinition == heliumDef 325 || particleDefinition == heliumDef
326 ) 326 )
327 327
328 lowLim = lowEnergyLimitOfModelForZ2; 328 lowLim = lowEnergyLimitOfModelForZ2;
329 329
330 G4double highLim = 0; 330 G4double highLim = 0;
331 G4double sigma=0; 331 G4double sigma=0;
332 332
333 G4double waterDensity = (*fpWaterDensity)[ma 333 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
334 334
335 const G4String& particleName = particleDefin 335 const G4String& particleName = particleDefinition->GetParticleName();
336 336
337 // SI - the following is useless since lowLi 337 // SI - the following is useless since lowLim is already defined
338 /* 338 /*
339 std::map< G4String,G4double,std::less<G4Stri 339 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
340 pos1 = lowEnergyLimit.find(particleName); 340 pos1 = lowEnergyLimit.find(particleName);
341 341
342 if (pos1 != lowEnergyLimit.end()) 342 if (pos1 != lowEnergyLimit.end())
343 { 343 {
344 lowLim = pos1->second; 344 lowLim = pos1->second;
345 } 345 }
346 */ 346 */
347 347
348 std::map< G4String,G4double,std::less<G4Stri 348 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
349 pos2 = highEnergyLimit.find(particleName); 349 pos2 = highEnergyLimit.find(particleName);
350 350
351 if (pos2 != highEnergyLimit.end()) 351 if (pos2 != highEnergyLimit.end())
352 { 352 {
353 highLim = pos2->second; 353 highLim = pos2->second;
354 } 354 }
355 355
356 if (k <= highLim) 356 if (k <= highLim)
357 { 357 {
358 //SI : XS must not be zero otherwise sampl 358 //SI : XS must not be zero otherwise sampling of secondaries method ignored
359 359
360 if (k < lowLim) k = lowLim; 360 if (k < lowLim) k = lowLim;
361 361
362 // 362 //
363 363
364 std::map< G4String,G4DNACrossSectionDataSe 364 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
365 pos = tableData.find(particleName); 365 pos = tableData.find(particleName);
366 366
367 if (pos != tableData.end()) 367 if (pos != tableData.end())
368 { 368 {
369 G4DNACrossSectionDataSet* table = pos->s 369 G4DNACrossSectionDataSet* table = pos->second;
370 if (table != nullptr) 370 if (table != nullptr)
371 { 371 {
372 sigma = table->FindValue(k); 372 sigma = table->FindValue(k);
373 } 373 }
374 } 374 }
375 else 375 else
376 { 376 {
377 G4Exception("G4DNARuddIonisationModel::C 377 G4Exception("G4DNARuddIonisationModel::CrossSectionPerVolume","em0002",
378 FatalException,"Model not applicable 378 FatalException,"Model not applicable to particle type.");
379 } 379 }
380 380
381 } 381 }
382 382
383 if (verboseLevel > 2) 383 if (verboseLevel > 2)
384 { 384 {
385 G4cout << "_______________________________ 385 G4cout << "__________________________________" << G4endl;
386 G4cout << "G4DNARuddIonisationModel - XS I 386 G4cout << "G4DNARuddIonisationModel - XS INFO START" << G4endl;
387 G4cout << "Kinetic energy(eV)=" << k/eV << 387 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
388 G4cout << "Cross section per water molecul 388 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
389 G4cout << "Cross section per water molecul 389 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
390 //G4cout << " - Cross section per water mo 390 //G4cout << " - Cross section per water molecule (cm^-1)="
391 //<< sigma*material->GetAtomicNumDensityVe 391 //<< sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
392 G4cout << "G4DNARuddIonisationModel - XS I 392 G4cout << "G4DNARuddIonisationModel - XS INFO END" << G4endl;
393 } 393 }
394 394
395 return sigma*waterDensity; 395 return sigma*waterDensity;
396 396
397 } 397 }
398 398
399 //....oooOO0OOooo........oooOO0OOooo........oo 399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 400
401 void G4DNARuddIonisationModel::SampleSecondari 401 void G4DNARuddIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
402 402 const G4MaterialCutsCouple* couple,
403 403 const G4DynamicParticle* particle,
404 404 G4double,
405 405 G4double)
406 { 406 {
407 if (verboseLevel > 3) 407 if (verboseLevel > 3)
408 { 408 {
409 G4cout << "Calling SampleSecondaries() of 409 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationModel"
410 << G4endl; 410 << G4endl;
411 } 411 }
412 412
413 G4double lowLim = 0; 413 G4double lowLim = 0;
414 G4double highLim = 0; 414 G4double highLim = 0;
415 415
416 if ( particle->GetDefinition() == protonDef 416 if ( particle->GetDefinition() == protonDef
417 || particle->GetDefinition() == hydrogen 417 || particle->GetDefinition() == hydrogenDef
418 ) 418 )
419 419
420 lowLim = killBelowEnergyForZ1; 420 lowLim = killBelowEnergyForZ1;
421 421
422 if ( particle->GetDefinition() == alphaPlusP 422 if ( particle->GetDefinition() == alphaPlusPlusDef
423 || particle->GetDefinition() == alphaPlu 423 || particle->GetDefinition() == alphaPlusDef
424 || particle->GetDefinition() == heliumDe 424 || particle->GetDefinition() == heliumDef
425 ) 425 )
426 426
427 lowLim = killBelowEnergyForZ2; 427 lowLim = killBelowEnergyForZ2;
428 428
429 G4double k = particle->GetKineticEnergy(); 429 G4double k = particle->GetKineticEnergy();
430 430
431 const G4String& particleName = particle->Get 431 const G4String& particleName = particle->GetDefinition()->GetParticleName();
432 432
433 // SI - the following is useless since lowLi 433 // SI - the following is useless since lowLim is already defined
434 /* 434 /*
435 std::map< G4String,G4double,std::less<G4Str 435 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
436 pos1 = lowEnergyLimit.find(particleName); 436 pos1 = lowEnergyLimit.find(particleName);
437 437
438 if (pos1 != lowEnergyLimit.end()) 438 if (pos1 != lowEnergyLimit.end())
439 { 439 {
440 lowLim = pos1->second; 440 lowLim = pos1->second;
441 } 441 }
442 */ 442 */
443 443
444 std::map< G4String,G4double,std::less<G4Stri 444 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
445 pos2 = highEnergyLimit.find(particleName); 445 pos2 = highEnergyLimit.find(particleName);
446 446
447 if (pos2 != highEnergyLimit.end()) 447 if (pos2 != highEnergyLimit.end())
448 { 448 {
449 highLim = pos2->second; 449 highLim = pos2->second;
450 } 450 }
451 451
452 if (k >= lowLim && k <= highLim) 452 if (k >= lowLim && k <= highLim)
453 { 453 {
454 G4ParticleDefinition* definition = particl 454 G4ParticleDefinition* definition = particle->GetDefinition();
455 G4ParticleMomentum primaryDirection = part 455 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
456 /* 456 /*
457 G4double particleMass = definition->GetPD 457 G4double particleMass = definition->GetPDGMass();
458 G4double totalEnergy = k + particleMass; 458 G4double totalEnergy = k + particleMass;
459 G4double pSquare = k*(totalEnergy+particl 459 G4double pSquare = k*(totalEnergy+particleMass);
460 G4double totalMomentum = std::sqrt(pSquar 460 G4double totalMomentum = std::sqrt(pSquare);
461 */ 461 */
462 462
463 G4int ionizationShell = RandomSelect(k,par 463 G4int ionizationShell = RandomSelect(k,particleName);
464 464
465 G4double bindingEnergy = 0; 465 G4double bindingEnergy = 0;
466 bindingEnergy = waterStructure.IonisationE 466 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
467 467
468 //SI: additional protection if tcs interpo 468 //SI: additional protection if tcs interpolation method is modified
469 if (k<bindingEnergy) return; 469 if (k<bindingEnergy) return;
470 // 470 //
471 471
472 // SI - For atom. deexc. tagging - 23/05/2 472 // SI - For atom. deexc. tagging - 23/05/2017
473 G4int Z = 8; 473 G4int Z = 8;
474 // 474 //
475 475
476 G4double secondaryKinetic = RandomizeEject 476 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
477 477
478 G4ThreeVector deltaDirection = 478 G4ThreeVector deltaDirection =
479 GetAngularDistribution()->SampleDirectionF 479 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
480 Z, ionizationShell, 480 Z, ionizationShell,
481 couple->GetMaterial()); 481 couple->GetMaterial());
482 482
483 auto dp = new G4DynamicParticle (G4Electr 483 auto dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
484 fvect->push_back(dp); 484 fvect->push_back(dp);
485 485
486 // Ignored for ions on electrons 486 // Ignored for ions on electrons
487 /* 487 /*
488 G4double deltaTotalMomentum = std::sqrt(s 488 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
489 489
490 G4double finalPx = totalMomentum*primaryD 490 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
491 G4double finalPy = totalMomentum*primaryD 491 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
492 G4double finalPz = totalMomentum*primaryD 492 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
493 G4double finalMomentum = std::sqrt(finalP 493 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
494 finalPx /= finalMomentum; 494 finalPx /= finalMomentum;
495 finalPy /= finalMomentum; 495 finalPy /= finalMomentum;
496 finalPz /= finalMomentum; 496 finalPz /= finalMomentum;
497 497
498 G4ThreeVector direction; 498 G4ThreeVector direction;
499 direction.set(finalPx,finalPy,finalPz); 499 direction.set(finalPx,finalPy,finalPz);
500 500
501 fParticleChangeForGamma->ProposeMomentumD 501 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
502 */ 502 */
503 503
504 fParticleChangeForGamma->ProposeMomentumDi 504 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
505 505
506 // sample deexcitation 506 // sample deexcitation
507 // here we assume that H_{2}O electronic l 507 // here we assume that H_{2}O electronic levels are the same of Oxigen.
508 // this can be considered true with a roug 508 // this can be considered true with a rough 10% error in energy on K-shell,
509 509
510 size_t secNumberInit = 0;// need to know a 510 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
511 size_t secNumberFinal = 0;// So I'll make 511 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
512 512
513 G4double scatteredEnergy = k-bindingEnergy 513 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
514 514
515 // SI: only atomic deexcitation from K she 515 // SI: only atomic deexcitation from K shell is considered
516 if((fAtomDeexcitation != nullptr) && ioniz 516 if((fAtomDeexcitation != nullptr) && ionizationShell == 4)
517 { 517 {
518 const G4AtomicShell* shell 518 const G4AtomicShell* shell
519 = fAtomDeexcitation->GetAtomicShell(Z, 519 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
520 secNumberInit = fvect->size(); 520 secNumberInit = fvect->size();
521 fAtomDeexcitation->GenerateParticles(fve 521 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
522 secNumberFinal = fvect->size(); 522 secNumberFinal = fvect->size();
523 523
524 if(secNumberFinal > secNumberInit) 524 if(secNumberFinal > secNumberInit)
525 { 525 {
526 for (size_t i=secNumberInit; i<secNumberFina 526 for (size_t i=secNumberInit; i<secNumberFinal; ++i)
527 { 527 {
528 //Check if there is enough residual 528 //Check if there is enough residual energy
529 if (bindingEnergy >= ((*fvect)[i])-> 529 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
530 { 530 {
531 //Ok, this is a valid secondary: 531 //Ok, this is a valid secondary: keep it
532 bindingEnergy -= ((*fvect)[i])->GetKine 532 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
533 } 533 }
534 else 534 else
535 { 535 {
536 //Invalid secondary: not enough energy 536 //Invalid secondary: not enough energy to create it!
537 //Keep its energy in the local deposit 537 //Keep its energy in the local deposit
538 delete (*fvect)[i]; 538 delete (*fvect)[i];
539 (*fvect)[i]=nullptr; 539 (*fvect)[i]=nullptr;
540 } 540 }
541 } 541 }
542 } 542 }
543 543
544 } 544 }
545 545
546 //This should never happen 546 //This should never happen
547 if(bindingEnergy < 0.0) 547 if(bindingEnergy < 0.0)
548 G4Exception("G4DNAEmfietzoglouIonisatioMo 548 G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
549 "em2050",FatalException,"Nega 549 "em2050",FatalException,"Negative local energy deposit");
550 550
551 //bindingEnergy has been decreased 551 //bindingEnergy has been decreased
552 //by the amount of energy taken away by de 552 //by the amount of energy taken away by deexc. products
553 if (!statCode) 553 if (!statCode)
554 { 554 {
555 fParticleChangeForGamma->SetProposedKine 555 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
556 fParticleChangeForGamma->ProposeLocalEne 556 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
557 } 557 }
558 else 558 else
559 { 559 {
560 fParticleChangeForGamma->SetProposedKine 560 fParticleChangeForGamma->SetProposedKineticEnergy(k);
561 fParticleChangeForGamma->ProposeLocalEne 561 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
562 } 562 }
563 563
564 // debug 564 // debug
565 // k-scatteredEnergy-secondaryKinetic-deex 565 // k-scatteredEnergy-secondaryKinetic-deexSecEnergy = k-(k-bindingEnergy-secondaryKinetic)-secondaryKinetic-deexSecEnergy =
566 // = k-k+bindingEnergy+secondaryKinetic-se 566 // = k-k+bindingEnergy+secondaryKinetic-secondaryKinetic-deexSecEnergy=
567 // = bindingEnergy-deexSecEnergy 567 // = bindingEnergy-deexSecEnergy
568 // SO deexSecEnergy=0 => LocalEnergyDeposi 568 // SO deexSecEnergy=0 => LocalEnergyDeposit = bindingEnergy
569 569
570 // TEST ////////////////////////// 570 // TEST //////////////////////////
571 // if (secondaryKinetic<0) abort(); 571 // if (secondaryKinetic<0) abort();
572 // if (scatteredEnergy<0) abort(); 572 // if (scatteredEnergy<0) abort();
573 // if (k-scatteredEnergy-secondaryKinetic- 573 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
574 // if (k-scatteredEnergy<0) abort(); 574 // if (k-scatteredEnergy<0) abort();
575 ///////////////////////////////// 575 /////////////////////////////////
576 576
577 const G4Track * theIncomingTrack = fPartic 577 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
578 G4DNAChemistryManager::Instance()->CreateW 578 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
579 ionizationShell, 579 ionizationShell,
580 theIncomingTrack); 580 theIncomingTrack);
581 } 581 }
582 582
583 // SI - not useful since low energy of model 583 // SI - not useful since low energy of model is 0 eV
584 584
585 if (k < lowLim) 585 if (k < lowLim)
586 { 586 {
587 fParticleChangeForGamma->SetProposedKineti 587 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
588 fParticleChangeForGamma->ProposeTrackStatu 588 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
589 fParticleChangeForGamma->ProposeLocalEnerg 589 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
590 } 590 }
591 591
592 } 592 }
593 593
594 //....oooOO0OOooo........oooOO0OOooo........oo 594 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
595 595
596 G4double G4DNARuddIonisationModel::RandomizeEj 596 G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
597 597 G4double k,
598 598 G4int shell)
599 { 599 {
600 G4double maximumKineticEnergyTransfer = 0.; 600 G4double maximumKineticEnergyTransfer = 0.;
601 601
602 if (particleDefinition == protonDef 602 if (particleDefinition == protonDef
603 || particleDefinition == hydrogenDef) 603 || particleDefinition == hydrogenDef)
604 { 604 {
605 maximumKineticEnergyTransfer = 4. * (elect 605 maximumKineticEnergyTransfer = 4. * (electron_mass_c2 / proton_mass_c2) * k;
606 } 606 }
607 607
608 else if (particleDefinition == heliumDef 608 else if (particleDefinition == heliumDef
609 || particleDefinition == alphaPlusDef 609 || particleDefinition == alphaPlusDef
610 || particleDefinition == alphaPlusPlusDe 610 || particleDefinition == alphaPlusPlusDef)
611 { 611 {
612 maximumKineticEnergyTransfer = 4. * (0.511 612 maximumKineticEnergyTransfer = 4. * (0.511 / 3728) * k;
613 } 613 }
614 614
615 G4double crossSectionMaximum = 0.; 615 G4double crossSectionMaximum = 0.;
616 616
617 for (G4double value = waterStructure.Ionisat 617 for (G4double value = waterStructure.IonisationEnergy(shell);
618 value <= 5. * waterStructure.IonisationE 618 value <= 5. * waterStructure.IonisationEnergy(shell) && k >= value;
619 value += 0.1 * eV) 619 value += 0.1 * eV)
620 { 620 {
621 G4double differentialCrossSection = 621 G4double differentialCrossSection =
622 DifferentialCrossSection(particleDefin 622 DifferentialCrossSection(particleDefinition, k, value, shell);
623 if (differentialCrossSection >= crossSecti 623 if (differentialCrossSection >= crossSectionMaximum)
624 crossSectionMaximum = differentialCrossS 624 crossSectionMaximum = differentialCrossSection;
625 } 625 }
626 626
627 G4double secElecKinetic = 0.; 627 G4double secElecKinetic = 0.;
628 628
629 do 629 do
630 { 630 {
631 secElecKinetic = G4UniformRand()* maximumK 631 secElecKinetic = G4UniformRand()* maximumKineticEnergyTransfer;
632 }while(G4UniformRand()*crossSectionMaximum > 632 }while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
633 k, 633 k,
634 secElecKinetic+waterStructure.Ionisa 634 secElecKinetic+waterStructure.IonisationEnergy(shell),
635 shell)); 635 shell));
636 636
637 return (secElecKinetic); 637 return (secElecKinetic);
638 } 638 }
639 639
640 //....oooOO0OOooo........oooOO0OOooo........oo 640 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
641 641
642 // The following section is not used anymore b 642 // The following section is not used anymore but is kept for memory
643 // GetAngularDistribution()->SampleDirectionFo 643 // GetAngularDistribution()->SampleDirectionForShell is used instead
644 644
645 /* 645 /*
646 void G4DNARuddIonisationModel::RandomizeEject 646 void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
647 G4double k, 647 G4double k,
648 G4double secKinetic, 648 G4double secKinetic,
649 G4double & cosTheta, 649 G4double & cosTheta,
650 G4double & phi ) 650 G4double & phi )
651 { 651 {
652 G4DNAGenericIonsManager *instance; 652 G4DNAGenericIonsManager *instance;
653 instance = G4DNAGenericIonsManager::Instance( 653 instance = G4DNAGenericIonsManager::Instance();
654 654
655 G4double maxSecKinetic = 0.; 655 G4double maxSecKinetic = 0.;
656 656
657 if (particleDefinition == protonDef 657 if (particleDefinition == protonDef
658 || particleDefinition == hydrogenDef) 658 || particleDefinition == hydrogenDef)
659 { 659 {
660 maxSecKinetic = 4.* (electron_mass_c2 / proto 660 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
661 } 661 }
662 662
663 else if (particleDefinition == heliumDef 663 else if (particleDefinition == heliumDef
664 || particleDefinition == alphaPlusDef 664 || particleDefinition == alphaPlusDef
665 || particleDefinition == alphaPlusPlusDef) 665 || particleDefinition == alphaPlusPlusDef)
666 { 666 {
667 maxSecKinetic = 4.* (0.511 / 3728) * k; 667 maxSecKinetic = 4.* (0.511 / 3728) * k;
668 } 668 }
669 669
670 phi = twopi * G4UniformRand(); 670 phi = twopi * G4UniformRand();
671 671
672 // cosTheta = std::sqrt(secKinetic / maxSecKi 672 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
673 673
674 // Restriction below 100 eV from Emfietzoglou 674 // Restriction below 100 eV from Emfietzoglou (2000)
675 675
676 if (secKinetic>100*eV) cosTheta = std::sqrt(s 676 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
677 else cosTheta = (2.*G4UniformRand())-1.; 677 else cosTheta = (2.*G4UniformRand())-1.;
678 678
679 } 679 }
680 */ 680 */
681 //....oooOO0OOooo........oooOO0OOooo........oo 681 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
682 G4double G4DNARuddIonisationModel::Differentia 682 G4double G4DNARuddIonisationModel::DifferentialCrossSection(G4ParticleDefinition* particleDefinition,
683 683 G4double k,
684 684 G4double energyTransfer,
685 685 G4int ionizationLevelIndex)
686 { 686 {
687 // Shells ids are 0 1 2 3 4 (4 is k shell) 687 // Shells ids are 0 1 2 3 4 (4 is k shell)
688 // !!Attention, "energyTransfer" here is the 688 // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
689 // that the secondary kinetic en 689 // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
690 // 690 //
691 // ds S F1(nu) + 691 // ds S F1(nu) + w * F2(nu)
692 // ---- = G(k) * ---- ----------------- 692 // ---- = G(k) * ---- -------------------------------------------
693 // dw Bj (1+w)^3 * [1 + e 693 // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
694 // 694 //
695 // w is the secondary electron kinetic Energ 695 // w is the secondary electron kinetic Energy in eV
696 // 696 //
697 // All the other parameters can be found in 697 // All the other parameters can be found in Rudd's Papers
698 // 698 //
699 // M.Eugene Rudd, 1988, User-Friendly model 699 // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
700 // electrons from protons or electron collis 700 // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
701 // 701 //
702 702
703 const G4int j = ionizationLevelIndex; 703 const G4int j = ionizationLevelIndex;
704 704
705 G4double A1; 705 G4double A1;
706 G4double B1; 706 G4double B1;
707 G4double C1; 707 G4double C1;
708 G4double D1; 708 G4double D1;
709 G4double E1; 709 G4double E1;
710 G4double A2; 710 G4double A2;
711 G4double B2; 711 G4double B2;
712 G4double C2; 712 G4double C2;
713 G4double D2; 713 G4double D2;
714 G4double alphaConst; 714 G4double alphaConst;
715 715
716 // const G4double Bj[5] = {12.61*eV, 14.73* 716 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
717 // The following values are provided by M. d 717 // The following values are provided by M. dingfelder (priv. comm)
718 const G4double Bj[5] = { 12.60 * eV, 14.70 * 718 const G4double Bj[5] = { 12.60 * eV, 14.70 * eV, 18.40 * eV, 32.20 * eV, 540
719 * eV }; 719 * eV };
720 720
721 if (j == 4) 721 if (j == 4)
722 { 722 {
723 //Data For Liquid Water K SHELL from Dingf 723 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
724 A1 = 1.25; 724 A1 = 1.25;
725 B1 = 0.5; 725 B1 = 0.5;
726 C1 = 1.00; 726 C1 = 1.00;
727 D1 = 1.00; 727 D1 = 1.00;
728 E1 = 3.00; 728 E1 = 3.00;
729 A2 = 1.10; 729 A2 = 1.10;
730 B2 = 1.30; 730 B2 = 1.30;
731 C2 = 1.00; 731 C2 = 1.00;
732 D2 = 0.00; 732 D2 = 0.00;
733 alphaConst = 0.66; 733 alphaConst = 0.66;
734 } else 734 } else
735 { 735 {
736 //Data For Liquid Water from Dingfelder (P 736 //Data For Liquid Water from Dingfelder (Protons in Water)
737 A1 = 1.02; 737 A1 = 1.02;
738 B1 = 82.0; 738 B1 = 82.0;
739 C1 = 0.45; 739 C1 = 0.45;
740 D1 = -0.80; 740 D1 = -0.80;
741 E1 = 0.38; 741 E1 = 0.38;
742 A2 = 1.07; 742 A2 = 1.07;
743 // Value provided by M. Dingfelder (priv. 743 // Value provided by M. Dingfelder (priv. comm)
744 B2 = 11.6; 744 B2 = 11.6;
745 // 745 //
746 C2 = 0.60; 746 C2 = 0.60;
747 D2 = 0.04; 747 D2 = 0.04;
748 alphaConst = 0.64; 748 alphaConst = 0.64;
749 } 749 }
750 750
751 const G4double n = 2.; 751 const G4double n = 2.;
752 const G4double Gj[5] = { 0.99, 1.11, 1.11, 0 752 const G4double Gj[5] = { 0.99, 1.11, 1.11, 0.52, 1. };
753 753
754 G4double wBig = (energyTransfer 754 G4double wBig = (energyTransfer
755 - waterStructure.IonisationEnergy(ioniza 755 - waterStructure.IonisationEnergy(ionizationLevelIndex));
756 if (wBig < 0) 756 if (wBig < 0)
757 return 0.; 757 return 0.;
758 758
759 G4double w = wBig / Bj[ionizationLevelIndex] 759 G4double w = wBig / Bj[ionizationLevelIndex];
760 // Note that the following (j==4) cases are 760 // Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
761 if (j == 4) 761 if (j == 4)
762 w = wBig / waterStructure.IonisationEnergy 762 w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
763 763
764 G4double Ry = 13.6 * eV; 764 G4double Ry = 13.6 * eV;
765 765
766 G4double tau = 0.; 766 G4double tau = 0.;
767 767
768 G4bool isProtonOrHydrogen = false; 768 G4bool isProtonOrHydrogen = false;
769 G4bool isHelium = false; 769 G4bool isHelium = false;
770 770
771 if (particleDefinition == protonDef 771 if (particleDefinition == protonDef
772 || particleDefinition == hydrogenDef) 772 || particleDefinition == hydrogenDef)
773 { 773 {
774 isProtonOrHydrogen = true; 774 isProtonOrHydrogen = true;
775 tau = (electron_mass_c2 / proton_mass_c2) 775 tau = (electron_mass_c2 / proton_mass_c2) * k;
776 } 776 }
777 777
778 else if (particleDefinition == heliumDef 778 else if (particleDefinition == heliumDef
779 || particleDefinition == alphaPlusDef 779 || particleDefinition == alphaPlusDef
780 || particleDefinition == alphaPlusPlusDe 780 || particleDefinition == alphaPlusPlusDef)
781 { 781 {
782 isHelium = true; 782 isHelium = true;
783 tau = (0.511 / 3728.) * k; 783 tau = (0.511 / 3728.) * k;
784 } 784 }
785 785
786 G4double S = 4. * pi * Bohr_radius * Bohr_ra 786 G4double S = 4. * pi * Bohr_radius * Bohr_radius * n
787 * gpow->powN((Ry / Bj[ionizationLevelInd 787 * gpow->powN((Ry / Bj[ionizationLevelIndex]), 2);
788 if (j == 4) 788 if (j == 4)
789 S = 4. * pi * Bohr_radius * Bohr_radius * 789 S = 4. * pi * Bohr_radius * Bohr_radius * n
790 * gpow->powN((Ry / waterStructure.Ioni 790 * gpow->powN((Ry / waterStructure.IonisationEnergy(ionizationLevelIndex)),
791 2); 791 2);
792 792
793 G4double v2 = tau / Bj[ionizationLevelIndex] 793 G4double v2 = tau / Bj[ionizationLevelIndex];
794 if (j == 4) 794 if (j == 4)
795 v2 = tau / waterStructure.IonisationEnergy 795 v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
796 796
797 G4double v = std::sqrt(v2); 797 G4double v = std::sqrt(v2);
798 G4double wc = 4. * v2 - 2. * v - (Ry / (4. * 798 G4double wc = 4. * v2 - 2. * v - (Ry / (4. * Bj[ionizationLevelIndex]));
799 if (j == 4) 799 if (j == 4)
800 wc = 4. * v2 - 2. * v 800 wc = 4. * v2 - 2. * v
801 - (Ry / (4. * waterStructure.Ionisatio 801 - (Ry / (4. * waterStructure.IonisationEnergy(ionizationLevelIndex)));
802 802
803 G4double L1 = (C1 * gpow->powA(v, (D1))) / ( 803 G4double L1 = (C1 * gpow->powA(v, (D1))) / (1. + E1 * gpow->powA(v, (D1 + 4.)));
804 G4double L2 = C2 * gpow->powA(v, (D2)); 804 G4double L2 = C2 * gpow->powA(v, (D2));
805 G4double H1 = (A1 * G4Log(1. + v2)) / (v2 + 805 G4double H1 = (A1 * G4Log(1. + v2)) / (v2 + (B1 / v2));
806 G4double H2 = (A2 / v2) + (B2 / (v2 * v2)); 806 G4double H2 = (A2 / v2) + (B2 / (v2 * v2));
807 807
808 G4double F1 = L1 + H1; 808 G4double F1 = L1 + H1;
809 G4double F2 = (L2 * H2) / (L2 + H2); 809 G4double F2 = (L2 * H2) / (L2 + H2);
810 810
811 G4double sigma = 811 G4double sigma =
812 CorrectionFactor(particleDefinition, k) 812 CorrectionFactor(particleDefinition, k) * Gj[j]
813 * (S / Bj[ionizationLevelIndex]) 813 * (S / Bj[ionizationLevelIndex])
814 * ((F1 + w * F2) 814 * ((F1 + w * F2)
815 / (gpow->powN((1. + w), 3) 815 / (gpow->powN((1. + w), 3)
816 * (1. + G4Exp(alphaConst * ( 816 * (1. + G4Exp(alphaConst * (w - wc) / v))));
817 817
818 if (j == 4) 818 if (j == 4)
819 sigma = CorrectionFactor(particleDefinitio 819 sigma = CorrectionFactor(particleDefinition, k) * Gj[j]
820 * (S / waterStructure.IonisationEnergy 820 * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
821 * ((F1 + w * F2) 821 * ((F1 + w * F2)
822 / (gpow->powN((1. + w), 3) 822 / (gpow->powN((1. + w), 3)
823 * (1. + G4Exp(alphaConst * (w 823 * (1. + G4Exp(alphaConst * (w - wc) / v))));
824 824
825 if ((particleDefinition == hydrogenDef) 825 if ((particleDefinition == hydrogenDef)
826 && (ionizationLevelIndex == 4)) 826 && (ionizationLevelIndex == 4))
827 { 827 {
828 // sigma = Gj[j] * (S/Bj[ionizationLeve 828 // sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
829 sigma = Gj[j] * (S / waterStructure.Ionisa 829 sigma = Gj[j] * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
830 * ((F1 + w * F2) 830 * ((F1 + w * F2)
831 / (gpow->powN((1. + w), 3) 831 / (gpow->powN((1. + w), 3)
832 * (1. + G4Exp(alphaConst * (w 832 * (1. + G4Exp(alphaConst * (w - wc) / v))));
833 } 833 }
834 834
835 // if ( particleDefinition == protonDe 835 // if ( particleDefinition == protonDef
836 // || particleDefinition == hydro 836 // || particleDefinition == hydrogenDef
837 // ) 837 // )
838 838
839 if (isProtonOrHydrogen) 839 if (isProtonOrHydrogen)
840 { 840 {
841 return (sigma); 841 return (sigma);
842 } 842 }
843 843
844 if (particleDefinition == alphaPlusPlusDef) 844 if (particleDefinition == alphaPlusPlusDef)
845 { 845 {
846 slaterEffectiveCharge[0] = 0.; 846 slaterEffectiveCharge[0] = 0.;
847 slaterEffectiveCharge[1] = 0.; 847 slaterEffectiveCharge[1] = 0.;
848 slaterEffectiveCharge[2] = 0.; 848 slaterEffectiveCharge[2] = 0.;
849 sCoefficient[0] = 0.; 849 sCoefficient[0] = 0.;
850 sCoefficient[1] = 0.; 850 sCoefficient[1] = 0.;
851 sCoefficient[2] = 0.; 851 sCoefficient[2] = 0.;
852 } 852 }
853 853
854 else if (particleDefinition == alphaPlusDef) 854 else if (particleDefinition == alphaPlusDef)
855 { 855 {
856 slaterEffectiveCharge[0] = 2.0; 856 slaterEffectiveCharge[0] = 2.0;
857 // The following values are provided by M. 857 // The following values are provided by M. Dingfelder (priv. comm)
858 slaterEffectiveCharge[1] = 2.0; 858 slaterEffectiveCharge[1] = 2.0;
859 slaterEffectiveCharge[2] = 2.0; 859 slaterEffectiveCharge[2] = 2.0;
860 // 860 //
861 sCoefficient[0] = 0.7; 861 sCoefficient[0] = 0.7;
862 sCoefficient[1] = 0.15; 862 sCoefficient[1] = 0.15;
863 sCoefficient[2] = 0.15; 863 sCoefficient[2] = 0.15;
864 } 864 }
865 865
866 else if (particleDefinition == heliumDef) 866 else if (particleDefinition == heliumDef)
867 { 867 {
868 slaterEffectiveCharge[0] = 1.7; 868 slaterEffectiveCharge[0] = 1.7;
869 slaterEffectiveCharge[1] = 1.15; 869 slaterEffectiveCharge[1] = 1.15;
870 slaterEffectiveCharge[2] = 1.15; 870 slaterEffectiveCharge[2] = 1.15;
871 sCoefficient[0] = 0.5; 871 sCoefficient[0] = 0.5;
872 sCoefficient[1] = 0.25; 872 sCoefficient[1] = 0.25;
873 sCoefficient[2] = 0.25; 873 sCoefficient[2] = 0.25;
874 } 874 }
875 875
876 // if ( particleDefinition == heliumDe 876 // if ( particleDefinition == heliumDef
877 // || particleDefinition == alpha 877 // || particleDefinition == alphaPlusDef
878 // || particleDefinition == alpha 878 // || particleDefinition == alphaPlusPlusDef
879 // ) 879 // )
880 if (isHelium) 880 if (isHelium)
881 { 881 {
882 sigma = Gj[j] * (S / Bj[ionizationLevelInd 882 sigma = Gj[j] * (S / Bj[ionizationLevelIndex])
883 * ((F1 + w * F2) 883 * ((F1 + w * F2)
884 / (gpow->powN((1. + w), 3) 884 / (gpow->powN((1. + w), 3)
885 * (1. + G4Exp(alphaConst * (w 885 * (1. + G4Exp(alphaConst * (w - wc) / v))));
886 886
887 if (j == 4) 887 if (j == 4)
888 sigma = Gj[j] 888 sigma = Gj[j]
889 * (S / waterStructure.IonisationEner 889 * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
890 * ((F1 + w * F2) 890 * ((F1 + w * F2)
891 / (gpow->powN((1. + w), 3) 891 / (gpow->powN((1. + w), 3)
892 * (1. + G4Exp(alphaConst * ( 892 * (1. + G4Exp(alphaConst * (w - wc) / v))));
893 893
894 G4double zEff = particleDefinition->GetPDG 894 G4double zEff = particleDefinition->GetPDGCharge() / eplus
895 + particleDefinition->GetLeptonNumber( 895 + particleDefinition->GetLeptonNumber();
896 896
897 zEff -= (sCoefficient[0] 897 zEff -= (sCoefficient[0]
898 * S_1s(k, energyTransfer, slaterEffect 898 * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.)
899 + sCoefficient[1] 899 + sCoefficient[1]
900 * S_2s(k, energyTransfer, slaterEf 900 * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.)
901 + sCoefficient[2] 901 + sCoefficient[2]
902 * S_2p(k, energyTransfer, slaterEf 902 * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.));
903 903
904 return zEff * zEff * sigma; 904 return zEff * zEff * sigma;
905 } 905 }
906 906
907 return 0; 907 return 0;
908 } 908 }
909 909
910 //....oooOO0OOooo........oooOO0OOooo........oo 910 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
911 911
912 G4double G4DNARuddIonisationModel::S_1s(G4doub 912 G4double G4DNARuddIonisationModel::S_1s(G4double t,
913 G4doub 913 G4double energyTransferred,
914 G4doub 914 G4double slaterEffectiveChg,
915 G4doub 915 G4double shellNumber)
916 { 916 {
917 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2) 917 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
918 // Dingfelder, in Chattanooga 2005 proceedin 918 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
919 919
920 G4double r = R(t, energyTransferred, slaterE 920 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
921 G4double value = 1. - G4Exp(-2 * r) * ((2. * 921 G4double value = 1. - G4Exp(-2 * r) * ((2. * r + 2.) * r + 1.);
922 922
923 return value; 923 return value;
924 } 924 }
925 925
926 //....oooOO0OOooo........oooOO0OOooo........oo 926 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
927 927
928 G4double G4DNARuddIonisationModel::S_2s(G4doub 928 G4double G4DNARuddIonisationModel::S_2s(G4double t,
929 G4doub 929 G4double energyTransferred,
930 G4doub 930 G4double slaterEffectiveChg,
931 G4doub 931 G4double shellNumber)
932 { 932 {
933 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4) 933 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
934 // Dingfelder, in Chattanooga 2005 proceedin 934 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
935 935
936 G4double r = R(t, energyTransferred, slaterE 936 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
937 G4double value = 1. 937 G4double value = 1.
938 - G4Exp(-2 * r) * (((2. * r * r + 2.) * 938 - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
939 939
940 return value; 940 return value;
941 941
942 } 942 }
943 943
944 //....oooOO0OOooo........oooOO0OOooo........oo 944 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
945 945
946 G4double G4DNARuddIonisationModel::S_2p(G4doub 946 G4double G4DNARuddIonisationModel::S_2p(G4double t,
947 G4doub 947 G4double energyTransferred,
948 G4doub 948 G4double slaterEffectiveChg,
949 G4doub 949 G4double shellNumber)
950 { 950 {
951 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^ 951 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
952 // Dingfelder, in Chattanooga 2005 proceedin 952 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
953 953
954 G4double r = R(t, energyTransferred, slaterE 954 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
955 G4double value = 1. 955 G4double value = 1.
956 - G4Exp(-2 * r) 956 - G4Exp(-2 * r)
957 * ((((2. / 3. * r + 4. / 3.) * r + 2 957 * ((((2. / 3. * r + 4. / 3.) * r + 2.) * r + 2.) * r + 1.);
958 958
959 return value; 959 return value;
960 } 960 }
961 961
962 //....oooOO0OOooo........oooOO0OOooo........oo 962 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
963 963
964 G4double G4DNARuddIonisationModel::R(G4double 964 G4double G4DNARuddIonisationModel::R(G4double t,
965 G4double 965 G4double energyTransferred,
966 G4double 966 G4double slaterEffectiveChg,
967 G4double 967 G4double shellNumber)
968 { 968 {
969 // tElectron = m_electron / m_alpha * t 969 // tElectron = m_electron / m_alpha * t
970 // Dingfelder, in Chattanooga 2005 proceedin 970 // Dingfelder, in Chattanooga 2005 proceedings, p 4
971 971
972 G4double tElectron = 0.511 / 3728. * t; 972 G4double tElectron = 0.511 / 3728. * t;
973 // The following values are provided by M. D 973 // The following values are provided by M. Dingfelder (priv. comm)
974 G4double H = 2. * 13.60569172 * eV; 974 G4double H = 2. * 13.60569172 * eV;
975 G4double value = std::sqrt(2. * tElectron / 975 G4double value = std::sqrt(2. * tElectron / H) / (energyTransferred / H)
976 * (slaterEffectiveChg / shellNumber); 976 * (slaterEffectiveChg / shellNumber);
977 977
978 return value; 978 return value;
979 } 979 }
980 980
981 //....oooOO0OOooo........oooOO0OOooo........oo 981 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
982 982
983 G4double G4DNARuddIonisationModel::CorrectionF 983 G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition,
984 984 G4double k)
985 { 985 {
986 986
987 if (particleDefinition == G4Proton::Proton() 987 if (particleDefinition == G4Proton::Proton())
988 { 988 {
989 return (1.); 989 return (1.);
990 } 990 }
991 if (particleDefinition == hydrogenDef) 991 if (particleDefinition == hydrogenDef)
992 { 992 {
993 G4double value = (G4Log(k / eV)/gpow->logZ 993 G4double value = (G4Log(k / eV)/gpow->logZ(10) - 4.2) / 0.5;
994 // The following values are provided by M. 994 // The following values are provided by M. Dingfelder (priv. comm)
995 return ((0.6 / (1 + G4Exp(value))) + 0.9); 995 return ((0.6 / (1 + G4Exp(value))) + 0.9);
996 } 996 }
997 return (1.); 997 return (1.);
998 } 998 }
999 999
1000 //....oooOO0OOooo........oooOO0OOooo........o 1000 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1001 1001
1002 G4int G4DNARuddIonisationModel::RandomSelect( 1002 G4int G4DNARuddIonisationModel::RandomSelect(G4double k,
1003 1003 const G4String& particle)
1004 { 1004 {
1005 1005
1006 // BEGIN PART 1/2 OF ELECTRON CORRECTION 1006 // BEGIN PART 1/2 OF ELECTRON CORRECTION
1007 1007
1008 // add ONE or TWO electron-water ionisation 1008 // add ONE or TWO electron-water ionisation for alpha+ and helium
1009 1009
1010 G4int level = 0; 1010 G4int level = 0;
1011 1011
1012 // Retrieve data table corresponding to the 1012 // Retrieve data table corresponding to the current particle type
1013 1013
1014 std::map<G4String, G4DNACrossSectionDataSet 1014 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1015 pos = tableData.find(particle); 1015 pos = tableData.find(particle);
1016 1016
1017 if (pos != tableData.end()) 1017 if (pos != tableData.end())
1018 { 1018 {
1019 G4DNACrossSectionDataSet* table = pos->se 1019 G4DNACrossSectionDataSet* table = pos->second;
1020 1020
1021 if (table != nullptr) 1021 if (table != nullptr)
1022 { 1022 {
1023 auto valuesBuffer = new G4double[table 1023 auto valuesBuffer = new G4double[table->NumberOfComponents()];
1024 1024
1025 const auto n = (G4int)table->NumberOfC 1025 const auto n = (G4int)table->NumberOfComponents();
1026 G4int i(n); 1026 G4int i(n);
1027 G4double value = 0.; 1027 G4double value = 0.;
1028 1028
1029 while (i > 0) 1029 while (i > 0)
1030 { 1030 {
1031 --i; 1031 --i;
1032 valuesBuffer[i] = table->GetComponent 1032 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1033 value += valuesBuffer[i]; 1033 value += valuesBuffer[i];
1034 } 1034 }
1035 1035
1036 value *= G4UniformRand(); 1036 value *= G4UniformRand();
1037 1037
1038 i = n; 1038 i = n;
1039 1039
1040 while (i > 0) 1040 while (i > 0)
1041 { 1041 {
1042 --i; 1042 --i;
1043 1043
1044 if (valuesBuffer[i] > value) 1044 if (valuesBuffer[i] > value)
1045 { 1045 {
1046 delete[] valuesBuffer; 1046 delete[] valuesBuffer;
1047 return i; 1047 return i;
1048 } 1048 }
1049 value -= valuesBuffer[i]; 1049 value -= valuesBuffer[i];
1050 } 1050 }
1051 1051
1052 1052
1053 delete[] valuesBuffer; 1053 delete[] valuesBuffer;
1054 1054
1055 } 1055 }
1056 } else 1056 } else
1057 { 1057 {
1058 G4Exception("G4DNARuddIonisationModel::Ra 1058 G4Exception("G4DNARuddIonisationModel::RandomSelect",
1059 "em0002", 1059 "em0002",
1060 FatalException, 1060 FatalException,
1061 "Model not applicable to part 1061 "Model not applicable to particle type.");
1062 } 1062 }
1063 1063
1064 return level; 1064 return level;
1065 } 1065 }
1066 1066
1067 //....oooOO0OOooo........oooOO0OOooo........o 1067 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1068 1068
1069 G4double G4DNARuddIonisationModel::PartialCro 1069 G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track)
1070 { 1070 {
1071 G4double sigma = 0.; 1071 G4double sigma = 0.;
1072 1072
1073 const G4DynamicParticle* particle = track.G 1073 const G4DynamicParticle* particle = track.GetDynamicParticle();
1074 G4double k = particle->GetKineticEnergy(); 1074 G4double k = particle->GetKineticEnergy();
1075 1075
1076 G4double lowLim = 0; 1076 G4double lowLim = 0;
1077 G4double highLim = 0; 1077 G4double highLim = 0;
1078 1078
1079 const G4String& particleName = particle->Ge 1079 const G4String& particleName = particle->GetDefinition()->GetParticleName();
1080 1080
1081 std::map<G4String, G4double, std::less<G4St 1081 std::map<G4String, G4double, std::less<G4String> >::iterator pos1;
1082 pos1 = lowEnergyLimit.find(particleName); 1082 pos1 = lowEnergyLimit.find(particleName);
1083 1083
1084 if (pos1 != lowEnergyLimit.end()) 1084 if (pos1 != lowEnergyLimit.end())
1085 { 1085 {
1086 lowLim = pos1->second; 1086 lowLim = pos1->second;
1087 } 1087 }
1088 1088
1089 std::map<G4String, G4double, std::less<G4St 1089 std::map<G4String, G4double, std::less<G4String> >::iterator pos2;
1090 pos2 = highEnergyLimit.find(particleName); 1090 pos2 = highEnergyLimit.find(particleName);
1091 1091
1092 if (pos2 != highEnergyLimit.end()) 1092 if (pos2 != highEnergyLimit.end())
1093 { 1093 {
1094 highLim = pos2->second; 1094 highLim = pos2->second;
1095 } 1095 }
1096 1096
1097 if (k >= lowLim && k <= highLim) 1097 if (k >= lowLim && k <= highLim)
1098 { 1098 {
1099 std::map<G4String, G4DNACrossSectionDataS 1099 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1100 pos = tableData.find(particleName); 1100 pos = tableData.find(particleName);
1101 1101
1102 if (pos != tableData.end()) 1102 if (pos != tableData.end())
1103 { 1103 {
1104 G4DNACrossSectionDataSet* table = pos-> 1104 G4DNACrossSectionDataSet* table = pos->second;
1105 if (table != nullptr) 1105 if (table != nullptr)
1106 { 1106 {
1107 sigma = table->FindValue(k); 1107 sigma = table->FindValue(k);
1108 } 1108 }
1109 } else 1109 } else
1110 { 1110 {
1111 G4Exception("G4DNARuddIonisationModel:: 1111 G4Exception("G4DNARuddIonisationModel::PartialCrossSection",
1112 "em0002", 1112 "em0002",
1113 FatalException, 1113 FatalException,
1114 "Model not applicable to pa 1114 "Model not applicable to particle type.");
1115 } 1115 }
1116 } 1116 }
1117 1117
1118 return sigma; 1118 return sigma;
1119 } 1119 }
1120 1120
1121 //....oooOO0OOooo........oooOO0OOooo........o 1121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1122 1122
1123 G4double G4DNARuddIonisationModel::Sum(G4doub 1123 G4double G4DNARuddIonisationModel::Sum(G4double /* energy */,
1124 const 1124 const G4String& /* particle */)
1125 { 1125 {
1126 return 0; 1126 return 0;
1127 } 1127 }
1128 1128
1129 1129