ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BetaDecayCorrections.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BetaDecayCorrections.cc
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(const G4int theZ, const 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 electron 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 
112  G4int loop = 0;
113  G4ExceptionDescription ed;
114  ed << " While count exceeded " << G4endl;
115  while (x > 1.0) { /* Loop checking, 01.09.2015, D.Wright */
116  fac *= x;
117  x -= 1.0;
118  loop++;
119  if (loop > 1000) {
120  G4Exception("G4BetaDecayCorrections::Gamma()", "HAD_RDM_100", JustWarning, ed);
121  break;
122  }
123  }
124 
125  // Calculation of Gamma function with real argument
126  // 0 < arg < 1 using polynomial from Abramowitz and Stegun
127  G4double sum = gc[0];
128  for (G4int i = 1; i < 6; i++) sum = sum*x + gc[i];
129 
130  return sum*fac;
131 }
132 
133 
134 G4double
135 G4BetaDecayCorrections::ShapeFactor(const G4BetaDecayType& bdt,
136  const G4double& p_e, const G4double& e_nu)
137 {
138  G4double twoPR = 2.*p_e*Rnuc;
139  G4double factor(1.);
140 
141  switch (bdt)
142  {
143  case (allowed) :
144  break;
145 
146  case (firstForbidden) :
147  {
148  // Parameters for 1st forbidden shape determined from 210Bi data
149  // Not valid for other 1st forbidden nuclei
150  G4double c1 = 0.578;
151  G4double c2 = 28.466;
152  G4double c3 = -0.658;
153 
154  G4double w = std::sqrt(1. + p_e*p_e);
155  factor = 1. + c1*w + c2/w + c3*w*w;
156  }
157  break;
158 
159  case (uniqueFirstForbidden) :
160  {
161  G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
162  G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
163  G4double gamterm1 = Gamma(2.*gamma0+1.)/Gamma(2.*gamma1+1.);
164  G4double term1 = e_nu*e_nu*(1. + gamma0)/6.;
165  G4double term2 = 12.*(2. + gamma1)*p_e*p_e
166  *std::pow(twoPR, 2.*(gamma1-gamma0-1) )
167  *gamterm1*gamterm1
168  *ModSquared(gamma1, eta)/ModSquared(gamma0, eta);
169  factor = term1 + term2;
170  }
171  break;
172 
173  case (secondForbidden) :
174  break;
175 
176  case (uniqueSecondForbidden) :
177  {
178  G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
179  G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
180  G4double gamma2 = std::sqrt(9. - alphaZ*alphaZ);
181  G4double gamterm0 = Gamma(2.*gamma0+1.);
182  G4double gamterm1 = gamterm0/Gamma(2.*gamma1+1.);
183  G4double gamterm2 = gamterm0/Gamma(2.*gamma2+1.);
184  G4double term1 = e_nu*e_nu*e_nu*e_nu*(1. + gamma0)/60.;
185 
186  G4double term2 = 4.*(2. + gamma1)*e_nu*e_nu*p_e*p_e
187  *std::pow(twoPR, 2.*(gamma1-gamma0-1.) )
188  *gamterm1*gamterm1
189  *ModSquared(gamma1, eta)/ModSquared(gamma0, eta);
190 
191  G4double term3 = 180.*(3.+gamma2)*p_e*p_e*p_e*p_e
192  *std::pow(twoPR, 2.*(gamma2-gamma0-2) )
193  *gamterm2*gamterm2
194  *ModSquared(gamma2, eta)/ModSquared(gamma0, eta);
195 
196  factor = term1 + term2 + term3;
197  }
198  break;
199 
200  case (thirdForbidden) :
201  break;
202 
203  case (uniqueThirdForbidden) :
204  {
205  G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
206  G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
207  G4double gamma2 = std::sqrt(9. - alphaZ*alphaZ);
208  G4double gamma3 = std::sqrt(16. - alphaZ*alphaZ);
209  G4double gamterm0 = Gamma(2.*gamma0+1.);
210  G4double gamterm1 = gamterm0/Gamma(2.*gamma1+1.);
211  G4double gamterm2 = gamterm0/Gamma(2.*gamma2+1.);
212  G4double gamterm3 = gamterm0/Gamma(2.*gamma3+1.);
213 
214  G4double term1 = e_nu*e_nu*e_nu*e_nu*e_nu*e_nu*(1. + gamma0)/1260.;
215 
216  G4double term2 = 2.*(2. + gamma1)*e_nu*e_nu*e_nu*e_nu*p_e*p_e
217  *std::pow(twoPR, 2.*(gamma1-gamma0-1.) )
218  *gamterm1*gamterm1
219  *ModSquared(gamma1, eta)/ModSquared(gamma0, eta)/5.;
220 
221  G4double term3 = 60.*(3.+gamma2)*p_e*p_e*p_e*p_e*e_nu*e_nu
222  *std::pow(twoPR, 2.*(gamma2-gamma0-2.) )
223  *gamterm2*gamterm2
224  *ModSquared(gamma2, eta)/ModSquared(gamma0, eta);
225 
226  G4double term4 = 2240.*p_e*p_e*p_e*p_e*p_e*p_e*(4. + gamma3)
227  *std::pow(twoPR, 2.*(gamma3-gamma0-3.) )
228  *gamterm3*gamterm3
229  *ModSquared(gamma3, eta)/ModSquared(gamma0, eta);
230 
231  factor = term1 + term2 + term3 + term4;
232  }
233  break;
234 
235  default:
236  G4Exception("G4BetaDecayCorrections::ShapeFactor()","HAD_RDM_010",
237  JustWarning,
238  "Transition not yet implemented - using allowed shape");
239  break;
240  }
241  return factor;
242 }
243 
244