G4ScreenedNuclearRecoil.cc File Reference

Implementation of the G4ScreenedNuclearRecoil class. More...

#include <stdio.h>
#include "globals.hh"
#include "G4ScreenedNuclearRecoil.hh"
#include "G4ParticleTypes.hh"
#include "G4ParticleTable.hh"
#include "G4IonTable.hh"
#include "G4VParticleChange.hh"
#include "G4ParticleChangeForLoss.hh"
#include "G4DataVector.hh"
#include "G4Track.hh"
#include "G4Step.hh"
#include "G4Material.hh"
#include "G4Element.hh"
#include "G4Isotope.hh"
#include "G4MaterialCutsCouple.hh"
#include "G4ElementVector.hh"
#include "G4IsotopeVector.hh"
#include "G4EmProcessSubType.hh"
#include "G4ParticleDefinition.hh"
#include "G4DynamicParticle.hh"
#include "G4ProcessManager.hh"
#include "G4StableIsotopes.hh"
#include "G4LindhardPartition.hh"
#include "Randomize.hh"
#include <iostream>
#include <iomanip>
#include "c2_factory.hh"
#include <vector>

Include dependency graph for G4ScreenedNuclearRecoil.cc:

Go to the source code of this file.

Typedefs
typedef c2_ptr< G4double >	c2p

