| Geant4 Cross Reference |
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 #include "globals.hh"
28 #include "G4PhysicalConstants.hh"
29 #include "G4BetaDecayType.hh"
30 #include "G4BetaDecayCorrections.hh"
31 #include "G4Pow.hh"
32
33 G4BetaDecayCorrections::G4BetaDecayCorrections(const G4int theZ, const G4int theA)
34 : Z(theZ), A(theA)
35 {
36 // alphaZ = fine_structure_const*std::abs(Z);
37 alphaZ = fine_structure_const*Z;
38
39 // Nuclear radius in units of hbar/m_e/c
40 G4double a13 = G4Pow::GetInstance()->Z13(A);
41 Rnuc = 0.5*fine_structure_const*a13;
42
43 // Electron screening potential in units of electron mass
44 V0 = 1.13*fine_structure_const*fine_structure_const
45 *std::pow(std::abs(Z), 4./3.);
46
47 gamma0 = std::sqrt(1. - alphaZ*alphaZ);
48
49 // Largest allowed value of im argument in ModSquared
50 // imMax = std::log(DBL_MAX)/pi;
51 imMax = 200.; // actual value = 225.931, but use 200 to be safe
52 // G4cout << " imMax = " << imMax << G4endl;
53
54 // Coefficients for gamma function with real argument
55 gc[0] = -0.1010678;
56 gc[1] = 0.4245549;
57 gc[2] = -0.6998588;
58 gc[3] = 0.9512363;
59 gc[4] = -0.5748646;
60 gc[5] = 1.0;
61 }
62
63
64 G4double G4BetaDecayCorrections::FermiFunction(const G4double& W)
65 {
66 // Calculate the relativistic Fermi function. Argument W is the
67 // total electron energy in units of electron mass.
68 // Ref: E. Feenberg, G. Trigg, Reviews of Modern Physics, 22(1950)
69
70 G4double Wprime;
71 if (Z < 0) {
72 Wprime = W + V0;
73 } else {
74 Wprime = W - V0;
75 // if (Wprime < 1.) Wprime = W;
76 if (Wprime <= 1.00001) Wprime = 1.00001;
77 }
78
79 G4double p_e = std::sqrt(Wprime*Wprime - 1.);
80 G4double eta = alphaZ*Wprime/p_e;
81 G4double epieta = std::exp(pi*eta);
82 G4double realGamma = Gamma(2.*gamma0+1);
83 G4double mod2Gamma = ModSquared(gamma0, eta);
84
85 // Fermi function
86 G4double factor1 = 2*(1+gamma0)*mod2Gamma/realGamma/realGamma;
87 G4double factor2 = epieta*std::pow(2*p_e*Rnuc, 2*(gamma0-1) );
88
89 // Electron screening factor
90 G4double factor3 = (Wprime/W)*std::sqrt( (Wprime*Wprime - 1.)/(W*W - 1.) );
91
92 return factor1*factor2*factor3;
93 }
94
95
96 G4double
97 G4BetaDecayCorrections::ModSquared(const G4double& re, G4double im)
98 {
99 // Calculate the squared modulus of the Gamma function
100 // with complex argument (re, im) using approximation B
101 // of Wilkinson, Nucl. Instr. & Meth. 82, 122 (1970).
102 // Here, choose N = 1 in Wilkinson's notation for approximation B
103
104 im = std::max(std::min(im, imMax), -imMax);
105 G4double factor1 = std::pow( (1+re)*(1+re) + im*im, re+0.5);
106 G4double factor2 = std::exp(2*im * std::atan(im/(1+re)));
107 G4double factor3 = std::exp(2*(1+re));
108 G4double factor4 = 2.*pi;
109 G4double factor5 = std::exp( (1+re)/( (1+re)*(1+re) + im*im)/6 );
110 G4double factor6 = re*re + im*im;
111 return factor1*factor4*factor5/factor2/factor3/factor6;
112 }
113
114
115 G4double G4BetaDecayCorrections::Gamma(const G4double& arg)
116 {
117 // Use recursion relation to get argument < 1
118 G4double fac = 1.0;
119 G4double x = arg - 1.;
120
121 G4int loop = 0;
122 G4ExceptionDescription ed;
123 ed << " While count exceeded " << G4endl;
124 while (x > 1.0) { /* Loop checking, 01.09.2015, D.Wright */
125 fac *= x;
126 x -= 1.0;
127 loop++;
128 if (loop > 1000) {
129 G4Exception("G4BetaDecayCorrections::Gamma()", "HAD_RDM_100", JustWarning, ed);
130 break;
131 }
132 }
133
134 // Calculation of Gamma function with real argument
135 // 0 < arg < 1 using polynomial from Abramowitz and Stegun
136 G4double sum = gc[0];
137 for (G4int i = 1; i < 6; i++) sum = sum*x + gc[i];
138
139 return sum*fac;
140 }
141
142
143 G4double
144 G4BetaDecayCorrections::ShapeFactor(const G4BetaDecayType& bdt,
145 const G4double& p_e, const G4double& e_nu)
146 {
147 G4double twoPR = 2.*p_e*Rnuc;
148 G4double factor(1.);
149
150 switch (bdt)
151 {
152 case (allowed) :
153 break;
154
155 case (firstForbidden) :
156 {
157 // Parameters for 1st forbidden shape determined from 210Bi data
158 // Not valid for other 1st forbidden nuclei
159 G4double c1 = 0.578;
160 G4double c2 = 28.466;
161 G4double c3 = -0.658;
162
163 G4double w = std::sqrt(1. + p_e*p_e);
164 factor = 1. + c1*w + c2/w + c3*w*w;
165 }
166 break;
167
168 case (uniqueFirstForbidden) :
169 {
170 G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
171 G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
172 G4double gamterm1 = Gamma(2.*gamma0+1.)/Gamma(2.*gamma1+1.);
173 G4double term1 = e_nu*e_nu*(1. + gamma0)/6.;
174 G4double term2 = 12.*(2. + gamma1)*p_e*p_e
175 *std::pow(twoPR, 2.*(gamma1-gamma0-1) )
176 *gamterm1*gamterm1
177 *ModSquared(gamma1, eta)/ModSquared(gamma0, eta);
178 factor = term1 + term2;
179 }
180 break;
181
182 case (secondForbidden) :
183 break;
184
185 case (uniqueSecondForbidden) :
186 {
187 G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
188 G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
189 G4double gamma2 = std::sqrt(9. - alphaZ*alphaZ);
190 G4double gamterm0 = Gamma(2.*gamma0+1.);
191 G4double gamterm1 = gamterm0/Gamma(2.*gamma1+1.);
192 G4double gamterm2 = gamterm0/Gamma(2.*gamma2+1.);
193 G4double term1 = e_nu*e_nu*e_nu*e_nu*(1. + gamma0)/60.;
194
195 G4double term2 = 4.*(2. + gamma1)*e_nu*e_nu*p_e*p_e
196 *std::pow(twoPR, 2.*(gamma1-gamma0-1.) )
197 *gamterm1*gamterm1
198 *ModSquared(gamma1, eta)/ModSquared(gamma0, eta);
199
200 G4double term3 = 180.*(3.+gamma2)*p_e*p_e*p_e*p_e
201 *std::pow(twoPR, 2.*(gamma2-gamma0-2) )
202 *gamterm2*gamterm2
203 *ModSquared(gamma2, eta)/ModSquared(gamma0, eta);
204
205 factor = term1 + term2 + term3;
206 }
207 break;
208
209 case (thirdForbidden) :
210 break;
211
212 case (uniqueThirdForbidden) :
213 {
214 G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
215 G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
216 G4double gamma2 = std::sqrt(9. - alphaZ*alphaZ);
217 G4double gamma3 = std::sqrt(16. - alphaZ*alphaZ);
218 G4double gamterm0 = Gamma(2.*gamma0+1.);
219 G4double gamterm1 = gamterm0/Gamma(2.*gamma1+1.);
220 G4double gamterm2 = gamterm0/Gamma(2.*gamma2+1.);
221 G4double gamterm3 = gamterm0/Gamma(2.*gamma3+1.);
222
223 G4double term1 = e_nu*e_nu*e_nu*e_nu*e_nu*e_nu*(1. + gamma0)/1260.;
224
225 G4double term2 = 2.*(2. + gamma1)*e_nu*e_nu*e_nu*e_nu*p_e*p_e
226 *std::pow(twoPR, 2.*(gamma1-gamma0-1.) )
227 *gamterm1*gamterm1
228 *ModSquared(gamma1, eta)/ModSquared(gamma0, eta)/5.;
229
230 G4double term3 = 60.*(3.+gamma2)*p_e*p_e*p_e*p_e*e_nu*e_nu
231 *std::pow(twoPR, 2.*(gamma2-gamma0-2.) )
232 *gamterm2*gamterm2
233 *ModSquared(gamma2, eta)/ModSquared(gamma0, eta);
234
235 G4double term4 = 2240.*p_e*p_e*p_e*p_e*p_e*p_e*(4. + gamma3)
236 *std::pow(twoPR, 2.*(gamma3-gamma0-3.) )
237 *gamterm3*gamterm3
238 *ModSquared(gamma3, eta)/ModSquared(gamma0, eta);
239
240 factor = term1 + term2 + term3 + term4;
241 }
242 break;
243
244 default:
245 G4Exception("G4BetaDecayCorrections::ShapeFactor()","HAD_RDM_010",
246 JustWarning,
247 "Transition not yet implemented - using allowed shape");
248 break;
249 }
250
251 return factor;
252 }
253
254
255