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