Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BetaDecayCorrections.cc
Go to the documentation of this file.
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 
32 G4BetaDecayCorrections::G4BetaDecayCorrections(G4int theZ, G4int theA)
33  : Z(theZ), A(theA)
34 {
35  // alphaZ = fine_structure_const*std::abs(Z);
36  alphaZ = fine_structure_const*Z;
37 
38  // Nuclear radius in units of hbar/m_e/c
39  Rnuc = 0.5*fine_structure_const*std::pow(A, 0.33333);
40 
41  // Electron screening potential in units of electrom mass
42  V0 = 1.13*fine_structure_const*fine_structure_const
43  *std::pow(std::abs(Z), 1.33333);
44 
45  gamma0 = std::sqrt(1. - alphaZ*alphaZ);
46 
47  // Coefficients for gamma function with real argument
48  gc[0] = -0.1010678;
49  gc[1] = 0.4245549;
50  gc[2] = -0.6998588;
51  gc[3] = 0.9512363;
52  gc[4] = -0.5748646;
53  gc[5] = 1.0;
54 }
55 
56 
57 G4double G4BetaDecayCorrections::FermiFunction(const G4double& W)
58 {
59  // Calculate the relativistic Fermi function. Argument W is the
60  // total electron energy in units of electron mass.
61 
62  G4double Wprime;
63  if (Z < 0) {
64  Wprime = W + V0;
65  } else {
66  Wprime = W - V0;
67 // if (Wprime < 1.) Wprime = W;
68  if (Wprime <= 1.00001) Wprime = 1.00001;
69  }
70 
71  G4double p_e = std::sqrt(Wprime*Wprime - 1.);
72  G4double eta = alphaZ*Wprime/p_e;
73  G4double epieta = std::exp(pi*eta);
74  G4double realGamma = Gamma(2.*gamma0+1);
75  G4double mod2Gamma = ModSquared(gamma0, eta);
76 
77  // Fermi function
78  G4double factor1 = 2*(1+gamma0)*mod2Gamma/realGamma/realGamma;
79  G4double factor2 = epieta*std::pow(2*p_e*Rnuc, 2*(gamma0-1) );
80 
81  // Electron screening factor
82  G4double factor3 = (Wprime/W)*std::sqrt( (Wprime*Wprime - 1.)/(W*W - 1.) );
83 
84  return factor1*factor2*factor3;
85 }
86 
87 
88 G4double
89 G4BetaDecayCorrections::ModSquared(const G4double& re, const G4double& im)
90 {
91  // Calculate the squared modulus of the Gamma function
92  // with complex argument (re, im) using approximation B
93  // of Wilkinson, Nucl. Instr. & Meth. 82, 122 (1970).
94  // Here, choose N = 1 in Wilkinson's notation for approximation B
95 
96  G4double factor1 = std::pow( (1+re)*(1+re) + im*im, re+0.5);
97  G4double factor2 = std::exp(2*im * std::atan(im/(1+re)));
98  G4double factor3 = std::exp(2*(1+re));
99  G4double factor4 = 2.*pi;
100  G4double factor5 = std::exp( (1+re)/( (1+re)*(1+re) + im*im)/6 );
101  G4double factor6 = re*re + im*im;
102  return factor1*factor4*factor5/factor2/factor3/factor6;
103 }
104 
105 
106 G4double G4BetaDecayCorrections::Gamma(const G4double& arg)
107 {
108  // Use recursion relation to get argument < 1
109  G4double fac = 1.0;
110  G4double x = arg - 1.;
111  while (x > 1.0) {
112  fac *= x;
113  x -= 1.0;
114  }
115 
116  // Calculation of Gamma function with real argument
117  // 0 < arg < 1 using polynomial from Abramowitz and Stegun
118  G4double sum = gc[0];
119  for (G4int i = 1; i < 6; i++) sum = sum*x + gc[i];
120 
121  return sum*fac;
122 }
123 
124 
125 G4double
126 G4BetaDecayCorrections::ShapeFactor(const G4BetaDecayType& bdt,
127  const G4double& p_e, const G4double& e_nu)
128 {
129  G4double twoPR = 2.*p_e*Rnuc;
130  G4double factor(1.);
131 
132  switch (bdt)
133  {
134  case (allowed) :
135  break;
136 
137  case (firstForbidden) :
138  {
139  // Parameters for 1st forbidden shape determined from 210Bi data
140  // Not valid for other 1st forbidden nuclei
141  G4double c1 = 0.578;
142  G4double c2 = 28.466;
143  G4double c3 = -0.658;
144 
145  G4double w = std::sqrt(1. + p_e*p_e);
146  factor = 1. + c1*w + c2/w + c3*w*w;
147  }
148  break;
149 
150  case (uniqueFirstForbidden) :
151  {
152  G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
153  G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
154  G4double gamterm1 = Gamma(2.*gamma0+1.)/Gamma(2.*gamma1+1.);
155  G4double term1 = e_nu*e_nu*(1. + gamma0)/6.;
156  G4double term2 = 12.*(2. + gamma1)*p_e*p_e
157  *std::pow(twoPR, 2.*(gamma1-gamma0-1) )
158  *gamterm1*gamterm1
159  *ModSquared(gamma1, eta)/ModSquared(gamma0, eta);
160  factor = term1 + term2;
161  }
162  break;
163 
164  case (uniqueSecondForbidden) :
165  {
166  G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
167  G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
168  G4double gamma2 = std::sqrt(9. - alphaZ*alphaZ);
169  G4double gamterm0 = Gamma(2.*gamma0+1.);
170  G4double gamterm1 = gamterm0/Gamma(2.*gamma1+1.);
171  G4double gamterm2 = gamterm0/Gamma(2.*gamma2+1.);
172  G4double term1 = e_nu*e_nu*e_nu*e_nu*(1. + gamma0)/60.;
173 
174  G4double term2 = 4.*(2. + gamma1)*e_nu*e_nu*p_e*p_e
175  *std::pow(twoPR, 2.*(gamma1-gamma0-1.) )
176  *gamterm1*gamterm1
177  *ModSquared(gamma1, eta)/ModSquared(gamma0, eta);
178 
179  G4double term3 = 180.*(3.+gamma2)*p_e*p_e*p_e*p_e
180  *std::pow(twoPR, 2.*(gamma2-gamma0-2) )
181  *gamterm2*gamterm2
182  *ModSquared(gamma2, eta)/ModSquared(gamma0, eta);
183 
184  factor = term1 + term2 + term3;
185  }
186  break;
187 
188  case (uniqueThirdForbidden) :
189  {
190  G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
191  G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
192  G4double gamma2 = std::sqrt(9. - alphaZ*alphaZ);
193  G4double gamma3 = std::sqrt(16. - alphaZ*alphaZ);
194  G4double gamterm0 = Gamma(2.*gamma0+1.);
195  G4double gamterm1 = gamterm0/Gamma(2.*gamma1+1.);
196  G4double gamterm2 = gamterm0/Gamma(2.*gamma2+1.);
197  G4double gamterm3 = gamterm0/Gamma(2.*gamma3+1.);
198 
199  G4double term1 = e_nu*e_nu*e_nu*e_nu*e_nu*e_nu*(1. + gamma0)/1260.;
200 
201  G4double term2 = 2.*(2. + gamma1)*e_nu*e_nu*e_nu*e_nu*p_e*p_e
202  *std::pow(twoPR, 2.*(gamma1-gamma0-1.) )
203  *gamterm1*gamterm1
204  *ModSquared(gamma1, eta)/ModSquared(gamma0, eta)/5.;
205 
206  G4double term3 = 60.*(3.+gamma2)*p_e*p_e*p_e*p_e*e_nu*e_nu
207  *std::pow(twoPR, 2.*(gamma2-gamma0-2.) )
208  *gamterm2*gamterm2
209  *ModSquared(gamma2, eta)/ModSquared(gamma0, eta);
210 
211  G4double term4 = 2240.*p_e*p_e*p_e*p_e*p_e*p_e*(4. + gamma3)
212  *std::pow(twoPR, 2.*(gamma3-gamma0-3.) )
213  *gamterm3*gamterm3
214  *ModSquared(gamma3, eta)/ModSquared(gamma0, eta);
215 
216  factor = term1 + term2 + term3 + term4;
217  }
218  break;
219 
220  default:
221  G4Exception("G4BetaDecayCorrections::ShapeFactor()","HAD_RDM_010",
222  JustWarning,
223  "Transition not yet implemented - using allowed shape");
224  break;
225  }
226  return factor;
227 }
228 
229 
Generated on Sat May 25 2013 14:34:07 for Geant4 by   doxygen 1.8.4