source: trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionExcitationMillerGreenPartial.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: 10.9 KB
RevLine 
[819]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: G4CrossSectionExcitationMillerGreenPartial.cc,v 1.1 2007/11/08 19:57:23 pia Exp $
28// GEANT4 tag $Name: $
29//
30// Contact Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
31//
32// Reference: TNS Geant4-DNA paper
33// Reference for implementation model: NIM. 155, pp. 145-156, 1978
34
35// History:
36// -----------
37// Date Name Modification
38// 28 Apr 2007 M.G. Pia Created in compliance with design described in TNS paper
39//
40// -------------------------------------------------------------------
41
42// Class description:
43// Geant4-DNA Cross total cross section for electron elastic scattering in water
44// Reference: TNS Geant4-DNA paper
45// S. Chauvie et al., Geant4 physics processes for microdosimetry simulation:
46// design foundation and implementation of the first set of models,
47// IEEE Trans. Nucl. Sci., vol. 54, no. 6, Dec. 2007.
48// Further documentation available from http://www.ge.infn.it/geant4/dna
49
50// -------------------------------------------------------------------
51
52
53#include "G4CrossSectionExcitationMillerGreenPartial.hh"
54#include "G4Track.hh"
55#include "G4DynamicParticle.hh"
56#include "G4ParticleDefinition.hh"
57#include "G4Proton.hh"
58#include "G4DNAGenericIonsManager.hh"
59#include "G4CrossSectionExcitationEmfietzoglouPartial.hh"
60#include "Randomize.hh"
61#include <deque>
62
63
64G4CrossSectionExcitationMillerGreenPartial::G4CrossSectionExcitationMillerGreenPartial()
65{
66
67 nLevels = waterExcitation.NumberOfLevels();
68 //G4cout << "Water excitation energy: number of levels = " << nLevels << G4endl;
69
70 //PROTON
71 kineticEnergyCorrection[0] = 1.;
72 slaterEffectiveCharge[0][0] = 0.;
73 slaterEffectiveCharge[1][0] = 0.;
74 slaterEffectiveCharge[2][0] = 0.;
75 sCoefficient[0][0] = 0.;
76 sCoefficient[1][0] = 0.;
77 sCoefficient[2][0] = 0.;
78
79 //ALPHA++
80 kineticEnergyCorrection[1] = 0.9382723/3.727417;
81 slaterEffectiveCharge[0][1]=0.;
82 slaterEffectiveCharge[1][1]=0.;
83 slaterEffectiveCharge[2][1]=0.;
84 sCoefficient[0][1]=0.;
85 sCoefficient[1][1]=0.;
86 sCoefficient[2][1]=0.;
87
88 // ALPHA+
89 kineticEnergyCorrection[2] = 0.9382723/3.727417;
90 slaterEffectiveCharge[0][2]=2.0;
91 slaterEffectiveCharge[1][2]=1.15;
92 slaterEffectiveCharge[2][2]=1.15;
93 sCoefficient[0][2]=0.7;
94 sCoefficient[1][2]=0.15;
95 sCoefficient[2][2]=0.15;
96
97 // HELIUM
98 kineticEnergyCorrection[3] = 0.9382723/3.727417;
99 slaterEffectiveCharge[0][3]=1.7;
100 slaterEffectiveCharge[1][3]=1.15;
101 slaterEffectiveCharge[2][3]=1.15;
102 sCoefficient[0][3]=0.5;
103 sCoefficient[1][3]=0.25;
104 sCoefficient[2][3]=0.25;
105
106}
107
108G4CrossSectionExcitationMillerGreenPartial::~G4CrossSectionExcitationMillerGreenPartial()
109{ }
110
111G4double G4CrossSectionExcitationMillerGreenPartial::CrossSection(G4double k, G4int excitationLevel,
112 const G4ParticleDefinition* particleDefinition)
113{
114 // ( ( z * aj ) ^ omegaj ) * ( t - ej ) ^ nu
115 // sigma(t) = zEff^2 * sigma0 * --------------------------------------------
116 // jj ^ ( omegaj + nu ) + t ^ ( omegaj + nu )
117 //
118 // where t is the kinetic energy corrected by Helium mass over proton mass for Helium ions
119 //
120 // zEff is:
121 // 1 for protons
122 // 2 for alpha++
123 // and 2 - c1 S_1s - c2 S_2s - c3 S_2p for alpha+ and He
124 //
125 // Dingfelder et al., RPC 59, 255-275, 2000 from Miller and Green (1973)
126 // Formula (34) and Table 2
127
128 const G4double sigma0(1.E+8 * barn);
129 const G4double nu(1.);
130 const G4double aj[]={876.*eV, 2084.* eV, 1373.*eV, 692.*eV, 900.*eV};
131 const G4double jj[]={19820.*eV, 23490.*eV, 27770.*eV, 30830.*eV, 33080.*eV};
132 const G4double omegaj[]={0.85, 0.88, 0.88, 0.78, 0.78};
133
134 G4int particleTypeIndex = 0;
135 G4DNAGenericIonsManager* instance;
136 instance = G4DNAGenericIonsManager::Instance();
137
138 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
139 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
140 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
141 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=3;
142
143 G4double tCorrected;
144 tCorrected = k * kineticEnergyCorrection[particleTypeIndex];
145
146 // Assume that the material is water; proper algorithm to calculate correctly for any material to be inserted here
147 G4int z = 10;
148
149 G4double numerator;
150 numerator = std::pow(z * aj[excitationLevel], omegaj[excitationLevel]) *
151 std::pow(tCorrected - waterExcitation.ExcitationEnergy(excitationLevel), nu);
152
153 G4double power;
154 power = omegaj[excitationLevel] + nu;
155
156 G4double denominator;
157 denominator = std::pow(jj[excitationLevel], power) + std::pow(tCorrected, power);
158
159 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
160
161 zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, waterExcitation.ExcitationEnergy(excitationLevel), slaterEffectiveCharge[0][particleTypeIndex], 1.) +
162 sCoefficient[1][particleTypeIndex] * S_2s(k, waterExcitation.ExcitationEnergy(excitationLevel), slaterEffectiveCharge[1][particleTypeIndex], 2.) +
163 sCoefficient[2][particleTypeIndex] * S_2p(k, waterExcitation.ExcitationEnergy(excitationLevel), slaterEffectiveCharge[2][particleTypeIndex], 2.) );
164
165 G4double cross = sigma0 * zEff * zEff * numerator / denominator;
166
167 return cross;
168}
169
170G4int G4CrossSectionExcitationMillerGreenPartial::RandomSelect(G4double k,
171 const G4ParticleDefinition* particle)
172{
173 G4int i = nLevels;
174 G4double value = 0.;
175 std::deque<double> values;
176
177 // ---- MGP ---- The following algorithm is wrong: it works is the cross section
178 // is a monotone increasing function.
179 // The algorithm should be corrected by building the cumulative function
180 // of the cross section and comparing a random number in the range 0-1 against
181 // the cumulative value at each bin
182
183 G4DNAGenericIonsManager *instance;
184 instance = G4DNAGenericIonsManager::Instance();
185
186 // ELECTRON CORRECTION
187
188 if ( particle == instance->GetIon("alpha++"))
189 { while (i > 0)
190 {
191 i--;
192 G4double partial = CrossSection(k,i,particle);
193 values.push_front(partial);
194 value += partial;
195 }
196
197 value *= G4UniformRand();
198
199 i = nLevels;
200
201 while (i > 0)
202 {
203 i--;
204 if (values[i] > value) return i;
205 value -= values[i];
206 }
207 }
208
209 //
210 // add ONE or TWO electron-water excitation for alpha+ and helium
211
212 if ( particle == instance->GetIon("alpha+")
213 ||
214 particle == instance->GetIon("helium")
215 )
216 {
217
218 while (i>0)
219 {
220 i--;
221
222 G4CrossSectionExcitationEmfietzoglouPartial* excitationXS =
223 new G4CrossSectionExcitationEmfietzoglouPartial();
224
225 G4double sigmaExcitation=0;
226 if (k*0.511/3728 > 7.4*eV && k*0.511/3728 < 10*keV) sigmaExcitation = excitationXS->CrossSection(k*0.511/3728,i);
227
228 G4double partial = CrossSection(k,i,particle);
229 if (particle == instance->GetIon("alpha+")) partial = CrossSection(k,i,particle) + sigmaExcitation;
230 if (particle == instance->GetIon("helium")) partial = CrossSection(k,i,particle) + 2*sigmaExcitation;
231 values.push_front(partial);
232 value += partial;
233 delete excitationXS;
234 }
235
236 value*=G4UniformRand();
237
238 i=5;
239 while (i>0)
240 {
241 i--;
242
243 if (values[i]>value) return i;
244
245 value-=values[i];
246 }
247 }
248 //
249 return 0;
250}
251
252G4double G4CrossSectionExcitationMillerGreenPartial::Sum(G4double k, const G4ParticleDefinition* particle)
253{
254 G4double totalCrossSection = 0.;
255
256 for (G4int i=0; i<nLevels; i++)
257 {
258 totalCrossSection += CrossSection(k,i,particle);
259 }
260 return totalCrossSection;
261}
262
263
264G4double G4CrossSectionExcitationMillerGreenPartial::S_1s(G4double t,
265 G4double energyTransferred,
266 G4double slaterEffectiveCharge,
267 G4double shellNumber)
268{
269 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
270 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
271
272 G4double r = R(t, energyTransferred, slaterEffectiveCharge, shellNumber);
273 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
274
275 return value;
276}
277
278
279G4double G4CrossSectionExcitationMillerGreenPartial::S_2s(G4double t,
280 G4double energyTransferred,
281 G4double slaterEffectiveCharge,
282 G4double shellNumber)
283{
284 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
285 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
286
287 G4double r = R(t, energyTransferred, slaterEffectiveCharge, shellNumber);
288 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
289
290 return value;
291
292}
293
294
295G4double G4CrossSectionExcitationMillerGreenPartial::S_2p(G4double t,
296 G4double energyTransferred,
297 G4double slaterEffectiveCharge,
298 G4double shellNumber)
299{
300 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
301 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
302
303 G4double r = R(t, energyTransferred, slaterEffectiveCharge, shellNumber);
304 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
305
306 return value;
307}
308
309
310G4double G4CrossSectionExcitationMillerGreenPartial::R(G4double t,
311 G4double energyTransferred,
312 G4double slaterEffectiveCharge,
313 G4double shellNumber)
314{
315 // tElectron = m_electron / m_alpha * t
316 // Hardcoded in Riccardo's implementation; to be corrected
317 // Dingfelder, in Chattanooga 2005 proceedings, p 4
318
319 G4double tElectron = 0.511/3728. * t;
320 G4double value = 2. * tElectron * slaterEffectiveCharge / (energyTransferred * shellNumber);
321
322 return value;
323}
324
Note: See TracBrowser for help on using the repository browser.