#include <G4ScreenedNuclearRecoil.hh>

Inheritance 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 159 of file G4ScreenedNuclearRecoil.hh.

Constructor & Destructor Documentation

G4ScreenedCoulombClassicalKinematics::G4ScreenedCoulombClassicalKinematics ( )

Definition at line 518 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 169 of file G4ScreenedNuclearRecoil.hh.

169 { }

Member Function Documentation

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

virtual

Implements G4ScreenedCollisionStage.

Definition at line 609 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) ;
         }        
 }

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

Definition at line 529 of file G4ScreenedNuclearRecoil.cc.

References G4CoulombKinematicsInfo::a1, G4CoulombKinematicsInfo::a2, G4ScreeningTables::au, G4ScreenedNuclearRecoil::CheckNuclearCollision(), G4CoulombKinematicsInfo::cosTheta, G4CoulombKinematicsInfo::cosZeta, diff, G4ScreeningTables::EMphiData, c2_function< float_type >::find_root(), G4cout, G4endl, c2_const_ptr< float_type >::get(), G4ScreenedNuclearRecoil::GetKinematics(), G4INCL::Math::min(), phifunc, c2_linear_p< float_type >::reset(), c2_const_plugin_function_p< float_type >::set_function(), G4CoulombKinematicsInfo::sinTheta, G4CoulombKinematicsInfo::sinZeta, c2_plugin_function_p< float_type >::unset_function(), test::x, c2_function< float_type >::xmax(), c2_function< float_type >::xmin(), and xovereps.

Referenced by DoCollisionStep().

 {
         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=CLHEP::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
 }