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