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

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
  • 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.11 2009/06/11 15:47:08 mantero Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-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.