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

Last change on this file since 1058 was 1055, checked in by garnier, 17 years ago

maj sur la beta de geant 4.9.3
  • Property svn:executable set to *
File size: 17.0 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.9 2009/05/02 15:07:47 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
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 G4cout << G4endl;
74 G4cout << "*******************************************************************************" << G4endl;
75 G4cout << "*******************************************************************************" << G4endl;
76 G4cout << " The class G4FinalStateIonisationRudd is NOT SUPPORTED ANYMORE. " << G4endl;
77 G4cout << " It will be REMOVED with the next major release of Geant4. " << G4endl;
78 G4cout << " Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
79 G4cout << "*******************************************************************************" << G4endl;
80 G4cout << "*******************************************************************************" << G4endl;
81 G4cout << G4endl;
82}
83
84//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85
86G4FinalStateIonisationRudd::~G4FinalStateIonisationRudd()
87{}
88
89//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90
91const G4FinalStateProduct& G4FinalStateIonisationRudd::GenerateFinalState(const G4Track& track, const G4Step& /* step */)
92{
93 product.Clear();
94
95 const G4DynamicParticle* particle = track.GetDynamicParticle();
96
97 G4double lowLim = lowEnergyLimitDefault;
98 G4double highLim = highEnergyLimitDefault;
99
100 G4double k = particle->GetKineticEnergy();
101
102 const G4String& particleName = particle->GetDefinition()->GetParticleName();
103
104 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
105 pos1 = lowEnergyLimit.find(particleName);
106
107 if (pos1 != lowEnergyLimit.end())
108 {
109 lowLim = pos1->second;
110 }
111
112 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
113 pos2 = highEnergyLimit.find(particleName);
114
115 if (pos2 != highEnergyLimit.end())
116 {
117 highLim = pos2->second;
118 }
119
120 if (k >= lowLim && k <= highLim)
121 {
122 G4ParticleDefinition* definition = particle->GetDefinition();
123 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
124 G4double particleMass = definition->GetPDGMass();
125 G4double totalEnergy = k + particleMass;
126 G4double pSquare = k*(totalEnergy+particleMass);
127 G4double totalMomentum = std::sqrt(pSquare);
128
129 const G4String& particleName = definition->GetParticleName();
130
131 G4int ionizationShell = cross.RandomSelect(k,particleName);
132
133 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
134
135 G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
136
137 G4double cosTheta = 0.;
138 G4double phi = 0.;
139 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi);
140
141 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
142 G4double dirX = sinTheta*std::cos(phi);
143 G4double dirY = sinTheta*std::sin(phi);
144 G4double dirZ = cosTheta;
145 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
146 deltaDirection.rotateUz(primaryDirection);
147
148 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
149
150 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
151 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
152 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
153 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
154 finalPx /= finalMomentum;
155 finalPy /= finalMomentum;
156 finalPz /= finalMomentum;
157
158 product.ModifyPrimaryParticle(finalPx,finalPy,finalPz,k-bindingEnergy-secondaryKinetic);
159 product.AddEnergyDeposit(bindingEnergy);
160
161 G4DynamicParticle* aElectron = new G4DynamicParticle(G4Electron::Electron(),deltaDirection,secondaryKinetic);
162 product.AddSecondary(aElectron);
163 }
164
165 if (k < lowLim)
166 {
167 product.KillPrimaryParticle();
168 }
169
170 return product;
171}
172
173//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
174
175G4double G4FinalStateIonisationRudd::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
176 G4double k,
177 G4int shell)
178{
179 G4double maximumKineticEnergyTransfer = 0.;
180
181 G4DNAGenericIonsManager *instance;
182 instance = G4DNAGenericIonsManager::Instance();
183
184 if (particleDefinition == G4Proton::ProtonDefinition()
185 || particleDefinition == instance->GetIon("hydrogen"))
186 {
187 maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
188 }
189
190 if (particleDefinition == instance->GetIon("helium")
191 || particleDefinition == instance->GetIon("alpha+")
192 || particleDefinition == instance->GetIon("alpha++"))
193 {
194 maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
195 }
196
197 G4double crossSectionMaximum = 0.;
198
199 for(G4double value=waterStructure.IonisationEnergy(shell); value<=4.*waterStructure.IonisationEnergy(shell) ; value+=0.1*eV)
200 {
201 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
202 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
203 }
204
205 G4double secElecKinetic = 0.;
206
207 do
208 {
209 secElecKinetic = G4UniformRand() * maximumKineticEnergyTransfer;
210 } while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
211 k,
212 secElecKinetic+waterStructure.IonisationEnergy(shell),
213 shell));
214
215 return(secElecKinetic);
216}
217
218//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
219
220
221void G4FinalStateIonisationRudd::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
222 G4double k,
223 G4double secKinetic,
224 G4double & cosTheta,
225 G4double & phi )
226{
227 G4DNAGenericIonsManager *instance;
228 instance = G4DNAGenericIonsManager::Instance();
229
230 G4double maxSecKinetic = 0.;
231
232 if (particleDefinition == G4Proton::ProtonDefinition()
233 || particleDefinition == instance->GetIon("hydrogen"))
234 {
235 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
236 }
237
238 if (particleDefinition == instance->GetIon("helium")
239 || particleDefinition == instance->GetIon("alpha+")
240 || particleDefinition == instance->GetIon("alpha++"))
241 {
242 maxSecKinetic = 4.* (0.511 / 3728) * k;
243 }
244
245 phi = twopi * G4UniformRand();
246 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
247}
248
249//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
250
251
252G4double G4FinalStateIonisationRudd::DifferentialCrossSection(G4ParticleDefinition* particleDefinition,
253 G4double k,
254 G4double energyTransfer,
255 G4int ionizationLevelIndex)
256{
257 // Shells ids are 0 1 2 3 4 (4 is k shell)
258 // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
259 // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
260 //
261 // ds S F1(nu) + w * F2(nu)
262 // ---- = G(k) * ---- -------------------------------------------
263 // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
264 //
265 // w is the secondary electron kinetic Energy in eV
266 //
267 // All the other parameters can be found in Rudd's Papers
268 //
269 // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
270 // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
271 //
272
273 const G4int j=ionizationLevelIndex;
274
275 G4double A1 ;
276 G4double B1 ;
277 G4double C1 ;
278 G4double D1 ;
279 G4double E1 ;
280 G4double A2 ;
281 G4double B2 ;
282 G4double C2 ;
283 G4double D2 ;
284 G4double alphaConst ;
285
286 if (j == 4)
287 {
288 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
289 A1 = 1.25;
290 B1 = 0.5;
291 C1 = 1.00;
292 D1 = 1.00;
293 E1 = 3.00;
294 A2 = 1.10;
295 B2 = 1.30;
296 C2 = 1.00;
297 D2 = 0.00;
298 alphaConst = 0.66;
299 }
300 else
301 {
302 //Data For Liquid Water from Dingfelder (Protons in Water)
303 A1 = 1.02;
304 B1 = 82.0;
305 C1 = 0.45;
306 D1 = -0.80;
307 E1 = 0.38;
308 A2 = 1.07;
309 B2 = 14.6;
310 C2 = 0.60;
311 D2 = 0.04;
312 alphaConst = 0.64;
313 }
314
315 const G4double n = 2.;
316 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
317
318 G4DNAGenericIonsManager* instance;
319 instance = G4DNAGenericIonsManager::Instance();
320
321 G4double wBig = (energyTransfer - waterStructure.IonisationEnergy(ionizationLevelIndex));
322 G4double w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
323 G4double Ry = 13.6*eV;
324
325 G4double tau = 0.;
326
327 if (particleDefinition == G4Proton::ProtonDefinition()
328 || particleDefinition == instance->GetIon("hydrogen"))
329 {
330 tau = (electron_mass_c2/proton_mass_c2) * k ;
331 }
332
333 if ( particleDefinition == instance->GetIon("helium")
334 || particleDefinition == instance->GetIon("alpha+")
335 || particleDefinition == instance->GetIon("alpha++"))
336 {
337 tau = (0.511/3728.) * k ;
338 }
339
340 G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
341 G4double v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
342 G4double v = std::sqrt(v2);
343 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.IonisationEnergy(ionizationLevelIndex)));
344
345 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
346 G4double L2 = C2*std::pow(v,(D2));
347 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
348 G4double H2 = (A2/v2) + (B2/(v2*v2));
349
350 G4double F1 = L1+H1;
351 G4double F2 = (L2*H2)/(L2+H2);
352
353 G4double sigma = CorrectionFactor(particleDefinition, k/eV)
354 * Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex))
355 * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
356
357 if ( particleDefinition == G4Proton::ProtonDefinition()
358 || particleDefinition == instance->GetIon("hydrogen")
359 )
360 {
361 return(sigma);
362 }
363
364 if (particleDefinition == instance->GetIon("alpha++") )
365 {
366 slaterEffectiveCharge[0]=0.;
367 slaterEffectiveCharge[1]=0.;
368 slaterEffectiveCharge[2]=0.;
369 sCoefficient[0]=0.;
370 sCoefficient[1]=0.;
371 sCoefficient[2]=0.;
372 }
373
374 if (particleDefinition == instance->GetIon("alpha+") )
375 {
376 slaterEffectiveCharge[0]=2.0;
377 slaterEffectiveCharge[1]=1.15;
378 slaterEffectiveCharge[2]=1.15;
379 sCoefficient[0]=0.7;
380 sCoefficient[1]=0.15;
381 sCoefficient[2]=0.15;
382 }
383
384 if (particleDefinition == instance->GetIon("helium") )
385 {
386 slaterEffectiveCharge[0]=1.7;
387 slaterEffectiveCharge[1]=1.15;
388 slaterEffectiveCharge[2]=1.15;
389 sCoefficient[0]=0.5;
390 sCoefficient[1]=0.25;
391 sCoefficient[2]=0.25;
392 }
393
394 if ( particleDefinition == instance->GetIon("helium")
395 || particleDefinition == instance->GetIon("alpha+")
396 || particleDefinition == instance->GetIon("alpha++")
397 )
398 {
399 sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
400
401 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
402
403 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
404 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
405 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
406
407 return zEff * zEff * sigma ;
408 }
409
410 return 0;
411}
412
413//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
414
415G4double G4FinalStateIonisationRudd::S_1s(G4double t,
416 G4double energyTransferred,
417 G4double slaterEffectiveChg,
418 G4double shellNumber)
419{
420 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
421 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
422
423 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
424 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
425
426 return value;
427}
428
429//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
430
431G4double G4FinalStateIonisationRudd::S_2s(G4double t,
432 G4double energyTransferred,
433 G4double slaterEffectiveChg,
434 G4double shellNumber)
435{
436 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
437 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
438
439 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
440 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
441
442 return value;
443
444}
445
446//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
447
448G4double G4FinalStateIonisationRudd::S_2p(G4double t,
449 G4double energyTransferred,
450 G4double slaterEffectiveChg,
451 G4double shellNumber)
452{
453 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
454 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
455
456 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
457 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
458
459 return value;
460}
461
462//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
463
464G4double G4FinalStateIonisationRudd::R(G4double t,
465 G4double energyTransferred,
466 G4double slaterEffectiveChg,
467 G4double shellNumber)
468{
469 // tElectron = m_electron / m_alpha * t
470 // Dingfelder, in Chattanooga 2005 proceedings, p 4
471
472 G4double tElectron = 0.511/3728. * t;
473 G4double value = 2. * tElectron * slaterEffectiveChg / (energyTransferred * shellNumber);
474
475 return value;
476}
477
478//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
479
480G4double G4FinalStateIonisationRudd::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k)
481{
482 G4DNAGenericIonsManager *instance;
483 instance = G4DNAGenericIonsManager::Instance();
484
485 if (particleDefinition == G4Proton::Proton())
486 {
487 return(1.);
488 }
489 else
490 if (particleDefinition == instance->GetIon("hydrogen"))
491 {
492 G4double value = (std::log(k/eV)-4.2)/0.5;
493 return((0.8/(1+std::exp(value))) + 0.9);
494 }
495 else
496 {
497 return(1.);
498 }
499}
Note: See TracBrowser for help on using the repository browser.