00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include "globals.hh"
00028 #include "G4PhysicalConstants.hh"
00029 #include "G4BetaDecayType.hh"
00030 #include "G4BetaDecayCorrections.hh"
00031
00032 G4BetaDecayCorrections::G4BetaDecayCorrections(G4int theZ, G4int theA)
00033 : Z(theZ), A(theA)
00034 {
00035
00036 alphaZ = fine_structure_const*Z;
00037
00038
00039 Rnuc = 0.5*fine_structure_const*std::pow(A, 0.33333);
00040
00041
00042 V0 = 1.13*fine_structure_const*fine_structure_const
00043 *std::pow(std::abs(Z), 1.33333);
00044
00045 gamma0 = std::sqrt(1. - alphaZ*alphaZ);
00046
00047
00048 gc[0] = -0.1010678;
00049 gc[1] = 0.4245549;
00050 gc[2] = -0.6998588;
00051 gc[3] = 0.9512363;
00052 gc[4] = -0.5748646;
00053 gc[5] = 1.0;
00054 }
00055
00056
00057 G4double G4BetaDecayCorrections::FermiFunction(const G4double& W)
00058 {
00059
00060
00061
00062 G4double Wprime;
00063 if (Z < 0) {
00064 Wprime = W + V0;
00065 } else {
00066 Wprime = W - V0;
00067
00068 if (Wprime <= 1.00001) Wprime = 1.00001;
00069 }
00070
00071 G4double p_e = std::sqrt(Wprime*Wprime - 1.);
00072 G4double eta = alphaZ*Wprime/p_e;
00073 G4double epieta = std::exp(pi*eta);
00074 G4double realGamma = Gamma(2.*gamma0+1);
00075 G4double mod2Gamma = ModSquared(gamma0, eta);
00076
00077
00078 G4double factor1 = 2*(1+gamma0)*mod2Gamma/realGamma/realGamma;
00079 G4double factor2 = epieta*std::pow(2*p_e*Rnuc, 2*(gamma0-1) );
00080
00081
00082 G4double factor3 = (Wprime/W)*std::sqrt( (Wprime*Wprime - 1.)/(W*W - 1.) );
00083
00084 return factor1*factor2*factor3;
00085 }
00086
00087
00088 G4double
00089 G4BetaDecayCorrections::ModSquared(const G4double& re, const G4double& im)
00090 {
00091
00092
00093
00094
00095
00096 G4double factor1 = std::pow( (1+re)*(1+re) + im*im, re+0.5);
00097 G4double factor2 = std::exp(2*im * std::atan(im/(1+re)));
00098 G4double factor3 = std::exp(2*(1+re));
00099 G4double factor4 = 2.*pi;
00100 G4double factor5 = std::exp( (1+re)/( (1+re)*(1+re) + im*im)/6 );
00101 G4double factor6 = re*re + im*im;
00102 return factor1*factor4*factor5/factor2/factor3/factor6;
00103 }
00104
00105
00106 G4double G4BetaDecayCorrections::Gamma(const G4double& arg)
00107 {
00108
00109 G4double fac = 1.0;
00110 G4double x = arg - 1.;
00111 while (x > 1.0) {
00112 fac *= x;
00113 x -= 1.0;
00114 }
00115
00116
00117
00118 G4double sum = gc[0];
00119 for (G4int i = 1; i < 6; i++) sum = sum*x + gc[i];
00120
00121 return sum*fac;
00122 }
00123
00124
00125 G4double
00126 G4BetaDecayCorrections::ShapeFactor(const G4BetaDecayType& bdt,
00127 const G4double& p_e, const G4double& e_nu)
00128 {
00129 G4double twoPR = 2.*p_e*Rnuc;
00130 G4double factor(1.);
00131
00132 switch (bdt)
00133 {
00134 case (allowed) :
00135 break;
00136
00137 case (firstForbidden) :
00138 {
00139
00140
00141 G4double c1 = 0.578;
00142 G4double c2 = 28.466;
00143 G4double c3 = -0.658;
00144
00145 G4double w = std::sqrt(1. + p_e*p_e);
00146 factor = 1. + c1*w + c2/w + c3*w*w;
00147 }
00148 break;
00149
00150 case (uniqueFirstForbidden) :
00151 {
00152 G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
00153 G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
00154 G4double gamterm1 = Gamma(2.*gamma0+1.)/Gamma(2.*gamma1+1.);
00155 G4double term1 = e_nu*e_nu*(1. + gamma0)/6.;
00156 G4double term2 = 12.*(2. + gamma1)*p_e*p_e
00157 *std::pow(twoPR, 2.*(gamma1-gamma0-1) )
00158 *gamterm1*gamterm1
00159 *ModSquared(gamma1, eta)/ModSquared(gamma0, eta);
00160 factor = term1 + term2;
00161 }
00162 break;
00163
00164 case (uniqueSecondForbidden) :
00165 {
00166 G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
00167 G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
00168 G4double gamma2 = std::sqrt(9. - alphaZ*alphaZ);
00169 G4double gamterm0 = Gamma(2.*gamma0+1.);
00170 G4double gamterm1 = gamterm0/Gamma(2.*gamma1+1.);
00171 G4double gamterm2 = gamterm0/Gamma(2.*gamma2+1.);
00172 G4double term1 = e_nu*e_nu*e_nu*e_nu*(1. + gamma0)/60.;
00173
00174 G4double term2 = 4.*(2. + gamma1)*e_nu*e_nu*p_e*p_e
00175 *std::pow(twoPR, 2.*(gamma1-gamma0-1.) )
00176 *gamterm1*gamterm1
00177 *ModSquared(gamma1, eta)/ModSquared(gamma0, eta);
00178
00179 G4double term3 = 180.*(3.+gamma2)*p_e*p_e*p_e*p_e
00180 *std::pow(twoPR, 2.*(gamma2-gamma0-2) )
00181 *gamterm2*gamterm2
00182 *ModSquared(gamma2, eta)/ModSquared(gamma0, eta);
00183
00184 factor = term1 + term2 + term3;
00185 }
00186 break;
00187
00188 case (uniqueThirdForbidden) :
00189 {
00190 G4double eta = alphaZ*std::sqrt(1. + p_e*p_e)/p_e;
00191 G4double gamma1 = std::sqrt(4. - alphaZ*alphaZ);
00192 G4double gamma2 = std::sqrt(9. - alphaZ*alphaZ);
00193 G4double gamma3 = std::sqrt(16. - alphaZ*alphaZ);
00194 G4double gamterm0 = Gamma(2.*gamma0+1.);
00195 G4double gamterm1 = gamterm0/Gamma(2.*gamma1+1.);
00196 G4double gamterm2 = gamterm0/Gamma(2.*gamma2+1.);
00197 G4double gamterm3 = gamterm0/Gamma(2.*gamma3+1.);
00198
00199 G4double term1 = e_nu*e_nu*e_nu*e_nu*e_nu*e_nu*(1. + gamma0)/1260.;
00200
00201 G4double term2 = 2.*(2. + gamma1)*e_nu*e_nu*e_nu*e_nu*p_e*p_e
00202 *std::pow(twoPR, 2.*(gamma1-gamma0-1.) )
00203 *gamterm1*gamterm1
00204 *ModSquared(gamma1, eta)/ModSquared(gamma0, eta)/5.;
00205
00206 G4double term3 = 60.*(3.+gamma2)*p_e*p_e*p_e*p_e*e_nu*e_nu
00207 *std::pow(twoPR, 2.*(gamma2-gamma0-2.) )
00208 *gamterm2*gamterm2
00209 *ModSquared(gamma2, eta)/ModSquared(gamma0, eta);
00210
00211 G4double term4 = 2240.*p_e*p_e*p_e*p_e*p_e*p_e*(4. + gamma3)
00212 *std::pow(twoPR, 2.*(gamma3-gamma0-3.) )
00213 *gamterm3*gamterm3
00214 *ModSquared(gamma3, eta)/ModSquared(gamma0, eta);
00215
00216 factor = term1 + term2 + term3 + term4;
00217 }
00218 break;
00219
00220 default:
00221 G4Exception("G4BetaDecayCorrections::ShapeFactor()","HAD_RDM_010",
00222 JustWarning,
00223 "Transition not yet implemented - using allowed shape");
00224 break;
00225 }
00226 return factor;
00227 }
00228
00229