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

Last change on this file since 1007 was 1007, checked in by garnier, 15 years ago

update to geant4.9.2
  • 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.