G4ScreeningMottCrossSection.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//      G4ScreeningMottCrossSection.cc
//-------------------------------------------------------------------
//
// GEANT4 Class header file
//
// File name:    G4ScreeningMottCrossSection
//
// Author:      Cristina Consolandi
//
// Creation date: 20.10.2011  
//
// Modifications:
// 27-05-2012 Added Analytic Fitting to the Mott Cross Section by means of G4MottCoefficients class.
//
//
// Class Description:
//      Computation of electron Coulomb Scattering Cross Section.
//      Suitable for high energy electrons and light target materials. 
//
//      Reference:
//      M.J. Boschini et al.
//     "Non Ionizing Energy Loss induced by Electrons in the Space Environment"
//      Proc. of the 13th International Conference on Particle Physics and Advanced Technology 
//      (13th ICPPAT, Como 3-7/10/2011), World Scientific (Singapore).
//      Available at: http://arxiv.org/abs/1111.4042v4
//
//      1) Mott Differential Cross Section Approximation:  
//         For Target material up to Z=92 (U):
//         As described in http://arxiv.org/abs/1111.4042v4 
//         par. 2.1 , eq. (16)-(17)
//         Else (Z>92):
//         W. A. McKinley and H. Fashbach, Phys. Rev. 74, (1948) 1759.
//      2) Screening coefficient: 
//      vomn G. Moliere, Z. Naturforsh A2 (1947), 133-145; A3 (1948), 78.
//      3) Nuclear Form Factor: 
//      A.V. Butkevich et al. Nucl. Instr. and Meth. in Phys. Res. A 488 (2002), 282-294.
//
// -------------------------------------------------------------------------------------
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

#include "G4ScreeningMottCrossSection.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4MottCoefficients.hh"
#include "Randomize.hh"
#include "G4Proton.hh"
#include "G4LossTableManager.hh"
#include "G4NucleiProperties.hh"
#include "G4Element.hh"
#include "G4UnitsTable.hh"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......


using namespace std;

