source: trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateElasticBrennerZaider.cc@ 1014

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

update to geant4.9.2
File size: 6.2 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: G4FinalStateElasticBrennerZaider.cc,v 1.8 2008/12/05 11:58:16 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28
29#include "G4FinalStateElasticBrennerZaider.hh"
30
31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33G4FinalStateElasticBrennerZaider::G4FinalStateElasticBrennerZaider()
34{
35 lowEnergyLimit = 8.23 * eV; // SI : i/o of 7.4 * eV;
36 highEnergyLimit = 200 * eV;
37
38 betaCoeff.push_back(7.51525);
39 betaCoeff.push_back(-0.41912);
40 betaCoeff.push_back(7.2017E-3);
41 betaCoeff.push_back(-4.646E-5);
42 betaCoeff.push_back(1.02897E-7);
43
44 deltaCoeff.push_back(2.9612);
45 deltaCoeff.push_back(-0.26376);
46 deltaCoeff.push_back(4.307E-3);
47 deltaCoeff.push_back(-2.6895E-5);
48 deltaCoeff.push_back(5.83505E-8);
49
50 gamma035_10Coeff.push_back(-1.7013);
51 gamma035_10Coeff.push_back(-1.48284);
52 gamma035_10Coeff.push_back(0.6331);
53 gamma035_10Coeff.push_back(-0.10911);
54 gamma035_10Coeff.push_back(8.358E-3);
55 gamma035_10Coeff.push_back(-2.388E-4);
56
57 gamma10_100Coeff.push_back(-3.32517);
58 gamma10_100Coeff.push_back(0.10996);
59 gamma10_100Coeff.push_back(-4.5255E-3);
60 gamma10_100Coeff.push_back(5.8372E-5);
61 gamma10_100Coeff.push_back(-2.4659E-7);
62
63 gamma100_200Coeff.push_back(2.4775E-2);
64 gamma100_200Coeff.push_back(-2.96264E-5);
65 gamma100_200Coeff.push_back(-1.20655E-7);
66}
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69
70G4FinalStateElasticBrennerZaider::~G4FinalStateElasticBrennerZaider()
71{}
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74
75const G4FinalStateProduct& G4FinalStateElasticBrennerZaider::GenerateFinalState(const G4Track& track, const G4Step& /* step */)
76{
77 product.Clear();
78
79 G4double k = track.GetDynamicParticle()->GetKineticEnergy();
80
81 if ( k>= lowEnergyLimit && k<=highEnergyLimit )
82 {
83 G4double cosTheta = RandomizeCosTheta(k);
84
85 G4double phi = 2. * pi * G4UniformRand();
86
87 G4ThreeVector zVers = track.GetDynamicParticle()->GetMomentumDirection();
88 G4ThreeVector xVers = zVers.orthogonal();
89 G4ThreeVector yVers = zVers.cross(xVers);
90
91 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
92 G4double yDir = xDir;
93 xDir *= std::cos(phi);
94 yDir *= std::sin(phi);
95
96 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
97
98 product.ModifyPrimaryParticle(zPrimeVers,k);
99 }
100
101 if (k<lowEnergyLimit)
102 {
103 product.KillPrimaryParticle();
104 }
105
106 return product;
107}
108
109//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110
111G4double G4FinalStateElasticBrennerZaider::RandomizeCosTheta(G4double k)
112{
113 // d sigma_el 1 beta(K)
114 // ------------ (K) ~ --------------------------------- + ---------------------------------
115 // d Omega (1 + 2 gamma(K) - cos(theta))^2 (1 + 2 delta(K) + cos(theta))^2
116 //
117 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
118 //
119 // Phys. Med. Biol. 29 N.4 (1983) 443-447
120
121 // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
122
123 k /= eV;
124
125 G4double beta = std::exp(CalculatePolynomial(k,betaCoeff));
126 G4double delta = std::exp(CalculatePolynomial(k,deltaCoeff));
127 G4double gamma;
128
129 if (k > 100.)
130 {
131 gamma = CalculatePolynomial(k, gamma100_200Coeff);
132 // Only in this case it is not the exponent of the polynomial
133 }
134 else
135 {
136 if (k>10)
137 {
138 gamma = std::exp(CalculatePolynomial(k, gamma10_100Coeff));
139 }
140 else
141 {
142 gamma = std::exp(CalculatePolynomial(k, gamma035_10Coeff));
143 }
144 }
145
146 G4double oneOverMax = 1. / (1./(4.*gamma*gamma) + beta/( (2.+2.*delta)*(2.+2.*delta) ));
147
148 G4double cosTheta = 0.;
149 G4double leftDenominator = 0.;
150 G4double rightDenominator = 0.;
151 G4double fCosTheta = 0.;
152
153 do
154 {
155 cosTheta = 2. * G4UniformRand() - 1.;
156 leftDenominator = (1. + 2.*gamma - cosTheta);
157 rightDenominator = (1. + 2.*delta + cosTheta);
158 if ( (leftDenominator * rightDenominator) != 0. )
159 {
160 fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
161 }
162 }
163 while (fCosTheta < G4UniformRand());
164
165 return cosTheta;
166}
167
168//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
169
170G4double G4FinalStateElasticBrennerZaider::CalculatePolynomial(G4double k, std::vector<G4double>& vec)
171{
172 // Sum_{i=0}^{size-1} vector_i k^i
173 //
174 // Phys. Med. Biol. 29 N.4 (1983) 443-447
175
176 G4double result = 0.;
177 size_t size = vec.size();
178
179 while (size>0)
180 {
181 size--;
182
183 result *= k;
184 result += vec[size];
185 }
186
187 return result;
188}
189
Note: See TracBrowser for help on using the repository browser.