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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 7.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: G4FinalStateElasticBrennerZaider.cc,v 1.11 2009/06/11 15:47:08 mantero Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
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 G4cout << G4endl;
68 G4cout << "*******************************************************************************" << G4endl;
69 G4cout << "*******************************************************************************" << G4endl;
70 G4cout << " The class G4FinalStateElasticBrennerZaider is NOT SUPPORTED ANYMORE. " << G4endl;
71 G4cout << " It will be REMOVED with the next major release of Geant4. " << G4endl;
72 G4cout << " Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
73 G4cout << "*******************************************************************************" << G4endl;
74 G4cout << "*******************************************************************************" << G4endl;
75 G4cout << G4endl;
76}
77
78//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79
80G4FinalStateElasticBrennerZaider::~G4FinalStateElasticBrennerZaider()
81{}
82
83//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84
85const G4FinalStateProduct& G4FinalStateElasticBrennerZaider::GenerateFinalState(const G4Track& track, const G4Step& /* step */)
86{
87 product.Clear();
88
89 G4double k = track.GetDynamicParticle()->GetKineticEnergy();
90
91 if ( k>= lowEnergyLimit && k<=highEnergyLimit )
92 {
93 G4double cosTheta = RandomizeCosTheta(k);
94
95 G4double phi = 2. * pi * G4UniformRand();
96
97 G4ThreeVector zVers = track.GetDynamicParticle()->GetMomentumDirection();
98 G4ThreeVector xVers = zVers.orthogonal();
99 G4ThreeVector yVers = zVers.cross(xVers);
100
101 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
102 G4double yDir = xDir;
103 xDir *= std::cos(phi);
104 yDir *= std::sin(phi);
105
106 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
107
108 product.ModifyPrimaryParticle(zPrimeVers,k);
109 }
110
111 if (k<lowEnergyLimit)
112 {
113 product.KillPrimaryParticle();
114 }
115
116 return product;
117}
118
119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120
121G4double G4FinalStateElasticBrennerZaider::RandomizeCosTheta(G4double k)
122{
123 // d sigma_el 1 beta(K)
124 // ------------ (K) ~ --------------------------------- + ---------------------------------
125 // d Omega (1 + 2 gamma(K) - cos(theta))^2 (1 + 2 delta(K) + cos(theta))^2
126 //
127 // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
128 //
129 // Phys. Med. Biol. 29 N.4 (1983) 443-447
130
131 // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
132
133 k /= eV;
134
135 G4double beta = std::exp(CalculatePolynomial(k,betaCoeff));
136 G4double delta = std::exp(CalculatePolynomial(k,deltaCoeff));
137 G4double gamma;
138
139 if (k > 100.)
140 {
141 gamma = CalculatePolynomial(k, gamma100_200Coeff);
142 // Only in this case it is not the exponent of the polynomial
143 }
144 else
145 {
146 if (k>10)
147 {
148 gamma = std::exp(CalculatePolynomial(k, gamma10_100Coeff));
149 }
150 else
151 {
152 gamma = std::exp(CalculatePolynomial(k, gamma035_10Coeff));
153 }
154 }
155
156 G4double oneOverMax = 1. / (1./(4.*gamma*gamma) + beta/( (2.+2.*delta)*(2.+2.*delta) ));
157
158 G4double cosTheta = 0.;
159 G4double leftDenominator = 0.;
160 G4double rightDenominator = 0.;
161 G4double fCosTheta = 0.;
162
163 do
164 {
165 cosTheta = 2. * G4UniformRand() - 1.;
166 leftDenominator = (1. + 2.*gamma - cosTheta);
167 rightDenominator = (1. + 2.*delta + cosTheta);
168 if ( (leftDenominator * rightDenominator) != 0. )
169 {
170 fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
171 }
172 }
173 while (fCosTheta < G4UniformRand());
174
175 return cosTheta;
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
179
180G4double G4FinalStateElasticBrennerZaider::CalculatePolynomial(G4double k, std::vector<G4double>& vec)
181{
182 // Sum_{i=0}^{size-1} vector_i k^i
183 //
184 // Phys. Med. Biol. 29 N.4 (1983) 443-447
185
186 G4double result = 0.;
187 size_t size = vec.size();
188
189 while (size>0)
190 {
191 size--;
192
193 result *= k;
194 result += vec[size];
195 }
196
197 return result;
198}
199
Note: See TracBrowser for help on using the repository browser.