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;
}