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