G4ScreeningMottCrossSection::G4ScreeningMottCrossSection():
   cosThetaMin(1.0),
   cosThetaMax(-1.0),
   alpha(fine_structure_const),
   htc2(hbarc_squared),
   e2(electron_mass_c2*classic_electr_radius) 
{

        TotalCross=0;

        fNistManager = G4NistManager::Instance();
        particle=0;

        spin = mass = mu_rel=0;
        tkinLab = momLab2 = invbetaLab2=0;
        tkin = mom2 = invbeta2=beta=gamma=0;

        Trec=targetZ = targetMass = As =0;
        etag = ecut = 0.0;

        targetA = 0;

        cosTetMinNuc=0;
        cosTetMaxNuc=0;

        for(G4int i=0 ; i<5; i++){
                for(G4int j=0; j< 6; j++){
                        coeffb[i][j]=0;
                }
        }

        mottcoeff = new G4MottCoefficients();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4ScreeningMottCrossSection::~G4ScreeningMottCrossSection()
{
  delete mottcoeff;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void G4ScreeningMottCrossSection::Initialise(const G4ParticleDefinition* p,
                                          G4double CosThetaLim)
{

        SetupParticle(p);
        tkin = targetZ = mom2 = DBL_MIN;
        ecut = etag = DBL_MAX;
        particle = p;
        cosThetaMin = CosThetaLim; 

}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void G4ScreeningMottCrossSection::SetScreeningCoefficient(){

        G4double alpha2=alpha*alpha;
        //Bohr radius 
        G4double a0=  Bohr_radius  ;//0.529e-8*cm;
        //Thomas-Fermi screening length
        G4double aU=0.88534*a0/pow(targetZ,1./3.);
        G4double twoR2=aU*aU;

        G4double factor= 1.13 + 3.76*targetZ*targetZ*invbeta2*alpha2;
        As=0.25*(htc2)/(twoR2*mom2)*factor;




//cout<<"0k .........................As  "<<As<<endl;

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4double G4ScreeningMottCrossSection::GetScreeningAngle(){

                SetScreeningCoefficient();

                //cout<<" As "<<As<<endl;
                if(As < 0.0) { As = 0.0; }
                else if(As > 1.0) { As = 1.0; }

                G4double screenangle=2.*asin(sqrt(As));

        //      cout<<"  screenangle  "<<  screenangle <<endl;
        
                if(screenangle>=pi) screenangle=pi;

        
return screenangle;

}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void G4ScreeningMottCrossSection::SetupKinematic(G4double ekin, G4double Z )
{

                //...Target
                G4int iz = G4int(Z);
                G4double A = fNistManager->GetAtomicMassAmu(iz);
                G4int ia = G4int(A);
                G4double mass2 = G4NucleiProperties::GetNuclearMass(ia, iz);

                targetZ = Z;
                targetA = fNistManager->GetAtomicMassAmu(iz);
                targetMass= mass2;

                mottcoeff->SetMottCoeff(targetZ, coeffb);

        //cout<<"......... targetA "<< targetA <<endl;
        //cout<<"......... targetMass "<< targetMass/MeV <<endl;

                // incident particle lab
                tkinLab = ekin;
                momLab2 = tkinLab*(tkinLab + 2.0*mass);
                invbetaLab2 = 1.0 +  mass*mass/momLab2;

                G4double etot = tkinLab + mass;
                G4double ptot = sqrt(momLab2);
                G4double m12  = mass*mass;
                

        // relativistic reduced mass from publucation
        // A.P. Martynenko, R.N. Faustov, Teoret. mat. Fiz. 64 (1985) 179
        
                //incident particle & target nucleus
                G4double Ecm=sqrt(m12 + mass2*mass2 + 2.0*etot*mass2);
                 mu_rel=mass*mass2/Ecm;
                G4double momCM= ptot*mass2/Ecm;
                // relative system
                mom2 = momCM*momCM;
                invbeta2 = 1.0 +  mu_rel*mu_rel/mom2;
                tkin = momCM*sqrt(invbeta2) - mu_rel;//Ekin of mu_rel
                        G4double  beta2=1./invbeta2;
                beta=std::sqrt(beta2) ;
                        G4double gamma2= 1./(1.-beta2);
                gamma=std::sqrt(gamma2);



               //.........................................................


                G4double screenangle=GetScreeningAngle()/10.;
                //cout<<" screenangle [rad] "<<screenangle/rad <<endl;

                cosTetMinNuc =min( cosThetaMin ,cos(screenangle));
                cosTetMaxNuc =cosThetaMax;
        
//cout<<"ok..............mu_rel "<<mu_rel<<endl;

}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4double G4ScreeningMottCrossSection::FormFactor2ExpHof(G4double angle){


       G4double M=targetMass; 
       G4double E=tkinLab;
       G4double Etot=E+mass;
       G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
       G4double T=Tmax*pow(sin(angle/2.),2.);
       G4double q2=T*(T+2.*M);
             q2/=htc2;//1/cm2
       G4double RN=1.27e-13*pow(targetA,0.27)*cm;
       G4double xN= (RN*RN*q2);
       G4double den=(1.+xN/12.);
       G4double FN=1./(den*den);
       G4double form2=(FN*FN);


  return form2;

//cout<<"..................... form2 "<< form2<<endl;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4double G4ScreeningMottCrossSection::McFcorrection(G4double angle ){


        G4double  beta2=1./invbeta2;
        G4double sintmezzi=std::sin(angle/2.);
        G4double sin2tmezzi = sintmezzi*sintmezzi;
        G4double R=1.-beta2*sin2tmezzi + targetZ*alpha*beta*pi*sintmezzi*(1.-sintmezzi);


return R;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double G4ScreeningMottCrossSection::RatioMottRutherford(G4double angle){


        G4double R=0;
        G4double fcost=std::sqrt((1. -cos(angle)));
        G4double a[5];
        G4double shift=0.7181228;
        G4double beta0= beta -shift;

              for(G4int j=0 ;j<=4;j++){
                a[j]=0;
                }

                for(G4int j=0 ;j<=4;j++){
                        for(G4int k=0;k<=5;k++ ){  
                                a[j]+=coeffb[j][k]*pow(beta0,k);
                        }
                }

        for(G4int j=0 ;j<=4 ;j++){
                R+=a[j]* pow(fcost,j);
                }



return R;

}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4double G4ScreeningMottCrossSection::NuclearCrossSection()
{
         if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;

        TotalCross=0;

        G4double anglemin =std::acos(cosTetMinNuc);
        G4double anglemax =std::acos(cosTetMaxNuc); 

        static const G4double limit = 1.e-9;
        if(anglemin < limit) {
          anglemin = GetScreeningAngle()/10.;
          if(anglemin < limit) { anglemin = limit; }
        }

        //cout<<" anglemin  "<< anglemin <<endl;

        G4double loganglemin=log10(anglemin);
        G4double loganglemax=log10(anglemax);
        G4double logdangle=0.01;

        G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);


        vector<G4double> angle;
        vector<G4double> tet;
        vector<G4double> dangle;
        vector<G4double> cross;

              for(G4int i=0; i<=bins; i++ ){
                     tet.push_back(0);
                     dangle.push_back(0);
                     angle.push_back(pow(10.,loganglemin+logdangle*i));
                     cross.push_back(0);
                        }


                G4int  dim = tet.size();
                //cout<<"dim--- "<<dim<<endl;


                for(G4int i=0; i<dim;i++){

                        if(i!=dim-1){
                        dangle[i]=(angle[i+1]-angle[i]);
                        tet[i]=(angle[i+1]+angle[i])/2.;
                        }else if(i==dim-1){
                                break;
                        }


                        G4double R=0;
                        G4double F2=FormFactor2ExpHof(tet[i]);
                        
                        if (coeffb[0][0]!=0){
                        //cout<<" Mott....targetZ "<< targetZ<<endl;    
                        R=RatioMottRutherford(tet[i]);
                        }

                        else if (coeffb[0][0]==0){
                       // cout<<" McF.... targetZ "<< targetZ<<endl;
                        R=McFcorrection(tet[i]);
                        }


//cout<<"----------------- R "<<R<<" F2 "<<F2<<endl;


//                cout<<"angle "<<tet[i] << " F2 "<<F2<<endl;

                        G4double e4=e2*e2;
                        G4double den=2.*As+2.*pow(sin(tet[i]/2.),2.);
                        G4double func=1./(den*den);

                        G4double fatt= targetZ/(mu_rel*gamma*beta*beta);
                        G4double sigma=e4*fatt*fatt*func;
                        cross[i]=F2*R*sigma;
                        G4double pi2sintet=2.*pi*sin(tet[i]);

                        TotalCross+=pi2sintet*cross[i]*dangle[i];
                        }//end integral


//cout<< "ok ......... TotalCross "<<TotalCross<<endl;


return TotalCross;

}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double G4ScreeningMottCrossSection::AngleDistribution(G4double anglein){


        G4double total=TotalCross ;
        G4double fatt= e2*targetZ/(mu_rel*gamma*beta*beta);
        G4double fatt2=fatt*fatt;
        total/=fatt2;

        G4double R=0;
                 if (coeffb[0][0]!=0){
                     //   cout<<" Mott....targetZ "<< targetZ<<endl;      
                        R=RatioMottRutherford(anglein);
                        }

                 else if (coeffb[0][0]==0){
                       // cout<<" McF.... targetZ "<< targetZ<<endl;
                        R=McFcorrection(anglein);
                        }


        G4double y=2.*pi*sin(anglein)*R*FormFactor2ExpHof(anglein)/
                                ((2*As+2.*pow(sin(anglein/2.),2.))*(2.*As+2.*pow(sin(anglein/2.),2.) ));

return y/total;

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4double G4ScreeningMottCrossSection::GetScatteringAngle()
{
  

//cout<<" G4ScreeningMottCrossSection::SampleCosineTheta ............."<<endl;  

        if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;

        G4double anglemin=std::acos(cosTetMinNuc);      
        G4double anglemax= std::acos(cosTetMaxNuc);

        static const G4double limit = 1.e-9;
        if(anglemin < limit) {
          anglemin = GetScreeningAngle()/10.;
          if(anglemin < limit) { anglemin = limit; }
        }

//      cout<<"................ tkinLab  "<< G4BestUnit(tkinLab,"Energy") << " anglemin=  "<<anglemin<<endl;
        //cout<<"anglemax=  "<<anglemax<<endl;
        G4double r =G4UniformRand();

        G4double loganglemin=log10(anglemin);
        G4double loganglemax=log10(anglemax);
        G4double logdangle=0.01;

        G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);

        std::vector<G4double> angle;
        std::vector<G4double> tet;
        std::vector<G4double> dangle;

              for(G4int i=0; i<=bins; i++ ){
                     tet.push_back(0);
                     dangle.push_back(0);
                     angle.push_back(pow(10.,loganglemin+logdangle*i));
                        }


                G4int  dim = tet.size();
                G4double scattangle=0;
                G4double y=0;
                G4double dy=0;
                G4double area=0;

                for(G4int i=0; i<dim;i++){

                        if(i!=dim-1){
                        dangle[i]=(angle[i+1]-angle[i]);
                        tet[i]=(angle[i+1]+angle[i])/2.;
                        }else if(i==dim-1){
                                break;
                        }

                        y+=AngleDistribution(tet[i])*dangle[i];
                        dy= y-area ;
                        area=y;

                        if(r >=y-dy && r<=y+dy ){       
                        scattangle= angle[i] +G4UniformRand()*dangle[i];
                        //cout<<"y "<<y <<" r  "<< r  << " y/r "<<y/r << endl;
                        break;

                        }       
                
                }


  return scattangle;


}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4ThreeVector G4ScreeningMottCrossSection::GetNewDirection(){

G4ThreeVector dir(0.0,0.0,1.0);

        
        G4double z1=GetScatteringAngle();

        G4double sint = sin(z1);
        G4double cost = sqrt(1.0 - sint*sint);
        G4double phi  = twopi* G4UniformRand();
        G4double dirx = sint*cos(phi);
        G4double diry = sint*sin(phi);
        G4double dirz = cost;


        //.......set Trc
        G4double etot=tkinLab+mass;
        G4double mass2=targetMass;
        Trec=(1.0 - cost)* mass2*(etot*etot - mass*mass )/
                              (mass*mass + mass2*mass2+ 2.*mass2*etot);
        
        dir.set(dirx,diry,dirz);

return dir;

}