G4ScreenedCoulombClassicalKinematics Class Reference

#include <G4ScreenedNuclearRecoil.hh>

Inheritance diagram for G4ScreenedCoulombClassicalKinematics:

Collaboration diagram for G4ScreenedCoulombClassicalKinematics:

Public Member Functions
	G4ScreenedCoulombClassicalKinematics ()
virtual void	DoCollisionStep (class G4ScreenedNuclearRecoil *master, const class G4Track &aTrack, const class G4Step &aStep)
G4bool	DoScreeningComputation (class G4ScreenedNuclearRecoil *master, const G4ScreeningTables *screen, G4double eps, G4double beta)
virtual	~G4ScreenedCoulombClassicalKinematics ()
Public Member Functions inherited from G4ScreenedCoulombCrossSectionInfo
	G4ScreenedCoulombCrossSectionInfo ()
	~G4ScreenedCoulombCrossSectionInfo ()
Public Member Functions inherited from G4ScreenedCollisionStage
virtual	~G4ScreenedCollisionStage ()

Protected Attributes
c2_const_plugin_function_p < G4double > &	phifunc
c2_linear_p< G4double > &	xovereps
G4_c2_ptr	diff

Additional Inherited Members
Static Public Member Functions inherited from G4ScreenedCoulombCrossSectionInfo
static const char *	CVSHeaderVers ()
static const char *	CVSFileVers ()

Detailed Description

Definition at line 173 of file G4ScreenedNuclearRecoil.hh.

Constructor & Destructor Documentation

G4ScreenedCoulombClassicalKinematics::G4ScreenedCoulombClassicalKinematics

(

)

Definition at line 579 of file G4ScreenedNuclearRecoil.cc.

                                                                           :
// instantiate all the needed functions statically, so no allocation is 
// done at run time
// we will be solving x^2 - x phi(x*au)/eps - beta^2 == 0.0
// or, for easier scaling, x'^2 - x' au phi(x')/eps - beta^2 au^2
// note that only the last of these gets deleted, since it owns the rest
  phifunc(c2.const_plugin_function()),
  xovereps(c2.linear(0., 0., 0.)), 
  // will fill this in with the right slope at run time
  diff(c2.quadratic(0., 0., 0., 1.)-xovereps*phifunc)
{
}

virtual G4ScreenedCoulombClassicalKinematics::~G4ScreenedCoulombClassicalKinematics ( )

inlinevirtual

Definition at line 184 of file G4ScreenedNuclearRecoil.hh.

{ }

Member Function Documentation

void G4ScreenedCoulombClassicalKinematics::DoCollisionStep ( class G4ScreenedNuclearRecoil *  master, const class G4Track &  aTrack, const class G4Step &  aStep  )

virtual

Implements G4ScreenedCollisionStage.

Definition at line 692 of file G4ScreenedNuclearRecoil.cc.

                                        {
        
        if(!master->GetValidCollision()) return;
        
        G4ParticleChange &aParticleChange=master->GetParticleChange();
        G4CoulombKinematicsInfo &kin=master->GetKinematics();
        
        const G4DynamicParticle* incidentParticle = aTrack.GetDynamicParticle();
        G4ParticleDefinition *baseParticle=aTrack.GetDefinition();
        
        G4double incidentEnergy = incidentParticle->GetKineticEnergy();
        
        // this adjustment to a1 gives the right results for soft 
        // (constant gamma)
        // relativistic collisions.  Hard collisions are wrong anyway, since the
        // Coulombic and hadronic terms interfere and cannot be added.
        G4double gamma=(1.0+incidentEnergy/baseParticle->GetPDGMass());
        G4double a1=kin.a1*gamma; // relativistic gamma correction
        
        G4ParticleDefinition *recoilIon=kin.recoilIon;
        G4double A=recoilIon->GetPDGMass()/amu_c2;
        G4int Z=(G4int)((recoilIon->GetPDGCharge()/eplus)+0.5);
        
        G4double Ec = incidentEnergy*(A/(A+a1)); 
        // energy in CM frame (non-relativistic!)
        const G4ScreeningTables *screen=kin.crossSection->GetScreening(Z);
        G4double au=screen->au; // screening length
        
        G4double beta = kin.impactParameter/au; 
        // dimensionless impact parameter
        G4double eps = Ec/(screen->z1*Z*elm_coupling/au); 
        // dimensionless energy
        
        G4bool ok=DoScreeningComputation(master, screen, eps, beta);        
        if(!ok) {
                master->SetValidCollision(0); // flag bad collision
                return; // just bail out without setting valid flag
        }
        
        G4double eRecoil=4*incidentEnergy*a1*A*kin.cosZeta*kin.cosZeta
          /((a1+A)*(a1+A));
        kin.eRecoil=eRecoil;
        
        if(incidentEnergy-eRecoil < master->GetRecoilCutoff()*a1) {
                aParticleChange.ProposeEnergy(0.0);
                master->DepositEnergy(int(screen->z1), a1, kin.targetMaterial, 
                                      incidentEnergy-eRecoil);
        } 
        
        if(master->GetEnableRecoils() && 
           eRecoil > master->GetRecoilCutoff() * kin.a2) {
                kin.recoilIon=recoilIon;                
        } else {
                kin.recoilIon=0; // this flags no recoil to be generated
                master->DepositEnergy(Z, A, kin.targetMaterial, eRecoil) ;
        }        
}

Here is the call graph for this function:

G4bool G4ScreenedCoulombClassicalKinematics::DoScreeningComputation	(	class G4ScreenedNuclearRecoil *	master,
		const G4ScreeningTables *	screen,
		G4double	eps,
		G4double	beta
	)

