| 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 "G4DNADingfelderChargeDecreaseModel.h 28 #include "G4DNADingfelderChargeDecreaseModel.hh"
29 #include "G4PhysicalConstants.hh" 29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh"
31 #include "G4DNAMolecularMaterial.hh" 31 #include "G4DNAMolecularMaterial.hh"
32 #include "G4DNAChemistryManager.hh" 32 #include "G4DNAChemistryManager.hh"
33 #include "G4Log.hh" <<
34 #include "G4Pow.hh" <<
35 #include "G4Alpha.hh" <<
36 <<
37 <<
38 static G4Pow * gpow = G4Pow::GetInstance(); <<
39 33
40 //....oooOO0OOooo........oooOO0OOooo........oo 34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 35
42 using namespace std; 36 using namespace std;
43 37
44 //....oooOO0OOooo........oooOO0OOooo........oo 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 39
46 G4DNADingfelderChargeDecreaseModel::G4DNADingf 40 G4DNADingfelderChargeDecreaseModel::G4DNADingfelderChargeDecreaseModel(const G4ParticleDefinition*,
47 41 const G4String& nam) :
48 G4VEmModel(nam) << 42 G4VEmModel(nam), isInitialised(false)
49 { 43 {
>> 44 fpMolWaterDensity = 0;
50 numberOfPartialCrossSections[0] = 0; 45 numberOfPartialCrossSections[0] = 0;
51 numberOfPartialCrossSections[1] = 0; 46 numberOfPartialCrossSections[1] = 0;
52 numberOfPartialCrossSections[2] = 0; 47 numberOfPartialCrossSections[2] = 0;
53 48
54 verboseLevel = 0; 49 verboseLevel = 0;
55 // Verbosity scale: 50 // Verbosity scale:
56 // 0 = nothing 51 // 0 = nothing
57 // 1 = warning for energy non-conservation 52 // 1 = warning for energy non-conservation
58 // 2 = details of energy budget 53 // 2 = details of energy budget
59 // 3 = calculation of cross sections, file o 54 // 3 = calculation of cross sections, file openings, sampling of atoms
60 // 4 = entering in methods 55 // 4 = entering in methods
61 56
62 if (verboseLevel > 0) 57 if (verboseLevel > 0)
63 { 58 {
64 G4cout << "Dingfelder charge decrease mode 59 G4cout << "Dingfelder charge decrease model is constructed " << G4endl;
65 } 60 }
>> 61 fParticleChangeForGamma = 0;
>> 62
66 // Selection of stationary mode 63 // Selection of stationary mode
67 64
68 statCode = false; 65 statCode = false;
69 } 66 }
70 67
71 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 69
>> 70 G4DNADingfelderChargeDecreaseModel::~G4DNADingfelderChargeDecreaseModel()
>> 71 {
>> 72 }
>> 73
>> 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 75
73 void G4DNADingfelderChargeDecreaseModel::Initi 76 void G4DNADingfelderChargeDecreaseModel::Initialise(const G4ParticleDefinition* particle,
74 77 const G4DataVector& /*cuts*/)
75 { 78 {
76 79
77 if (verboseLevel > 3) 80 if (verboseLevel > 3)
78 { 81 {
79 G4cout << "Calling G4DNADingfelderChargeDe 82 G4cout << "Calling G4DNADingfelderChargeDecreaseModel::Initialise()"
80 << G4endl; 83 << G4endl;
81 } 84 }
82 85
83 // Energy limits 86 // Energy limits
84 87
85 G4DNAGenericIonsManager *instance; 88 G4DNAGenericIonsManager *instance;
86 instance = G4DNAGenericIonsManager::Instance 89 instance = G4DNAGenericIonsManager::Instance();
87 protonDef = G4Proton::ProtonDefinition(); << 90 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
88 alphaPlusPlusDef = G4Alpha::Alpha(); << 91 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
89 alphaPlusDef = instance->GetIon("alpha+"); << 92 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
90 hydrogenDef = instance->GetIon("hydrogen"); <<
91 heliumDef = instance->GetIon("helium"); <<
92 93
93 G4String proton; 94 G4String proton;
94 G4String alphaPlusPlus; 95 G4String alphaPlusPlus;
95 G4String alphaPlus; 96 G4String alphaPlus;
96 97
97 // Limits 98 // Limits
98 99
99 proton = protonDef->GetParticleName(); 100 proton = protonDef->GetParticleName();
100 lowEnergyLimit[proton] = 100. * eV; 101 lowEnergyLimit[proton] = 100. * eV;
101 highEnergyLimit[proton] = 100. * MeV; 102 highEnergyLimit[proton] = 100. * MeV;
102 103
103 alphaPlusPlus = alphaPlusPlusDef->GetParticl 104 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
104 lowEnergyLimit[alphaPlusPlus] = 1. * keV; 105 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
105 highEnergyLimit[alphaPlusPlus] = 400. * MeV; 106 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
106 107
107 alphaPlus = alphaPlusDef->GetParticleName(); 108 alphaPlus = alphaPlusDef->GetParticleName();
108 lowEnergyLimit[alphaPlus] = 1. * keV; 109 lowEnergyLimit[alphaPlus] = 1. * keV;
109 highEnergyLimit[alphaPlus] = 400. * MeV; 110 highEnergyLimit[alphaPlus] = 400. * MeV;
110 111
111 // 112 //
112 113
113 if (particle==protonDef) 114 if (particle==protonDef)
114 { 115 {
115 SetLowEnergyLimit(lowEnergyLimit[proton]); 116 SetLowEnergyLimit(lowEnergyLimit[proton]);
116 SetHighEnergyLimit(highEnergyLimit[proton] 117 SetHighEnergyLimit(highEnergyLimit[proton]);
117 } 118 }
118 119
119 if (particle==alphaPlusPlusDef) 120 if (particle==alphaPlusPlusDef)
120 { 121 {
121 SetLowEnergyLimit(lowEnergyLimit[alphaPlus 122 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
122 SetHighEnergyLimit(highEnergyLimit[alphaPl 123 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
123 } 124 }
124 125
125 if (particle==alphaPlusDef) 126 if (particle==alphaPlusDef)
126 { 127 {
127 SetLowEnergyLimit(lowEnergyLimit[alphaPlus 128 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
128 SetHighEnergyLimit(highEnergyLimit[alphaPl 129 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
129 } 130 }
130 131
131 // Final state 132 // Final state
132 133
133 //PROTON 134 //PROTON
134 f0[0][0]=1.; 135 f0[0][0]=1.;
135 a0[0][0]=-0.180; 136 a0[0][0]=-0.180;
136 a1[0][0]=-3.600; 137 a1[0][0]=-3.600;
137 b0[0][0]=-18.22; 138 b0[0][0]=-18.22;
138 b1[0][0]=-1.997; 139 b1[0][0]=-1.997;
139 c0[0][0]=0.215; 140 c0[0][0]=0.215;
140 d0[0][0]=3.550; 141 d0[0][0]=3.550;
141 x0[0][0]=3.450; 142 x0[0][0]=3.450;
142 x1[0][0]=5.251; 143 x1[0][0]=5.251;
143 144
144 numberOfPartialCrossSections[0] = 1; 145 numberOfPartialCrossSections[0] = 1;
145 146
146 //ALPHA++ 147 //ALPHA++
147 f0[0][1]=1.; a0[0][1]=0.95; 148 f0[0][1]=1.; a0[0][1]=0.95;
148 a1[0][1]=-2.75; 149 a1[0][1]=-2.75;
149 b0[0][1]=-23.00; 150 b0[0][1]=-23.00;
150 c0[0][1]=0.215; 151 c0[0][1]=0.215;
151 d0[0][1]=2.95; 152 d0[0][1]=2.95;
152 x0[0][1]=3.50; 153 x0[0][1]=3.50;
153 154
154 f0[1][1]=1.; 155 f0[1][1]=1.;
155 a0[1][1]=0.95; 156 a0[1][1]=0.95;
156 a1[1][1]=-2.75; 157 a1[1][1]=-2.75;
157 b0[1][1]=-23.73; 158 b0[1][1]=-23.73;
158 c0[1][1]=0.250; 159 c0[1][1]=0.250;
159 d0[1][1]=3.55; 160 d0[1][1]=3.55;
160 x0[1][1]=3.72; 161 x0[1][1]=3.72;
161 162
162 x1[0][1]=-1.; 163 x1[0][1]=-1.;
163 b1[0][1]=-1.; 164 b1[0][1]=-1.;
164 165
165 x1[1][1]=-1.; 166 x1[1][1]=-1.;
166 b1[1][1]=-1.; 167 b1[1][1]=-1.;
167 168
168 numberOfPartialCrossSections[1] = 2; 169 numberOfPartialCrossSections[1] = 2;
169 170
170 // ALPHA+ 171 // ALPHA+
171 f0[0][2]=1.; 172 f0[0][2]=1.;
172 a0[0][2]=0.65; 173 a0[0][2]=0.65;
173 a1[0][2]=-2.75; 174 a1[0][2]=-2.75;
174 b0[0][2]=-21.81; 175 b0[0][2]=-21.81;
175 c0[0][2]=0.232; 176 c0[0][2]=0.232;
176 d0[0][2]=2.95; 177 d0[0][2]=2.95;
177 x0[0][2]=3.53; 178 x0[0][2]=3.53;
178 179
179 x1[0][2]=-1.; 180 x1[0][2]=-1.;
180 b1[0][2]=-1.; 181 b1[0][2]=-1.;
181 182
182 numberOfPartialCrossSections[2] = 1; 183 numberOfPartialCrossSections[2] = 1;
183 184
184 // 185 //
185 186
186 if( verboseLevel>0 ) 187 if( verboseLevel>0 )
187 { 188 {
188 G4cout << "Dingfelder charge decrease mode 189 G4cout << "Dingfelder charge decrease model is initialized " << G4endl
189 << "Energy range: " 190 << "Energy range: "
190 << LowEnergyLimit() / keV << " keV - " 191 << LowEnergyLimit() / keV << " keV - "
191 << HighEnergyLimit() / MeV << " MeV for " 192 << HighEnergyLimit() / MeV << " MeV for "
192 << particle->GetParticleName() 193 << particle->GetParticleName()
193 << G4endl; 194 << G4endl;
194 } 195 }
195 196
196 // Initialize water density pointer 197 // Initialize water density pointer
197 fpMolWaterDensity = G4DNAMolecularMaterial:: 198 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
198 199
199 if (isInitialised) 200 if (isInitialised)
200 { return;} 201 { return;}
201 fParticleChangeForGamma = GetParticleChangeF 202 fParticleChangeForGamma = GetParticleChangeForGamma();
202 isInitialised = true; 203 isInitialised = true;
203 204
204 } 205 }
205 206
206 //....oooOO0OOooo........oooOO0OOooo........oo 207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
207 208
208 G4double G4DNADingfelderChargeDecreaseModel::C 209 G4double G4DNADingfelderChargeDecreaseModel::CrossSectionPerVolume(const G4Material* material,
209 210 const G4ParticleDefinition* particleDefinition,
210 211 G4double k,
211 212 G4double,
212 213 G4double)
213 { 214 {
214 if (verboseLevel > 3) 215 if (verboseLevel > 3)
215 { 216 {
216 G4cout 217 G4cout
217 << "Calling CrossSectionPerVolume() of 218 << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeDecreaseModel"
218 << G4endl; 219 << G4endl;
219 } 220 }
220 221
221 // Calculate total cross section for model 222 // Calculate total cross section for model
>> 223
>> 224 G4DNAGenericIonsManager *instance;
>> 225 instance = G4DNAGenericIonsManager::Instance();
>> 226
222 if ( 227 if (
223 particleDefinition != protonDef << 228 particleDefinition != G4Proton::ProtonDefinition()
224 && 229 &&
225 particleDefinition != alphaPlusPlusDef << 230 particleDefinition != instance->GetIon("alpha++")
226 && 231 &&
227 particleDefinition != alphaPlusDef << 232 particleDefinition != instance->GetIon("alpha+")
228 ) 233 )
229 234
230 return 0; 235 return 0;
231 236
232 G4double lowLim = 0; 237 G4double lowLim = 0;
233 G4double highLim = 0; 238 G4double highLim = 0;
234 G4double crossSection = 0.; 239 G4double crossSection = 0.;
235 240
236 G4double waterDensity = (*fpMolWaterDensity) 241 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
237 242
238 const G4String& particleName = particleDefin 243 const G4String& particleName = particleDefinition->GetParticleName();
239 244
240 std::map< G4String,G4double,std::less<G4Stri 245 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
241 pos1 = lowEnergyLimit.find(particleName); 246 pos1 = lowEnergyLimit.find(particleName);
242 247
243 if (pos1 != lowEnergyLimit.end()) 248 if (pos1 != lowEnergyLimit.end())
244 { 249 {
245 lowLim = pos1->second; 250 lowLim = pos1->second;
246 } 251 }
247 252
248 std::map< G4String,G4double,std::less<G4Stri 253 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
249 pos2 = highEnergyLimit.find(particleName); 254 pos2 = highEnergyLimit.find(particleName);
250 255
251 if (pos2 != highEnergyLimit.end()) 256 if (pos2 != highEnergyLimit.end())
252 { 257 {
253 highLim = pos2->second; 258 highLim = pos2->second;
254 } 259 }
255 260
256 if (k >= lowLim && k <= highLim) 261 if (k >= lowLim && k <= highLim)
257 { 262 {
258 crossSection = Sum(k,particleDefinition); 263 crossSection = Sum(k,particleDefinition);
259 } 264 }
260 265
261 if (verboseLevel > 2) 266 if (verboseLevel > 2)
262 { 267 {
263 G4cout << "_______________________________ 268 G4cout << "_______________________________________" << G4endl;
264 G4cout << "G4DNADingfelderChargeDecreaeMod 269 G4cout << "G4DNADingfelderChargeDecreaeModel" << G4endl;
265 G4cout << "Kinetic energy(eV)=" << k/eV << 270 G4cout << "Kinetic energy(eV)=" << k/eV << "particle :" << particleName << G4endl;
266 G4cout << "Cross section per water molecul 271 G4cout << "Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
267 G4cout << "Cross section per water molecul 272 G4cout << "Cross section per water molecule (cm^-1)=" << crossSection*
268 waterDensity/(1./cm) << G4endl; 273 waterDensity/(1./cm) << G4endl;
269 // material->GetAtomicNumDensityVector()[1 274 // material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
270 } 275 }
271 276
272 return crossSection*waterDensity; 277 return crossSection*waterDensity;
273 //return crossSection*material->GetAtomicNum 278 //return crossSection*material->GetAtomicNumDensityVector()[1];
274 279
275 } 280 }
276 281
277 //....oooOO0OOooo........oooOO0OOooo........oo 282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
278 283
279 void G4DNADingfelderChargeDecreaseModel::Sampl 284 void G4DNADingfelderChargeDecreaseModel::SampleSecondaries(std::vector<
280 285 G4DynamicParticle*>* fvect,
281 286 const G4MaterialCutsCouple* /*couple*/,
282 287 const G4DynamicParticle* aDynamicParticle,
283 288 G4double,
284 289 G4double)
285 { 290 {
286 if (verboseLevel > 3) 291 if (verboseLevel > 3)
287 { 292 {
288 G4cout 293 G4cout
289 << "Calling SampleSecondaries() of G4D 294 << "Calling SampleSecondaries() of G4DNADingfelderChargeDecreaseModel"
290 << G4endl; 295 << G4endl;
291 } 296 }
292 297
293 G4double inK = aDynamicParticle->GetKineticE 298 G4double inK = aDynamicParticle->GetKineticEnergy();
294 299
295 G4ParticleDefinition* definition = aDynamicP 300 G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
296 301
297 G4double particleMass = definition->GetPDGMa 302 G4double particleMass = definition->GetPDGMass();
298 303
299 G4int finalStateIndex = RandomSelect(inK,def 304 G4int finalStateIndex = RandomSelect(inK,definition);
300 305
301 G4int n = NumberOfFinalStates(definition, fi 306 G4int n = NumberOfFinalStates(definition, finalStateIndex);
302 G4double waterBindingEnergy = WaterBindingEn 307 G4double waterBindingEnergy = WaterBindingEnergyConstant(definition, finalStateIndex);
303 G4double outgoingParticleBindingEnergy = Out 308 G4double outgoingParticleBindingEnergy = OutgoingParticleBindingEnergyConstant(definition, finalStateIndex);
304 309
305 G4double outK = 0.; 310 G4double outK = 0.;
306 311
307 if (definition==G4Proton::Proton()) 312 if (definition==G4Proton::Proton())
308 { 313 {
309 if (!statCode) outK = inK - n*(inK*electron 314 if (!statCode) outK = inK - n*(inK*electron_mass_c2/proton_mass_c2)
310 - waterBindingEne 315 - waterBindingEnergy + outgoingParticleBindingEnergy;
311 else outK = inK; 316 else outK = inK;
312 } 317 }
313 318
314 else 319 else
315 { 320 {
316 if (!statCode) outK = inK - n*(inK*electron 321 if (!statCode) outK = inK - n*(inK*electron_mass_c2/particleMass)
317 - waterBindingEne 322 - waterBindingEnergy + outgoingParticleBindingEnergy;
318 else outK = inK; 323 else outK = inK;
319 } 324 }
320 325
321 if (outK<0) 326 if (outK<0)
322 { 327 {
323 G4Exception("G4DNADingfelderChargeDecrease 328 G4Exception("G4DNADingfelderChargeDecreaseModel::SampleSecondaries","em0004",
324 FatalException,"Final kinetic energy i 329 FatalException,"Final kinetic energy is negative.");
325 } 330 }
326 331
327 fParticleChangeForGamma->ProposeTrackStatus( 332 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
328 333
329 if (!statCode) fParticleChangeForGamma->Prop 334 if (!statCode) fParticleChangeForGamma->ProposeLocalEnergyDeposit(waterBindingEnergy);
330 335
331 else 336 else
332 337
333 { 338 {
334 if (definition==G4Proton::Proton()) 339 if (definition==G4Proton::Proton())
335 fParticleChangeForGamma->ProposeLocalEner 340 fParticleChangeForGamma->ProposeLocalEnergyDeposit(n*(inK*electron_mass_c2/proton_mass_c2)
336 + waterBindingEnergy - outgoingParticleBi 341 + waterBindingEnergy - outgoingParticleBindingEnergy);
337 else 342 else
338 fParticleChangeForGamma->ProposeLocalEner 343 fParticleChangeForGamma->ProposeLocalEnergyDeposit(n*(inK*electron_mass_c2/particleMass)
339 + waterBindingEnergy - outgoingParticleBi 344 + waterBindingEnergy - outgoingParticleBindingEnergy);
340 } 345 }
341 346
342 auto dp = new G4DynamicParticle (OutgoingPa << 347 G4DynamicParticle* dp = new G4DynamicParticle (OutgoingParticleDefinition(definition, finalStateIndex),
343 aDynamicParticle->GetMomentumDirection() 348 aDynamicParticle->GetMomentumDirection(),
344 outK); 349 outK);
345 fvect->push_back(dp); 350 fvect->push_back(dp);
346 351
347 const G4Track * theIncomingTrack = fParticle 352 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
348 G4DNAChemistryManager::Instance()->Creat 353 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
349 1, 354 1,
350 theIncomingTrack); 355 theIncomingTrack);
351 356
352 } 357 }
353 358
354 //....oooOO0OOooo........oooOO0OOooo........oo 359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
355 360
356 G4int G4DNADingfelderChargeDecreaseModel::Numb 361 G4int G4DNADingfelderChargeDecreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
357 362 G4int finalStateIndex)
358 363
359 { 364 {
360 if (particleDefinition == G4Proton::Proton() 365 if (particleDefinition == G4Proton::Proton())
361 return 1; 366 return 1;
362 367
363 if (particleDefinition == alphaPlusPlusDef) << 368 G4DNAGenericIonsManager*instance;
>> 369 instance = G4DNAGenericIonsManager::Instance();
>> 370
>> 371 if (particleDefinition == instance->GetIon("alpha++"))
364 { 372 {
365 if (finalStateIndex == 0) 373 if (finalStateIndex == 0)
366 return 1; 374 return 1;
367 return 2; 375 return 2;
368 } 376 }
369 377
370 if (particleDefinition == alphaPlusDef) << 378 if (particleDefinition == instance->GetIon("alpha+"))
371 return 1; 379 return 1;
372 380
373 return 0; 381 return 0;
374 } 382 }
375 383
376 //....oooOO0OOooo........oooOO0OOooo........oo 384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377 385
378 G4ParticleDefinition* G4DNADingfelderChargeDec 386 G4ParticleDefinition* G4DNADingfelderChargeDecreaseModel::OutgoingParticleDefinition(G4ParticleDefinition* particleDefinition,
379 387 G4int finalStateIndex)
380 { 388 {
>> 389 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
>> 390
381 if (particleDefinition == G4Proton::Proton() 391 if (particleDefinition == G4Proton::Proton())
382 return hydrogenDef; << 392 return instance->GetIon("hydrogen");
383 393
384 if (particleDefinition == alphaPlusPlusDef) << 394 if (particleDefinition == instance->GetIon("alpha++"))
385 { 395 {
386 if (finalStateIndex == 0) 396 if (finalStateIndex == 0)
387 return alphaPlusDef; << 397 return instance->GetIon("alpha+");
388 return heliumDef; << 398 return instance->GetIon("helium");
389 } 399 }
390 400
391 if (particleDefinition == alphaPlusDef) << 401 if (particleDefinition == instance->GetIon("alpha+"))
392 return heliumDef; << 402 return instance->GetIon("helium");
393 403
394 return nullptr; << 404 return 0;
395 } 405 }
396 406
397 //....oooOO0OOooo........oooOO0OOooo........oo 407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
398 408
399 G4double G4DNADingfelderChargeDecreaseModel::W 409 G4double G4DNADingfelderChargeDecreaseModel::WaterBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
400 410 G4int finalStateIndex)
401 { 411 {
402 // Ionization energy of first water shell 412 // Ionization energy of first water shell
403 // Rad. Phys. Chem. 59 p.267 by Dingf. et al 413 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
404 // W + 10.79 eV -> W+ + e- 414 // W + 10.79 eV -> W+ + e-
405 415
>> 416 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
>> 417
406 if (particleDefinition == G4Proton::Proton() 418 if (particleDefinition == G4Proton::Proton())
407 return 10.79 * eV; 419 return 10.79 * eV;
408 420
409 if (particleDefinition == alphaPlusPlusDef) << 421 if (particleDefinition == instance->GetIon("alpha++"))
410 { 422 {
411 // Binding energy for W+ -> W++ + e- 423 // Binding energy for W+ -> W++ + e- 10.79 eV
412 // Binding energy for W -> W+ + e- 424 // Binding energy for W -> W+ + e- 10.79 eV
413 // 425 //
414 // Ionization energy of first water shell 426 // Ionization energy of first water shell
415 // Rad. Phys. Chem. 59 p.267 by Dingf. et 427 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
416 428
417 if (finalStateIndex == 0) 429 if (finalStateIndex == 0)
418 return 10.79 * eV; 430 return 10.79 * eV;
419 431
420 return 10.79 * 2 * eV; 432 return 10.79 * 2 * eV;
421 } 433 }
422 434
423 if (particleDefinition == alphaPlusDef) << 435 if (particleDefinition == instance->GetIon("alpha+"))
424 { 436 {
425 // Binding energy for W+ -> W++ + e- 437 // Binding energy for W+ -> W++ + e- 10.79 eV
426 // Binding energy for W -> W+ + e- 438 // Binding energy for W -> W+ + e- 10.79 eV
427 // 439 //
428 // Ionization energy of first water shell 440 // Ionization energy of first water shell
429 // Rad. Phys. Chem. 59 p.267 by Dingf. et 441 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
430 442
431 return 10.79 * eV; 443 return 10.79 * eV;
432 } 444 }
433 445
434 return 0; 446 return 0;
435 } 447 }
436 448
437 //....oooOO0OOooo........oooOO0OOooo........oo 449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
438 450
439 G4double G4DNADingfelderChargeDecreaseModel::O 451 G4double G4DNADingfelderChargeDecreaseModel::OutgoingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
440 452 G4int finalStateIndex)
441 { 453 {
>> 454 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
>> 455
442 if (particleDefinition == G4Proton::Proton() 456 if (particleDefinition == G4Proton::Proton())
443 return 13.6 * eV; 457 return 13.6 * eV;
444 458
445 if (particleDefinition == alphaPlusPlusDef) << 459 if (particleDefinition == instance->GetIon("alpha++"))
446 { 460 {
447 // Binding energy for He+ -> He++ + e- 461 // Binding energy for He+ -> He++ + e- 54.509 eV
448 // Binding energy for He -> He+ + e- 462 // Binding energy for He -> He+ + e- 24.587 eV
449 463
450 if (finalStateIndex == 0) 464 if (finalStateIndex == 0)
451 return 54.509 * eV; 465 return 54.509 * eV;
452 466
453 return (54.509 + 24.587) * eV; 467 return (54.509 + 24.587) * eV;
454 } 468 }
455 469
456 if (particleDefinition == alphaPlusDef) << 470 if (particleDefinition == instance->GetIon("alpha+"))
457 { 471 {
458 // Binding energy for He+ -> He++ + e- 472 // Binding energy for He+ -> He++ + e- 54.509 eV
459 // Binding energy for He -> He+ + e- 473 // Binding energy for He -> He+ + e- 24.587 eV
460 474
461 return 24.587 * eV; 475 return 24.587 * eV;
462 } 476 }
463 477
464 return 0; 478 return 0;
465 } 479 }
466 480
467 //....oooOO0OOooo........oooOO0OOooo........oo 481 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
468 482
469 G4double G4DNADingfelderChargeDecreaseModel::P 483 G4double G4DNADingfelderChargeDecreaseModel::PartialCrossSection(G4double k,
470 484 G4int index,
471 485 const G4ParticleDefinition* particleDefinition)
472 { 486 {
473 G4int particleTypeIndex = 0; 487 G4int particleTypeIndex = 0;
>> 488 G4DNAGenericIonsManager* instance;
>> 489 instance = G4DNAGenericIonsManager::Instance();
474 490
475 if (particleDefinition == protonDef) << 491 if (particleDefinition == G4Proton::ProtonDefinition())
476 particleTypeIndex = 0; 492 particleTypeIndex = 0;
477 493
478 if (particleDefinition == alphaPlusPlusDef) << 494 if (particleDefinition == instance->GetIon("alpha++"))
479 particleTypeIndex = 1; 495 particleTypeIndex = 1;
480 496
481 if (particleDefinition == alphaPlusDef) << 497 if (particleDefinition == instance->GetIon("alpha+"))
482 particleTypeIndex = 2; 498 particleTypeIndex = 2;
483 499
484 // 500 //
485 // sigma(T) = f0 10 ^ y(log10(T/eV)) 501 // sigma(T) = f0 10 ^ y(log10(T/eV))
486 // 502 //
487 // / a0 x + b0 x 503 // / a0 x + b0 x < x0
488 // | 504 // |
489 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x 505 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1
490 // | 506 // |
491 // \ a1 x + b1 x 507 // \ a1 x + b1 x >= x1
492 // 508 //
493 // 509 //
494 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are pa 510 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
495 // 511 //
496 // f0 has been added to the code in order to 512 // f0 has been added to the code in order to manage partial (shell-dependent) cross sections (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1)
497 // 513 //
498 // From Rad. Phys. and Chem. 59 (2000) 255-2 514 // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
499 // Inelastic-collision cross sections of liq 515 // Inelastic-collision cross sections of liquid water for interactions of energetic proton
500 // 516 //
501 517
502 if (x1[index][particleTypeIndex] < x0[index] 518 if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex])
503 { 519 {
504 // 520 //
505 // if x1 < x0 means that x1 and b1 will be 521 // if x1 < x0 means that x1 and b1 will be calculated with the following formula (this piece of code is run on all alphas and not on protons)
506 // 522 //
507 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / 523 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
508 // 524 //
509 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x 525 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
510 // 526 //
511 527
512 x1[index][particleTypeIndex] = x0[index][p 528 x1[index][particleTypeIndex] = x0[index][particleTypeIndex]
513 + gpow->powA((a0[index][particleTypeIn << 529 + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
514 / (c0[index][particleTy 530 / (c0[index][particleTypeIndex]
515 * d0[index][particl 531 * d0[index][particleTypeIndex]),
516 1. / (d0[index][particleTyp 532 1. / (d0[index][particleTypeIndex] - 1.));
517 b1[index][particleTypeIndex] = (a0[index][ 533 b1[index][particleTypeIndex] = (a0[index][particleTypeIndex]
518 - a1[index][particleTypeIndex]) * x1[i 534 - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
519 + b0[index][particleTypeIndex] 535 + b0[index][particleTypeIndex]
520 - c0[index][particleTypeIndex] 536 - c0[index][particleTypeIndex]
521 * gpow->powA(x1[index][particleTyp << 537 * std::pow(x1[index][particleTypeIndex]
522 - x0[index][particl 538 - x0[index][particleTypeIndex],
523 d0[index][particleTypeI 539 d0[index][particleTypeIndex]);
524 } 540 }
525 541
526 G4double x(G4Log(k / eV)/gpow->logZ(10)); << 542 G4double x(std::log10(k / eV));
527 G4double y; 543 G4double y;
528 544
529 if (x < x0[index][particleTypeIndex]) 545 if (x < x0[index][particleTypeIndex])
530 y = a0[index][particleTypeIndex] * x + b0[ 546 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
531 else if (x < x1[index][particleTypeIndex]) 547 else if (x < x1[index][particleTypeIndex])
532 y = a0[index][particleTypeIndex] * x + b0[ 548 y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]
533 - c0[index][particleTypeIndex] 549 - c0[index][particleTypeIndex]
534 * gpow->powA(x - x0[index][particl << 550 * std::pow(x - x0[index][particleTypeIndex],
535 d0[index][particleTypeI 551 d0[index][particleTypeIndex]);
536 else 552 else
537 y = a1[index][particleTypeIndex] * x + b1[ 553 y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
538 554
539 return f0[index][particleTypeIndex] * gpow-> << 555 return f0[index][particleTypeIndex] * std::pow(10., y) * m * m;
540 556
541 } 557 }
542 558
543 G4int G4DNADingfelderChargeDecreaseModel::Rand 559 G4int G4DNADingfelderChargeDecreaseModel::RandomSelect(G4double k,
544 560 const G4ParticleDefinition* particleDefinition)
545 { 561 {
546 G4int particleTypeIndex = 0; 562 G4int particleTypeIndex = 0;
>> 563 G4DNAGenericIonsManager *instance;
>> 564 instance = G4DNAGenericIonsManager::Instance();
547 565
548 if (particleDefinition == protonDef) << 566 if (particleDefinition == G4Proton::ProtonDefinition())
549 particleTypeIndex = 0; 567 particleTypeIndex = 0;
550 568
551 if (particleDefinition == alphaPlusPlusDef) << 569 if (particleDefinition == instance->GetIon("alpha++"))
552 particleTypeIndex = 1; 570 particleTypeIndex = 1;
553 571
554 if (particleDefinition == alphaPlusDef) << 572 if (particleDefinition == instance->GetIon("alpha+"))
555 particleTypeIndex = 2; 573 particleTypeIndex = 2;
556 574
557 const G4int n = numberOfPartialCrossSections 575 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
558 auto values(new G4double[n]); << 576 G4double* values(new G4double[n]);
559 G4double value(0); 577 G4double value(0);
560 G4int i = n; 578 G4int i = n;
561 579
562 while (i > 0) 580 while (i > 0)
563 { 581 {
564 i--; 582 i--;
565 values[i] = PartialCrossSection(k, i, part 583 values[i] = PartialCrossSection(k, i, particleDefinition);
566 value += values[i]; 584 value += values[i];
567 } 585 }
568 586
569 value *= G4UniformRand(); 587 value *= G4UniformRand();
570 588
571 i = n; 589 i = n;
572 while (i > 0) 590 while (i > 0)
573 { 591 {
574 i--; 592 i--;
575 593
576 if (values[i] > value) 594 if (values[i] > value)
577 break; 595 break;
578 596
579 value -= values[i]; 597 value -= values[i];
580 } 598 }
581 599
582 delete[] values; 600 delete[] values;
583 601
584 return i; 602 return i;
585 } 603 }
586 604
587 //....oooOO0OOooo........oooOO0OOooo........oo 605 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
588 606
589 G4double G4DNADingfelderChargeDecreaseModel::S 607 G4double G4DNADingfelderChargeDecreaseModel::Sum(G4double k,
590 608 const G4ParticleDefinition* particleDefinition)
591 { 609 {
592 G4int particleTypeIndex = 0; 610 G4int particleTypeIndex = 0;
>> 611 G4DNAGenericIonsManager* instance;
>> 612 instance = G4DNAGenericIonsManager::Instance();
593 613
594 if (particleDefinition == protonDef) << 614 if (particleDefinition == G4Proton::ProtonDefinition())
595 particleTypeIndex = 0; 615 particleTypeIndex = 0;
596 616
597 if (particleDefinition == alphaPlusPlusDef) << 617 if (particleDefinition == instance->GetIon("alpha++"))
598 particleTypeIndex = 1; 618 particleTypeIndex = 1;
599 619
600 if (particleDefinition == alphaPlusDef) << 620 if (particleDefinition == instance->GetIon("alpha+"))
601 particleTypeIndex = 2; 621 particleTypeIndex = 2;
602 622
603 G4double totalCrossSection = 0.; 623 G4double totalCrossSection = 0.;
604 624
605 for (G4int i = 0; i < numberOfPartialCrossSe 625 for (G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++)
606 { 626 {
607 totalCrossSection += PartialCrossSection(k 627 totalCrossSection += PartialCrossSection(k, i, particleDefinition);
608 } 628 }
609 629
610 return totalCrossSection; 630 return totalCrossSection;
611 } 631 }
612 632
613 633