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

Last change on this file since 830 was 819, checked in by garnier, 17 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.