Functions
G4_c2_function &	ZBLScreening (G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
G4_c2_function &	MoliereScreening (G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
G4_c2_function &	LJScreening (G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
G4_c2_function &	LJZBLScreening (G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval)
static G4double	cm_energy (G4double a1, G4double a2, G4double t0)
static G4double	thetac (G4double m1, G4double mass2, G4double eratio)

Variables
static c2_factory< G4double >	c2

Detailed Description

Implementation of the G4ScreenedNuclearRecoil class.

Definition in file G4ScreenedNuclearRecoil.cc.

Typedef Documentation

typedef c2_ptr<G4double> c2p

Definition at line 127 of file G4ScreenedNuclearRecoil.cc.

Function Documentation

static G4double cm_energy ( G4double  a1, G4double  a2, G4double  t0  )

inlinestatic

Definition at line 979 of file G4ScreenedNuclearRecoil.cc.

                                                                        {
        // "relativistically correct energy in CM frame"
        G4double m1=a1*amu_c2, mass2=a2*amu_c2;
        G4double mc2=(m1+mass2);
        G4double f=2.0*mass2*t0/(mc2*mc2);
        // old way: return (f < 1e-6) ?  0.5*mc2*f : mc2*(std::sqrt(1.0+f)-1.0);
        // formally equivalent to previous, but numerically stable for all 
        // f without conditional
        // uses identity (sqrt(1+x) - 1)(sqrt(1+x) + 1) = x
        return mc2*f/(std::sqrt(1.0+f)+1.0); 
}

Here is the caller graph for this function:

G4_c2_function& LJScreening	(	G4int	z1,
		G4int	z2,
		size_t	npoints,
		G4double	rMax,
		G4double *	auval
	)

Definition at line 899 of file G4ScreenedNuclearRecoil.cc.

{        
//from Loftager, Besenbacher, Jensen & Sorensen
//PhysRev A20, 1443++, 1979
        G4double au=0.8853*0.529*angstrom/std::sqrt(std::pow(z1, 0.6667)
                                                    +std::pow(z2,0.6667));
        std::vector<G4double> r(npoints), phi(npoints);
        
        for(size_t i=0; i<npoints; i++) {
                G4double rr=(float)i/(float)(npoints-1);
                r[i]=rr*rr*rMax; 
                // use quadratic r scale to make sampling fine near the center

                G4double y=std::sqrt(9.67*r[i]/au);
                G4double ysq=y*y;
                G4double phipoly=1+y+0.3344*ysq+0.0485*y*ysq+0.002647*ysq*ysq;
                phi[i]=phipoly*std::exp(-y);
                // G4cout << r[i] << " " << phi[i] << G4endl;
        }

        // compute the derivative at the origin for the spline
        G4double logphiprime0=(9.67/2.0)*(2*0.3344-1.0); 
        // #avoid 0/0 on first element
        logphiprime0 *= (1.0/au); // #put back in natural units
        
        *auval=au;
        return c2.lin_log_interpolating_function().load(r, phi, false, 
                                                        logphiprime0*phi[0],
                                                        true,0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

G4_c2_function& LJZBLScreening	(	G4int	z1,
		G4int	z2,
		size_t	npoints,
		G4double	rMax,
		G4double *	auval
	)

connector in between. These numbers are selected so the switchover

Definition at line 931 of file G4ScreenedNuclearRecoil.cc.

{        
// hybrid of LJ and ZBL, uses LJ if x < 0.25*auniv, ZBL if x > 1.5*auniv, and 
// is very near the point where the functions naturally cross.
        G4double auzbl, aulj;
        
        c2p zbl=ZBLScreening(z1, z2, npoints, rMax, &auzbl);
        c2p lj=LJScreening(z1, z2, npoints, rMax, &aulj);

        G4double au=(auzbl+aulj)*0.5;
        lj->set_domain(lj->xmin(), 0.25*au);
        zbl->set_domain(1.5*au,zbl->xmax());
        
        c2p conn=
          c2.connector_function(lj->xmax(), lj, zbl->xmin(), zbl, true,0);
        c2_piecewise_function_p<G4double> &pw=c2.piecewise_function();
        c2p keepit(pw);
        pw.append_function(lj);
        pw.append_function(conn);
        pw.append_function(zbl);
        
        *auval=au;
        keepit.release_for_return();
        return pw;
}

Here is the call graph for this function:

Here is the caller graph for this function:

G4_c2_function& MoliereScreening	(	G4int	z1,
		G4int	z2,
		size_t	npoints,
		G4double	rMax,
		G4double *	auval
	)

Definition at line 867 of file G4ScreenedNuclearRecoil.cc.

{        
        static const size_t ncoef=3;
        static G4double scales[ncoef]={-6.0, -1.2, -0.3};
        static G4double coefs[ncoef]={0.10, 0.55, 0.35};
        
        G4double au=0.8853*0.529*angstrom/std::sqrt(std::pow(z1, 0.6667)
                                                    +std::pow(z2,0.6667));
        std::vector<G4double> r(npoints), phi(npoints);
        
        for(size_t i=0; i<npoints; i++) {
                G4double rr=(float)i/(float)(npoints-1);
                r[i]=rr*rr*rMax; 
                // use quadratic r scale to make sampling fine near the center
                G4double sum=0.0;
                for(size_t j=0; j<ncoef; j++) 
                  sum+=coefs[j]*std::exp(scales[j]*r[i]/au);
                phi[i]=sum;
        }
        
        // compute the derivative at the origin for the spline
        G4double phiprime0=0.0;
        for(size_t j=0; j<ncoef; j++) 
          phiprime0+=scales[j]*coefs[j]*std::exp(scales[j]*r[0]/au);
        phiprime0*=(1.0/au); // put back in natural units;
        
        *auval=au;
        return c2.lin_log_interpolating_function().load(r, phi, false, 
                                                        phiprime0,true,0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

static G4double thetac ( G4double  m1, G4double  mass2, G4double  eratio  )

inlinestatic

Definition at line 991 of file G4ScreenedNuclearRecoil.cc.

                                                                             {
        G4double s2th2=eratio*( (m1+mass2)*(m1+mass2)/(4.0*m1*mass2) );
        G4double sth2=std::sqrt(s2th2);
        return 2.0*std::asin(sth2);
}

Here is the caller graph for this function:

G4_c2_function& ZBLScreening	(	G4int	z1,
		G4int	z2,
		size_t	npoints,
		G4double	rMax,
		G4double *	auval
	)

Definition at line 835 of file G4ScreenedNuclearRecoil.cc.

{
        static const size_t ncoef=4;
        static G4double scales[ncoef]={-3.2, -0.9432, -0.4028, -0.2016};
        static G4double coefs[ncoef]={0.1818,0.5099,0.2802,0.0281};
        
        G4double au=
          0.8854*angstrom*0.529/(std::pow(z1, 0.23)+std::pow(z2,0.23));
        std::vector<G4double> r(npoints), phi(npoints);
        
        for(size_t i=0; i<npoints; i++) {
                G4double rr=(float)i/(float)(npoints-1);
                r[i]=rr*rr*rMax; 
                // use quadratic r scale to make sampling fine near the center
                G4double sum=0.0;
                for(size_t j=0; j<ncoef; j++) 
                  sum+=coefs[j]*std::exp(scales[j]*r[i]/au);
                phi[i]=sum;
        }

        // compute the derivative at the origin for the spline
        G4double phiprime0=0.0;
        for(size_t j=0; j<ncoef; j++) 
          phiprime0+=scales[j]*coefs[j]*std::exp(scales[j]*r[0]/au);
        phiprime0*=(1.0/au); // put back in natural units;
        
        *auval=au;
        return c2.lin_log_interpolating_function().load(r, phi, false, 
                                                        phiprime0,true,0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

c2_factory<G4double> c2

static

Definition at line 126 of file G4ScreenedNuclearRecoil.cc.