G4BetaDecayCorrections.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
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   // alphaZ = fine_structure_const*std::abs(Z);
00036   alphaZ = fine_structure_const*Z;
00037 
00038   // Nuclear radius in units of hbar/m_e/c
00039   Rnuc = 0.5*fine_structure_const*std::pow(A, 0.33333);
00040 
00041   // Electron screening potential in units of electrom mass
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   // Coefficients for gamma function with real argument
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   // Calculate the relativistic Fermi function.  Argument W is the
00060   // total electron energy in units of electron mass.
00061 
00062   G4double Wprime;
00063   if (Z < 0) {
00064     Wprime = W + V0;
00065   } else {
00066     Wprime = W - V0;
00067 //    if (Wprime < 1.) Wprime = W;
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   // Fermi function
00078   G4double factor1 = 2*(1+gamma0)*mod2Gamma/realGamma/realGamma;
00079   G4double factor2 = epieta*std::pow(2*p_e*Rnuc, 2*(gamma0-1) );
00080 
00081   // Electron screening factor
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   // Calculate the squared modulus of the Gamma function 
00092   // with complex argument (re, im) using approximation B 
00093   // of Wilkinson, Nucl. Instr. & Meth. 82, 122 (1970).
00094   // Here, choose N = 1 in Wilkinson's notation for approximation B
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   // Use recursion relation to get argument < 1
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   // Calculation of Gamma function with real argument
00117   // 0 < arg < 1 using polynomial from Abramowitz and Stegun
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         // Parameters for 1st forbidden shape determined from 210Bi data
00140         // Not valid for other 1st forbidden nuclei
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  

Generated on Mon May 27 17:47:44 2013 for Geant4 by  doxygen 1.4.7