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

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

import all except CVS
  • Property svn:executable set to *
File size: 17.4 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//
27// $Id: G4FinalStateIonisationRudd.cc,v 1.5 2007/11/26 17:27:09 pia Exp $
28// GEANT4 tag $Name:  $
29//
30// Contact Author: Sebastien Incerti (incerti@cenbg.in2p3.fr)
31//                 Maria Grazia Pia  (Maria.Grazia.Pia@cern.ch)
32//
33///
34// Reference: TNS Geant4-DNA paper
35// Reference for implementation model: NIM. 155, pp. 145-156, 1978
36//
37// History:
38// -----------
39// Date         Name              Modification
40// 28 Apr 2007  M.G. Pia          Created in compliance with design described in TNS paper
41//    Nov 2007  S. Incerti        Implementation
42// 26 Nov 2007  MGP               Cleaned up std::
43//
44// -------------------------------------------------------------------
45
46// Class description:
47// Reference: TNS Geant4-DNA paper
48// S. Chauvie et al., Geant4 physics processes for microdosimetry simulation:
49// design foundation and implementation of the first set of models,
50// IEEE Trans. Nucl. Sci., vol. 54, no. 6, Dec. 2007.
51// Further documentation available from http://www.ge.infn.it/geant4/dna
52
53// -------------------------------------------------------------------
54
55
56#include "G4FinalStateIonisationRudd.hh"
57#include "G4Track.hh"
58#include "G4Step.hh"
59#include "G4DynamicParticle.hh"
60#include "Randomize.hh"
61
62#include "G4ParticleTypes.hh"
63#include "G4ParticleDefinition.hh"
64#include "G4Electron.hh"
65#include "G4Proton.hh"
66#include "G4SystemOfUnits.hh"
67#include "G4ParticleMomentum.hh"
68#include "G4DNAGenericIonsManager.hh"
69
70
71G4FinalStateIonisationRudd::G4FinalStateIonisationRudd()
72{
73  name = "IonisationBorn";
74  // Default energy limits (defined for protection against anomalous behaviour only)
75  lowEnergyLimitDefault = 100 * eV;
76  highEnergyLimitDefault = 100 * MeV;
77
78  G4DNAGenericIonsManager *instance;
79  instance = G4DNAGenericIonsManager::Instance();
80
81  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
82  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
83  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
84  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
85  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
86
87  G4String proton;
88  G4String hydrogen;
89  G4String alphaPlusPlus;
90  G4String alphaPlus;
91  G4String helium;
92
93  proton = protonDef->GetParticleName();
94  lowEnergyLimit[proton] = 100. * eV;
95  highEnergyLimit[proton] = 500. * keV;
96
97  hydrogen = hydrogenDef->GetParticleName();
98  lowEnergyLimit[hydrogen] = 100. * eV;
99  highEnergyLimit[hydrogen] = 100. * MeV;
100
101  alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
102  lowEnergyLimit[alphaPlusPlus] = 1. * keV;
103  highEnergyLimit[alphaPlusPlus] = 10. * MeV;
104
105  alphaPlus = alphaPlusDef->GetParticleName();
106  lowEnergyLimit[alphaPlus] = 1. * keV;
107  highEnergyLimit[alphaPlus] = 10. * MeV;
108
109  helium = heliumDef->GetParticleName();
110  lowEnergyLimit[helium] = 1. * keV;
111  highEnergyLimit[helium] = 10. * MeV;
112}
113
114
115G4FinalStateIonisationRudd::~G4FinalStateIonisationRudd()
116{ }
117
118
119
120const G4FinalStateProduct& G4FinalStateIonisationRudd::GenerateFinalState(const G4Track& track, const G4Step& /* step */)
121{
122  // Clear previous secondaries, energy deposit and particle kill status
123  product.Clear();
124
125  const G4DynamicParticle* particle = track.GetDynamicParticle();
126
127  G4double lowLim = lowEnergyLimitDefault;
128  G4double highLim = highEnergyLimitDefault;
129
130  G4double k = particle->GetKineticEnergy();
131
132  const G4String& particleName = particle->GetDefinition()->GetParticleName();
133
134  // Retrieve energy limits for the current particle type
135
136  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
137  pos1 = lowEnergyLimit.find(particleName);
138
139  // Lower limit
140  if (pos1 != lowEnergyLimit.end())
141    {
142      lowLim = pos1->second;
143    }
144
145  // Upper limit
146  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
147  pos2 = highEnergyLimit.find(particleName);
148
149  if (pos2 != highEnergyLimit.end())
150    {
151      highLim = pos2->second;
152    }
153
154  // Verify that the current track is within the energy limits of validity of the cross section model
155
156  if (k >= lowLim && k <= highLim)
157    {
158      // Kinetic energy of primary particle
159
160      G4ParticleDefinition* definition = particle->GetDefinition();
161      G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
162      G4double particleMass = definition->GetPDGMass();
163      G4double totalEnergy = k + particleMass;
164      G4double pSquare = k*(totalEnergy+particleMass);
165      G4double totalMomentum = std::sqrt(pSquare);
166
167      const G4String& particleName = definition->GetParticleName();
168 
169      G4int ionizationShell = cross.RandomSelect(k,particleName);
170 
171      G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
172 
173      G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
174
175      G4double cosTheta = 0.;
176      G4double phi = 0.; 
177      RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi);
178
179      G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
180      G4double dirX = sinTheta*std::cos(phi);
181      G4double dirY = sinTheta*std::sin(phi);
182      G4double dirZ = cosTheta;
183      G4ThreeVector deltaDirection(dirX,dirY,dirZ);
184      deltaDirection.rotateUz(primaryDirection);
185
186      G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
187
188      // Primary Particle Direction
189      G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
190      G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
191      G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
192      G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
193      finalPx /= finalMomentum;
194      finalPy /= finalMomentum;
195      finalPz /= finalMomentum;
196
197      product.ModifyPrimaryParticle(finalPx,finalPy,finalPz,k-bindingEnergy-secondaryKinetic);
198      product.AddEnergyDeposit(bindingEnergy);
199
200      G4DynamicParticle* aElectron = new G4DynamicParticle(G4Electron::Electron(),deltaDirection,secondaryKinetic);
201      product.AddSecondary(aElectron);
202    }
203
204  if (k < lowLim) {product.KillPrimaryParticle();product.AddEnergyDeposit(k);}
205 
206  return product;
207}
208
209
210
211G4double G4FinalStateIonisationRudd::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition, 
212                                                                    G4double k, 
213                                                                    G4int shell)
214{
215  G4double maximumKineticEnergyTransfer = 0.;
216 
217  G4DNAGenericIonsManager *instance;
218  instance = G4DNAGenericIonsManager::Instance();
219
220  if (particleDefinition == G4Proton::ProtonDefinition() 
221      || particleDefinition == instance->GetIon("hydrogen"))
222         
223    { 
224      maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
225    }
226
227  if (particleDefinition == instance->GetIon("helium") 
228      || particleDefinition == instance->GetIon("alpha+")
229      || particleDefinition == instance->GetIon("alpha++"))
230    { 
231      maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
232    }
233
234  G4double crossSectionMaximum = 0.;
235  for(G4double value=waterStructure.IonisationEnergy(shell); value<=4.*waterStructure.IonisationEnergy(shell) ; value+=0.1*eV)
236    {
237      G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
238      if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
239    }
240  G4double secElecKinetic = 0.;
241  do{
242    secElecKinetic = G4UniformRand() * maximumKineticEnergyTransfer;
243  } while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition, 
244                                                                         k,
245                                                                         secElecKinetic+waterStructure.IonisationEnergy(shell),
246                                                                         shell));
247
248  return(secElecKinetic);
249}
250
251
252void G4FinalStateIonisationRudd::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition, 
253                                                                   G4double k, 
254                                                                   G4double secKinetic, 
255                                                                   G4double cosTheta, 
256                                                                   G4double phi )
257{
258  G4DNAGenericIonsManager *instance;
259  instance = G4DNAGenericIonsManager::Instance();
260
261  G4double maxSecKinetic = 0.;
262 
263  if (particleDefinition == G4Proton::ProtonDefinition() 
264      || particleDefinition == instance->GetIon("hydrogen")) 
265    { 
266      maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
267    }
268 
269  if (particleDefinition == instance->GetIon("helium") 
270      || particleDefinition == instance->GetIon("alpha+")
271      || particleDefinition == instance->GetIon("alpha++"))
272    { 
273      maxSecKinetic = 4.* (0.511 / 3728) * k;
274    }
275 
276  phi = twopi * G4UniformRand();
277  cosTheta = std::sqrt(secKinetic / maxSecKinetic);
278}
279
280
281G4double G4FinalStateIonisationRudd::DifferentialCrossSection(G4ParticleDefinition* particleDefinition, 
282                                                              G4double k, 
283                                                              G4double energyTransfer, 
284                                                              G4int ionizationLevelIndex)
285{
286  // Shells ids are 0 1 2 3 4 (4 is k shell)
287  // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
288  //             that the secondary kinetic energy is w = energyTransfer - bindingEnergy
289  //
290  //   ds            S                F1(nu) + w * F2(nu)
291  //  ---- = G(k) * ----     -------------------------------------------
292  //   dw            Bj       (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
293  //
294  // w is the secondary electron kinetic Energy in eV
295  //
296  // All the other parameters can be found in Rudd's Papers
297  //
298  // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
299  // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
300  //
301
302  const G4int j=ionizationLevelIndex;
303
304  G4double A1 ; 
305  G4double B1 ; 
306  G4double C1 ; 
307  G4double D1 ; 
308  G4double E1 ;
309  G4double A2 ; 
310  G4double B2 ; 
311  G4double C2 ; 
312  G4double D2 ; 
313  G4double alphaConst ;
314
315  if (j == 4) 
316    {
317      //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
318      A1 = 1.25; 
319      B1 = 0.5; 
320      C1 = 1.00; 
321      D1 = 1.00; 
322      E1 = 3.00; 
323      A2 = 1.10; 
324      B2 = 1.30;
325      C2 = 1.00; 
326      D2 = 0.00; 
327      alphaConst = 0.66;
328    }
329  else 
330    {
331      //Data For Liquid Water from Dingfelder (Protons in Water)
332      A1 = 1.02; 
333      B1 = 82.0; 
334      C1 = 0.45; 
335      D1 = -0.80; 
336      E1 = 0.38; 
337      A2 = 1.07; 
338      B2 = 14.6;
339      C2 = 0.60; 
340      D2 = 0.04; 
341      alphaConst = 0.64;
342    }
343 
344  const G4double n = 2.;
345  const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
346
347  //const G4double I[5]={12.61*eV, 14.73*eV, 18.55*eV, 32.2*eV, 539.7*eV}; // for water Vapor
348  //const G4double energyConstant[]={10.79*eV, 13.39*eV, 16.05*eV, 32.30*eV, 539.*eV};
349
350  G4DNAGenericIonsManager* instance;
351  instance = G4DNAGenericIonsManager::Instance();
352
353  G4double wBig = (energyTransfer - waterStructure.IonisationEnergy(ionizationLevelIndex));
354  G4double w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
355  G4double Ry = 13.6*eV;
356
357  G4double tau = 0.;
358
359  if (particleDefinition == G4Proton::ProtonDefinition() 
360      || particleDefinition == instance->GetIon("hydrogen")) 
361    {
362      tau = (electron_mass_c2/proton_mass_c2) * k ;
363    }
364   
365  if ( particleDefinition == instance->GetIon("helium") 
366       || particleDefinition == instance->GetIon("alpha+")
367       || particleDefinition == instance->GetIon("alpha++")) 
368    {
369      tau = (0.511/3728.) * k ;
370    }
371 
372  G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
373  G4double v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
374  G4double v = std::sqrt(v2);
375  G4double wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.IonisationEnergy(ionizationLevelIndex)));
376
377  G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
378  G4double L2 = C2*std::pow(v,(D2));
379  G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
380  G4double H2 = (A2/v2) + (B2/(v2*v2));
381
382  G4double F1 = L1+H1;
383  G4double F2 = (L2*H2)/(L2+H2);
384
385  G4double sigma = CorrectionFactor(particleDefinition, k/eV) 
386    * Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) 
387    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
388   
389  if (    particleDefinition == G4Proton::ProtonDefinition() 
390          || particleDefinition == instance->GetIon("hydrogen")
391          ) 
392    {
393      return(sigma);
394    }
395
396  // ------------
397 
398  if (particleDefinition == instance->GetIon("alpha++") ) 
399    {
400      slaterEffectiveCharge[0]=0.;
401      slaterEffectiveCharge[1]=0.;
402      slaterEffectiveCharge[2]=0.;
403      sCoefficient[0]=0.;
404      sCoefficient[1]=0.;
405      sCoefficient[2]=0.;
406    }
407
408  if (particleDefinition == instance->GetIon("alpha+") ) 
409    {
410      slaterEffectiveCharge[0]=2.0;
411      slaterEffectiveCharge[1]=1.15;
412      slaterEffectiveCharge[2]=1.15;
413      sCoefficient[0]=0.7;
414      sCoefficient[1]=0.15;
415      sCoefficient[2]=0.15;
416    }
417
418  if (particleDefinition == instance->GetIon("helium") ) 
419    {
420      slaterEffectiveCharge[0]=1.7;
421      slaterEffectiveCharge[1]=1.15;
422      slaterEffectiveCharge[2]=1.15;
423      sCoefficient[0]=0.5;
424      sCoefficient[1]=0.25;
425      sCoefficient[2]=0.25;
426    }
427 
428  if (    particleDefinition == instance->GetIon("helium") 
429          || particleDefinition == instance->GetIon("alpha+")
430          || particleDefinition == instance->GetIon("alpha++")
431          ) 
432    {
433      sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
434   
435      G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
436 
437      zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
438                sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
439                sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
440           
441      return zEff * zEff * sigma ;
442    } 
443 
444  return 0;
445}
446
447G4double G4FinalStateIonisationRudd::S_1s(G4double t, 
448                                          G4double energyTransferred, 
449                                          G4double slaterEffectiveChg, 
450                                          G4double shellNumber)
451{
452  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
453  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
454 
455  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
456  G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
457 
458  return value;
459}
460
461
462
463G4double G4FinalStateIonisationRudd::S_2s(G4double t,
464                                          G4double energyTransferred, 
465                                          G4double slaterEffectiveChg, 
466                                          G4double shellNumber)
467{
468  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
469  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
470
471  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
472  G4double value =  1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
473
474  return value;
475 
476}
477
478
479
480G4double G4FinalStateIonisationRudd::S_2p(G4double t, 
481                                          G4double energyTransferred,
482                                          G4double slaterEffectiveChg, 
483                                          G4double shellNumber)
484{
485  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
486  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
487
488  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
489  G4double value =  1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r  + 1.);
490
491  return value;
492}
493
494
495
496G4double G4FinalStateIonisationRudd::R(G4double t,
497                                       G4double energyTransferred,
498                                       G4double slaterEffectiveChg,
499                                       G4double shellNumber) 
500{
501  // tElectron = m_electron / m_alpha * t
502  // Hardcoded in Riccardo's implementation; to be corrected
503  // Dingfelder, in Chattanooga 2005 proceedings, p 4
504
505  G4double tElectron = 0.511/3728. * t;
506  G4double value = 2. * tElectron * slaterEffectiveChg / (energyTransferred * shellNumber);
507 
508  return value;
509}
510
511
512
513G4double G4FinalStateIonisationRudd::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k) 
514{
515  G4DNAGenericIonsManager *instance;
516  instance = G4DNAGenericIonsManager::Instance();
517
518  if (particleDefinition == G4Proton::Proton()) 
519    {
520      return(1.);
521    }
522  else 
523    if (particleDefinition == instance->GetIon("hydrogen")) 
524      { 
525        G4double value = (std::log(k/eV)-4.2)/0.5;
526        return((0.8/(1+std::exp(value))) + 0.9);
527      }
528    else 
529      {   
530        return(1.);
531      }
532}
Note: See TracBrowser for help on using the repository browser.