source: trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateIonisationRudd.cc@ 992

Last change on this file since 992 was 991, checked in by garnier, 17 years ago

update
  • Property svn:executable set to *
File size: 16.3 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4FinalStateIonisationRudd.cc,v 1.8 2008/08/20 14:51:48 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28
29#include "G4FinalStateIonisationRudd.hh"
30
31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33G4FinalStateIonisationRudd::G4FinalStateIonisationRudd()
34{
35 lowEnergyLimitDefault = 100 * eV;
36 highEnergyLimitDefault = 100 * MeV;
37
38 G4DNAGenericIonsManager *instance;
39 instance = G4DNAGenericIonsManager::Instance();
40
41 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
42 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
43 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
44 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
45 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
46
47 G4String proton;
48 G4String hydrogen;
49 G4String alphaPlusPlus;
50 G4String alphaPlus;
51 G4String helium;
52
53 proton = protonDef->GetParticleName();
54 lowEnergyLimit[proton] = 100. * eV;
55 highEnergyLimit[proton] = 500. * keV;
56
57 hydrogen = hydrogenDef->GetParticleName();
58 lowEnergyLimit[hydrogen] = 100. * eV;
59 highEnergyLimit[hydrogen] = 100. * MeV;
60
61 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
62 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
63 highEnergyLimit[alphaPlusPlus] = 10. * MeV;
64
65 alphaPlus = alphaPlusDef->GetParticleName();
66 lowEnergyLimit[alphaPlus] = 1. * keV;
67 highEnergyLimit[alphaPlus] = 10. * MeV;
68
69 helium = heliumDef->GetParticleName();
70 lowEnergyLimit[helium] = 1. * keV;
71 highEnergyLimit[helium] = 10. * MeV;
72}
73
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75
76G4FinalStateIonisationRudd::~G4FinalStateIonisationRudd()
77{}
78
79//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80
81const G4FinalStateProduct& G4FinalStateIonisationRudd::GenerateFinalState(const G4Track& track, const G4Step& /* step */)
82{
83 product.Clear();
84
85 const G4DynamicParticle* particle = track.GetDynamicParticle();
86
87 G4double lowLim = lowEnergyLimitDefault;
88 G4double highLim = highEnergyLimitDefault;
89
90 G4double k = particle->GetKineticEnergy();
91
92 const G4String& particleName = particle->GetDefinition()->GetParticleName();
93
94 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
95 pos1 = lowEnergyLimit.find(particleName);
96
97 if (pos1 != lowEnergyLimit.end())
98 {
99 lowLim = pos1->second;
100 }
101
102 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
103 pos2 = highEnergyLimit.find(particleName);
104
105 if (pos2 != highEnergyLimit.end())
106 {
107 highLim = pos2->second;
108 }
109
110 if (k >= lowLim && k <= highLim)
111 {
112 G4ParticleDefinition* definition = particle->GetDefinition();
113 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
114 G4double particleMass = definition->GetPDGMass();
115 G4double totalEnergy = k + particleMass;
116 G4double pSquare = k*(totalEnergy+particleMass);
117 G4double totalMomentum = std::sqrt(pSquare);
118
119 const G4String& particleName = definition->GetParticleName();
120
121 G4int ionizationShell = cross.RandomSelect(k,particleName);
122
123 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
124
125 G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
126
127 G4double cosTheta = 0.;
128 G4double phi = 0.;
129 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi);
130
131 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
132 G4double dirX = sinTheta*std::cos(phi);
133 G4double dirY = sinTheta*std::sin(phi);
134 G4double dirZ = cosTheta;
135 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
136 deltaDirection.rotateUz(primaryDirection);
137
138 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
139
140 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
141 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
142 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
143 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
144 finalPx /= finalMomentum;
145 finalPy /= finalMomentum;
146 finalPz /= finalMomentum;
147
148 product.ModifyPrimaryParticle(finalPx,finalPy,finalPz,k-bindingEnergy-secondaryKinetic);
149 product.AddEnergyDeposit(bindingEnergy);
150
151 G4DynamicParticle* aElectron = new G4DynamicParticle(G4Electron::Electron(),deltaDirection,secondaryKinetic);
152 product.AddSecondary(aElectron);
153 }
154
155 if (k < lowLim)
156 {
157 product.KillPrimaryParticle();
158 }
159
160 return product;
161}
162
163//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164
165G4double G4FinalStateIonisationRudd::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
166 G4double k,
167 G4int shell)
168{
169 G4double maximumKineticEnergyTransfer = 0.;
170
171 G4DNAGenericIonsManager *instance;
172 instance = G4DNAGenericIonsManager::Instance();
173
174 if (particleDefinition == G4Proton::ProtonDefinition()
175 || particleDefinition == instance->GetIon("hydrogen"))
176 {
177 maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
178 }
179
180 if (particleDefinition == instance->GetIon("helium")
181 || particleDefinition == instance->GetIon("alpha+")
182 || particleDefinition == instance->GetIon("alpha++"))
183 {
184 maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
185 }
186
187 G4double crossSectionMaximum = 0.;
188
189 for(G4double value=waterStructure.IonisationEnergy(shell); value<=4.*waterStructure.IonisationEnergy(shell) ; value+=0.1*eV)
190 {
191 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
192 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
193 }
194
195 G4double secElecKinetic = 0.;
196
197 do
198 {
199 secElecKinetic = G4UniformRand() * maximumKineticEnergyTransfer;
200 } while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
201 k,
202 secElecKinetic+waterStructure.IonisationEnergy(shell),
203 shell));
204
205 return(secElecKinetic);
206}
207
208//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209
210
211void G4FinalStateIonisationRudd::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
212 G4double k,
213 G4double secKinetic,
214 G4double & cosTheta,
215 G4double & phi )
216{
217 G4DNAGenericIonsManager *instance;
218 instance = G4DNAGenericIonsManager::Instance();
219
220 G4double maxSecKinetic = 0.;
221
222 if (particleDefinition == G4Proton::ProtonDefinition()
223 || particleDefinition == instance->GetIon("hydrogen"))
224 {
225 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
226 }
227
228 if (particleDefinition == instance->GetIon("helium")
229 || particleDefinition == instance->GetIon("alpha+")
230 || particleDefinition == instance->GetIon("alpha++"))
231 {
232 maxSecKinetic = 4.* (0.511 / 3728) * k;
233 }
234
235 phi = twopi * G4UniformRand();
236 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
237}
238
239//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
240
241
242G4double G4FinalStateIonisationRudd::DifferentialCrossSection(G4ParticleDefinition* particleDefinition,
243 G4double k,
244 G4double energyTransfer,
245 G4int ionizationLevelIndex)
246{
247 // Shells ids are 0 1 2 3 4 (4 is k shell)
248 // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
249 // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
250 //
251 // ds S F1(nu) + w * F2(nu)
252 // ---- = G(k) * ---- -------------------------------------------
253 // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
254 //
255 // w is the secondary electron kinetic Energy in eV
256 //
257 // All the other parameters can be found in Rudd's Papers
258 //
259 // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
260 // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
261 //
262
263 const G4int j=ionizationLevelIndex;
264
265 G4double A1 ;
266 G4double B1 ;
267 G4double C1 ;
268 G4double D1 ;
269 G4double E1 ;
270 G4double A2 ;
271 G4double B2 ;
272 G4double C2 ;
273 G4double D2 ;
274 G4double alphaConst ;
275
276 if (j == 4)
277 {
278 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
279 A1 = 1.25;
280 B1 = 0.5;
281 C1 = 1.00;
282 D1 = 1.00;
283 E1 = 3.00;
284 A2 = 1.10;
285 B2 = 1.30;
286 C2 = 1.00;
287 D2 = 0.00;
288 alphaConst = 0.66;
289 }
290 else
291 {
292 //Data For Liquid Water from Dingfelder (Protons in Water)
293 A1 = 1.02;
294 B1 = 82.0;
295 C1 = 0.45;
296 D1 = -0.80;
297 E1 = 0.38;
298 A2 = 1.07;
299 B2 = 14.6;
300 C2 = 0.60;
301 D2 = 0.04;
302 alphaConst = 0.64;
303 }
304
305 const G4double n = 2.;
306 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
307
308 G4DNAGenericIonsManager* instance;
309 instance = G4DNAGenericIonsManager::Instance();
310
311 G4double wBig = (energyTransfer - waterStructure.IonisationEnergy(ionizationLevelIndex));
312 G4double w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
313 G4double Ry = 13.6*eV;
314
315 G4double tau = 0.;
316
317 if (particleDefinition == G4Proton::ProtonDefinition()
318 || particleDefinition == instance->GetIon("hydrogen"))
319 {
320 tau = (electron_mass_c2/proton_mass_c2) * k ;
321 }
322
323 if ( particleDefinition == instance->GetIon("helium")
324 || particleDefinition == instance->GetIon("alpha+")
325 || particleDefinition == instance->GetIon("alpha++"))
326 {
327 tau = (0.511/3728.) * k ;
328 }
329
330 G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
331 G4double v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
332 G4double v = std::sqrt(v2);
333 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.IonisationEnergy(ionizationLevelIndex)));
334
335 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
336 G4double L2 = C2*std::pow(v,(D2));
337 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
338 G4double H2 = (A2/v2) + (B2/(v2*v2));
339
340 G4double F1 = L1+H1;
341 G4double F2 = (L2*H2)/(L2+H2);
342
343 G4double sigma = CorrectionFactor(particleDefinition, k/eV)
344 * Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
345 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
346
347 if ( particleDefinition == G4Proton::ProtonDefinition()
348 || particleDefinition == instance->GetIon("hydrogen")
349 )
350 {
351 return(sigma);
352 }
353
354 if (particleDefinition == instance->GetIon("alpha++") )
355 {
356 slaterEffectiveCharge[0]=0.;
357 slaterEffectiveCharge[1]=0.;
358 slaterEffectiveCharge[2]=0.;
359 sCoefficient[0]=0.;
360 sCoefficient[1]=0.;
361 sCoefficient[2]=0.;
362 }
363
364 if (particleDefinition == instance->GetIon("alpha+") )
365 {
366 slaterEffectiveCharge[0]=2.0;
367 slaterEffectiveCharge[1]=1.15;
368 slaterEffectiveCharge[2]=1.15;
369 sCoefficient[0]=0.7;
370 sCoefficient[1]=0.15;
371 sCoefficient[2]=0.15;
372 }
373
374 if (particleDefinition == instance->GetIon("helium") )
375 {
376 slaterEffectiveCharge[0]=1.7;
377 slaterEffectiveCharge[1]=1.15;
378 slaterEffectiveCharge[2]=1.15;
379 sCoefficient[0]=0.5;
380 sCoefficient[1]=0.25;
381 sCoefficient[2]=0.25;
382 }
383
384 if ( particleDefinition == instance->GetIon("helium")
385 || particleDefinition == instance->GetIon("alpha+")
386 || particleDefinition == instance->GetIon("alpha++")
387 )
388 {
389 sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
390
391 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
392
393 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
394 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
395 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
396
397 return zEff * zEff * sigma ;
398 }
399
400 return 0;
401}
402
403//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404
405G4double G4FinalStateIonisationRudd::S_1s(G4double t,
406 G4double energyTransferred,
407 G4double slaterEffectiveChg,
408 G4double shellNumber)
409{
410 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
411 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
412
413 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
414 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
415
416 return value;
417}
418
419//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
420
421G4double G4FinalStateIonisationRudd::S_2s(G4double t,
422 G4double energyTransferred,
423 G4double slaterEffectiveChg,
424 G4double shellNumber)
425{
426 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
427 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
428
429 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
430 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
431
432 return value;
433
434}
435
436//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
437
438G4double G4FinalStateIonisationRudd::S_2p(G4double t,
439 G4double energyTransferred,
440 G4double slaterEffectiveChg,
441 G4double shellNumber)
442{
443 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
444 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
445
446 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
447 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
448
449 return value;
450}
451
452//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
453
454G4double G4FinalStateIonisationRudd::R(G4double t,
455 G4double energyTransferred,
456 G4double slaterEffectiveChg,
457 G4double shellNumber)
458{
459 // tElectron = m_electron / m_alpha * t
460 // Dingfelder, in Chattanooga 2005 proceedings, p 4
461
462 G4double tElectron = 0.511/3728. * t;
463 G4double value = 2. * tElectron * slaterEffectiveChg / (energyTransferred * shellNumber);
464
465 return value;
466}
467
468//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
469
470G4double G4FinalStateIonisationRudd::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k)
471{
472 G4DNAGenericIonsManager *instance;
473 instance = G4DNAGenericIonsManager::Instance();
474
475 if (particleDefinition == G4Proton::Proton())
476 {
477 return(1.);
478 }
479 else
480 if (particleDefinition == instance->GetIon("hydrogen"))
481 {
482 G4double value = (std::log(k/eV)-4.2)/0.5;
483 return((0.8/(1+std::exp(value))) + 0.9);
484 }
485 else
486 {
487 return(1.);
488 }
489}
Note: See TracBrowser for help on using the repository browser.