Definition at line 592 of file G4ScreenedNuclearRecoil.cc.

{
  G4double au=screen->au;
  G4CoulombKinematicsInfo &kin=master->GetKinematics();
  G4double A=kin.a2;
  G4double a1=kin.a1;
        
  G4double xx0; // first estimate of closest approach
  if(eps < 5.0) {
    G4double y=std::log(eps);
    G4double mlrho4=((((3.517e-4*y+1.401e-2)*y+2.393e-1)*y+2.734)*y+2.220);
    G4double rho4=std::exp(-mlrho4); // W&M eq. 18
    G4double bb2=0.5*beta*beta;
    xx0=std::sqrt(bb2+std::sqrt(bb2*bb2+rho4)); // W&M eq. 17
  } else {
    G4double ee=1.0/(2.0*eps);
    xx0=ee+std::sqrt(ee*ee+beta*beta); // W&M eq. 15 (Rutherford value)  
    if(master->CheckNuclearCollision(A, a1, xx0*au)) return 0; 
    // nuclei too close
                
  }
        
  // we will be solving x^2 - x phi(x*au)/eps - beta^2 == 0.0
  // or, for easier scaling, x'^2 - x' au phi(x')/eps - beta^2 au^2
  xovereps.reset(0., 0.0, au/eps); // slope of x*au/eps term
  phifunc.set_function(&(screen->EMphiData.get())); 
  // install interpolating table
  G4double xx1, phip, phip2;
  G4int root_error;        
  xx1=diff->find_root(phifunc.xmin(), std::min(10*xx0*au,phifunc.xmax()), 
                      std::min(xx0*au, phifunc.xmax()), beta*beta*au*au, 
                      &root_error, &phip, &phip2)/au; 
        
        if(root_error) {
                G4cout << "Screened Coulomb Root Finder Error" << G4endl;
                G4cout << "au " << au << " A " << A << " a1 " << a1 
                       << " xx1 " << xx1 << " eps " << eps 
                       << " beta " << beta << G4endl;
                G4cout << " xmin " << phifunc.xmin() << " xmax " 
                       << std::min(10*xx0*au,phifunc.xmax()) ;
                G4cout << " f(xmin) " << phifunc(phifunc.xmin()) 
                       <<  " f(xmax) " 
                       << phifunc(std::min(10*xx0*au,phifunc.xmax())) ;
                G4cout << " xstart " << std::min(xx0*au, phifunc.xmax()) 
                       <<  " target " <<  beta*beta*au*au ;
                G4cout << G4endl;
                throw c2_exception("Failed root find");
        }
        
        // phiprime is scaled by one factor of au because phi is evaluated 
        // at (xx0*au),
        G4double phiprime=phip*au;
        
        //lambda0 is from W&M 19
        G4double lambda0=1.0/std::sqrt(0.5+beta*beta/(2.0*xx1*xx1)
                                       -phiprime/(2.0*eps));
        
        // compute the 6-term Lobatto integral alpha (per W&M 21, with 
        // different coefficients)
        // this is probably completely un-needed but gives the highest 
        // quality results,
        G4double alpha=(1.0+ lambda0)/30.0;
        G4double xvals[]={0.98302349, 0.84652241, 0.53235309, 0.18347974};
        G4double weights[]={0.03472124, 0.14769029, 0.23485003, 0.18602489};
        for(G4int k=0; k<4; k++) {
                G4double x, ff;
                x=xx1/xvals[k];
                ff=1.0/std::sqrt(1.0-phifunc(x*au)/(x*eps)-beta*beta/(x*x));
                alpha+=weights[k]*ff;
        }
        
        phifunc.unset_function(); 
        // throws an exception if used without setting again

        G4double thetac1=pi*beta*alpha/xx1; 
        // complement of CM scattering angle
        G4double sintheta=std::sin(thetac1); //note sin(pi-theta)=sin(theta)
        G4double costheta=-std::cos(thetac1); // note cos(pi-theta)=-cos(theta)
        // G4double psi=std::atan2(sintheta, costheta+a1/A); 
        // lab scattering angle (M&T 3rd eq. 8.69)
        
        // numerics note:  because we checked above for reasonable values 
        // of beta which give real recoils,
        // we don't have to look too closely for theta -> 0 here 
        // (which would cause sin(theta)
        // and 1-cos(theta) to both vanish and make the atan2 ill behaved).
        G4double zeta=std::atan2(sintheta, 1-costheta); 
        // lab recoil angle (M&T 3rd eq. 8.73)
        G4double coszeta=std::cos(zeta);
        G4double sinzeta=std::sin(zeta);

        kin.sinTheta=sintheta;
        kin.cosTheta=costheta;
        kin.sinZeta=sinzeta;
        kin.cosZeta=coszeta;
        return 1; // all OK, collision is valid
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

G4_c2_ptr G4ScreenedCoulombClassicalKinematics::diff

protected

Definition at line 190 of file G4ScreenedNuclearRecoil.hh.

c2_const_plugin_function_p<G4double>& G4ScreenedCoulombClassicalKinematics::phifunc

protected

Definition at line 188 of file G4ScreenedNuclearRecoil.hh.

c2_linear_p<G4double>& G4ScreenedCoulombClassicalKinematics::xovereps

protected

Definition at line 189 of file G4ScreenedNuclearRecoil.hh.

---

The documentation for this class was generated from the following files:

• source/geant4.10.03.p03/examples/extended/electromagnetic/TestEm7/include/G4ScreenedNuclearRecoil.hh
• source/geant4.10.03.p03/examples/extended/electromagnetic/TestEm7/src/G4ScreenedNuclearRecoil.cc