G4ScreeningMottCrossSection Class Reference

#include <G4ScreeningMottCrossSection.hh>


Public Member Functions

 G4ScreeningMottCrossSection ()
virtual ~G4ScreeningMottCrossSection ()
void Initialise (const G4ParticleDefinition *, G4double cosThetaLim)
G4double GetScreeningAngle ()
void SetScreeningCoefficient ()
void SetupParticle (const G4ParticleDefinition *)
void SetupKinematic (G4double kinEnergy, G4double Z)
G4double NuclearCrossSection ()
G4ThreeVector GetNewDirection ()
G4double GetMom2CM () const
G4double GetMom2Lab () const
G4double GetTrec () const
G4double GetScreeningCoefficient () const
G4double GetTotalCross () const
G4double McFcorrection (G4double)
G4double RatioMottRutherford (G4double)
G4double FormFactor2ExpHof (G4double)
G4double GetScatteringAngle ()
G4double AngleDistribution (G4double)

Detailed Description

Definition at line 86 of file G4ScreeningMottCrossSection.hh.

Constructor & Destructor Documentation

G4ScreeningMottCrossSection::G4ScreeningMottCrossSection (  ) 

Definition at line 83 of file G4ScreeningMottCrossSection.cc.

References G4NistManager::Instance().

00083                                                         :
00084    cosThetaMin(1.0),
00085    cosThetaMax(-1.0),
00086    alpha(fine_structure_const),
00087    htc2(hbarc_squared),
00088    e2(electron_mass_c2*classic_electr_radius) 
00089 {
00090 
00091         TotalCross=0;
00092 
00093         fNistManager = G4NistManager::Instance();
00094         particle=0;
00095 
00096         spin = mass = mu_rel=0;
00097         tkinLab = momLab2 = invbetaLab2=0;
00098         tkin = mom2 = invbeta2=beta=gamma=0;
00099 
00100         Trec=targetZ = targetMass = As =0;
00101         etag = ecut = 0.0;
00102 
00103         targetA = 0;
00104 
00105         cosTetMinNuc=0;
00106         cosTetMaxNuc=0;
00107 
00108         for(G4int i=0 ; i<5; i++){
00109                 for(G4int j=0; j< 6; j++){
00110                         coeffb[i][j]=0;
00111                 }
00112         }
00113 
00114 
00115 
00116 
00117 }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4ScreeningMottCrossSection::~G4ScreeningMottCrossSection (  )  [virtual]

Definition at line 120 of file G4ScreeningMottCrossSection.cc.

00121 {}

Member Function Documentation

G4double G4ScreeningMottCrossSection::AngleDistribution ( G4double   ) 

Definition at line 393 of file G4ScreeningMottCrossSection.cc.

References FormFactor2ExpHof(), McFcorrection(), G4INCL::Math::pi, and RatioMottRutherford().

Referenced by GetScatteringAngle().

00393                                                                        {
00394 
00395 
00396         G4double total=TotalCross ;
00397         G4double fatt= e2*targetZ/(mu_rel*gamma*beta*beta);
00398         G4double fatt2=fatt*fatt;
00399         total/=fatt2;
00400 
00401         G4double R=0;
00402                  if (coeffb[0][0]!=0){
00403                      //   cout<<" Mott....targetZ "<< targetZ<<endl;      
00404                         R=RatioMottRutherford(anglein);
00405                         }
00406 
00407                  else if (coeffb[0][0]==0){
00408                        // cout<<" McF.... targetZ "<< targetZ<<endl;
00409                         R=McFcorrection(anglein);
00410                         }
00411 
00412 
00413         G4double y=2.*pi*sin(anglein)*R*FormFactor2ExpHof(anglein)/
00414                                 ((2*As+2.*pow(sin(anglein/2.),2.))*(2.*As+2.*pow(sin(anglein/2.),2.) ));
00415 
00416 return y/total;
00417 
00418 }

G4double G4ScreeningMottCrossSection::FormFactor2ExpHof ( G4double   ) 

Definition at line 237 of file G4ScreeningMottCrossSection.cc.

Referenced by AngleDistribution(), and NuclearCrossSection().

00237                                                                      {
00238 
00239 
00240        G4double M=targetMass; 
00241        G4double E=tkinLab;
00242        G4double Etot=E+mass;
00243        G4double Tmax=2.*M*E*(E+2.*mass)/(mass*mass+M*M+2.*M*Etot);
00244        G4double T=Tmax*pow(sin(angle/2.),2.);
00245        G4double q2=T*(T+2.*M);
00246              q2/=htc2;//1/cm2
00247        G4double RN=1.27e-13*pow(targetA,0.27)*cm;
00248        G4double xN= (RN*RN*q2);
00249        G4double den=(1.+xN/12.);
00250        G4double FN=1./(den*den);
00251        G4double form2=(FN*FN);
00252 
00253 
00254   return form2;
00255 
00256 //cout<<"..................... form2 "<< form2<<endl;
00257 }

G4double G4ScreeningMottCrossSection::GetMom2CM (  )  const [inline]

Definition at line 195 of file G4ScreeningMottCrossSection.hh.

00196 {
00197         return mom2;
00198 }

G4double G4ScreeningMottCrossSection::GetMom2Lab (  )  const [inline]

Definition at line 202 of file G4ScreeningMottCrossSection.hh.

Referenced by G4eSingleCoulombScatteringModel::SampleSecondaries().

00203 {
00204         return momLab2;
00205 }

G4ThreeVector G4ScreeningMottCrossSection::GetNewDirection (  ) 

Definition at line 496 of file G4ScreeningMottCrossSection.cc.

References G4UniformRand, and GetScatteringAngle().

Referenced by G4eSingleCoulombScatteringModel::SampleSecondaries().

00496                                                           {
00497 
00498 G4ThreeVector dir(0.0,0.0,1.0);
00499 
00500         
00501         G4double z1=GetScatteringAngle();
00502 
00503         G4double cost = cos(z1);
00504         G4double sint = sin(z1);
00505         G4double phi  = twopi* G4UniformRand();
00506         G4double dirx = sint*cos(phi);
00507         G4double diry = sint*sin(phi);
00508         G4double dirz = cost;
00509 
00510 
00511         //.......set Trc
00512         G4double etot=tkinLab+mass;
00513         G4double mass2=targetMass;
00514         Trec=(1.0 - cost)* mass2*(etot*etot - mass*mass )/
00515                               (mass*mass + mass2*mass2+ 2.*mass2*etot);
00516         
00517 
00518 
00519         dir.set(dirx,diry,dirz);
00520 
00521 
00522 
00523 
00524 return dir;
00525 
00526 }

G4double G4ScreeningMottCrossSection::GetScatteringAngle (  ) 

Definition at line 422 of file G4ScreeningMottCrossSection.cc.

References AngleDistribution(), G4UniformRand, and GetScreeningAngle().

Referenced by GetNewDirection().

00423 {
00424   
00425 
00426 //cout<<" G4ScreeningMottCrossSection::SampleCosineTheta ............."<<endl;  
00427 
00428         if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
00429 
00430         G4double anglemin=std::acos(cosTetMinNuc);      
00431         G4double anglemax= std::acos(cosTetMaxNuc);
00432 
00433         const G4double limit = 1.e-9;
00434         if(anglemin < limit) {
00435           anglemin = GetScreeningAngle()/10.;
00436           if(anglemin < limit) { anglemin = limit; }
00437         }
00438 
00439 //      cout<<"................ tkinLab  "<< G4BestUnit(tkinLab,"Energy") << " anglemin=  "<<anglemin<<endl;
00440         //cout<<"anglemax=  "<<anglemax<<endl;
00441         G4double r =G4UniformRand();
00442 
00443         G4double loganglemin=log10(anglemin);
00444         G4double loganglemax=log10(anglemax);
00445         G4double logdangle=0.01;
00446 
00447         G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);
00448 
00449         std::vector<G4double> angle;
00450         std::vector<G4double> tet;
00451         std::vector<G4double> dangle;
00452 
00453               for(G4int i=0; i<=bins; i++ ){
00454                      tet.push_back(0);
00455                      dangle.push_back(0);
00456                      angle.push_back(pow(10.,loganglemin+logdangle*i));
00457                         }
00458 
00459 
00460                 G4int  dim = tet.size();
00461                 G4double scattangle=0;
00462                 G4double y=0;
00463                 G4double dy=0;
00464                 G4double area=0;
00465 
00466                 for(G4int i=0; i<dim;i++){
00467 
00468                         if(i!=dim-1){
00469                         dangle[i]=(angle[i+1]-angle[i]);
00470                         tet[i]=(angle[i+1]+angle[i])/2.;
00471                         }else if(i==dim-1){
00472                                 break;
00473                         }
00474 
00475                         y+=AngleDistribution(tet[i])*dangle[i];
00476                         dy= y-area ;
00477                         area=y;
00478 
00479                         if(r >=y-dy && r<=y+dy ){       
00480                         scattangle= angle[i] +G4UniformRand()*dangle[i];
00481                         //cout<<"y "<<y <<" r  "<< r  << " y/r "<<y/r << endl;
00482                         break;
00483 
00484                         }       
00485                 
00486                 }
00487 
00488 
00489   return scattangle;
00490 
00491 
00492 }

G4double G4ScreeningMottCrossSection::GetScreeningAngle (  ) 

Definition at line 157 of file G4ScreeningMottCrossSection.cc.

References G4INCL::Math::pi, and SetScreeningCoefficient().

Referenced by GetScatteringAngle(), NuclearCrossSection(), and SetupKinematic().

00157                                                        {
00158 
00159                 SetScreeningCoefficient();
00160 
00161                 //cout<<" As "<<As<<endl;
00162                 if(As < 0.0) { As = 0.0; }
00163                 else if(As > 1.0) { As = 1.0; }
00164 
00165                 G4double screenangle=2.*asin(sqrt(As));
00166 
00167         //      cout<<"  screenangle  "<<  screenangle <<endl;
00168         
00169                 if(screenangle>=pi) screenangle=pi;
00170 
00171         
00172 return screenangle;
00173 
00174 }

G4double G4ScreeningMottCrossSection::GetScreeningCoefficient (  )  const [inline]

Definition at line 217 of file G4ScreeningMottCrossSection.hh.

00218 {
00219         return As;
00220 }

G4double G4ScreeningMottCrossSection::GetTotalCross (  )  const [inline]

Definition at line 225 of file G4ScreeningMottCrossSection.hh.

Referenced by G4eSingleCoulombScatteringModel::SampleSecondaries().

00226 {
00227         return TotalCross;
00228 }

G4double G4ScreeningMottCrossSection::GetTrec (  )  const [inline]

Definition at line 210 of file G4ScreeningMottCrossSection.hh.

Referenced by G4eSingleCoulombScatteringModel::SampleSecondaries().

00211 {
00212         return Trec;
00213 }

void G4ScreeningMottCrossSection::Initialise ( const G4ParticleDefinition *  ,
G4double  cosThetaLim 
)

Definition at line 124 of file G4ScreeningMottCrossSection.cc.

References DBL_MAX, DBL_MIN, and SetupParticle().

Referenced by G4eSingleCoulombScatteringModel::Initialise().

00126 {
00127 
00128         SetupParticle(p);
00129         tkin = targetZ = mom2 = DBL_MIN;
00130         ecut = etag = DBL_MAX;
00131         particle = p;
00132         cosThetaMin = CosThetaLim; 
00133 
00134 }

G4double G4ScreeningMottCrossSection::McFcorrection ( G4double   ) 

Definition at line 261 of file G4ScreeningMottCrossSection.cc.

References G4INCL::Math::pi.

Referenced by AngleDistribution(), and NuclearCrossSection().

00261                                                                   {
00262 
00263 
00264         G4double  beta2=1./invbeta2;
00265         G4double sintmezzi=std::sin(angle/2.);
00266         G4double sin2tmezzi = sintmezzi*sintmezzi;
00267         G4double R=1.-beta2*sin2tmezzi + targetZ*alpha*beta*pi*sintmezzi*(1.-sintmezzi);
00268 
00269 
00270 return R;
00271 }

G4double G4ScreeningMottCrossSection::NuclearCrossSection (  ) 

Definition at line 303 of file G4ScreeningMottCrossSection.cc.

References FormFactor2ExpHof(), GetScreeningAngle(), McFcorrection(), G4INCL::Math::pi, and RatioMottRutherford().

Referenced by G4eSingleCoulombScatteringModel::ComputeCrossSectionPerAtom().

00304 {
00305          if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
00306 
00307         TotalCross=0;
00308 
00309         G4double anglemin =std::acos(cosTetMinNuc);
00310         G4double anglemax =std::acos(cosTetMaxNuc); 
00311 
00312         const G4double limit = 1.e-9;
00313         if(anglemin < limit) {
00314           anglemin = GetScreeningAngle()/10.;
00315           if(anglemin < limit) { anglemin = limit; }
00316         }
00317 
00318         //cout<<" anglemin  "<< anglemin <<endl;
00319 
00320         G4double loganglemin=log10(anglemin);
00321         G4double loganglemax=log10(anglemax);
00322         G4double logdangle=0.01;
00323 
00324         G4int bins=(G4int)((loganglemax-loganglemin)/logdangle);
00325 
00326 
00327         vector<G4double> angle;
00328         vector<G4double> tet;
00329         vector<G4double> dangle;
00330         vector<G4double> cross;
00331 
00332               for(G4int i=0; i<=bins; i++ ){
00333                      tet.push_back(0);
00334                      dangle.push_back(0);
00335                      angle.push_back(pow(10.,loganglemin+logdangle*i));
00336                      cross.push_back(0);
00337                         }
00338 
00339 
00340                 G4int  dim = tet.size();
00341                 //cout<<"dim--- "<<dim<<endl;
00342 
00343 
00344                 for(G4int i=0; i<dim;i++){
00345 
00346                         if(i!=dim-1){
00347                         dangle[i]=(angle[i+1]-angle[i]);
00348                         tet[i]=(angle[i+1]+angle[i])/2.;
00349                         }else if(i==dim-1){
00350                                 break;
00351                         }
00352 
00353 
00354                         G4double R=0;
00355                         G4double F2=FormFactor2ExpHof(tet[i]);
00356                         
00357                         if (coeffb[0][0]!=0){
00358                         //cout<<" Mott....targetZ "<< targetZ<<endl;    
00359                         R=RatioMottRutherford(tet[i]);
00360                         }
00361 
00362                         else if (coeffb[0][0]==0){
00363                        // cout<<" McF.... targetZ "<< targetZ<<endl;
00364                         R=McFcorrection(tet[i]);
00365                         }
00366 
00367 
00368 //cout<<"----------------- R "<<R<<" F2 "<<F2<<endl;
00369 
00370 
00371 //                cout<<"angle "<<tet[i] << " F2 "<<F2<<endl;
00372 
00373                         G4double e4=e2*e2;
00374                         G4double den=2.*As+2.*pow(sin(tet[i]/2.),2.);
00375                         G4double func=1./(den*den);
00376 
00377                         G4double fatt= targetZ/(mu_rel*gamma*beta*beta);
00378                         G4double sigma=e4*fatt*fatt*func;
00379                         cross[i]=F2*R*sigma;
00380                         G4double pi2sintet=2.*pi*sin(tet[i]);
00381 
00382                         TotalCross+=pi2sintet*cross[i]*dangle[i];
00383                         }//end integral
00384 
00385 
00386 //cout<< "ok ......... TotalCross "<<TotalCross<<endl;
00387 
00388 
00389 return TotalCross;
00390 
00391 }

G4double G4ScreeningMottCrossSection::RatioMottRutherford ( G4double   ) 

Definition at line 273 of file G4ScreeningMottCrossSection.cc.

Referenced by AngleDistribution(), and NuclearCrossSection().

00273                                                                        {
00274 
00275 
00276         G4double R=0;
00277         G4double fcost=std::sqrt((1. -cos(angle)));
00278         G4double a[5];
00279         G4double shift=0.7181228;
00280         G4double beta0= beta -shift;
00281 
00282               for(G4int j=0 ;j<=4;j++){
00283                 a[j]=0;
00284                 }
00285 
00286                 for(G4int j=0 ;j<=4;j++){
00287                         for(G4int k=0;k<=5;k++ ){  
00288                                 a[j]+=coeffb[j][k]*pow(beta0,k);
00289                         }
00290                 }
00291 
00292         for(G4int j=0 ;j<=4 ;j++){
00293                 R+=a[j]* pow(fcost,j);
00294                 }
00295 
00296 
00297 
00298 return R;
00299 
00300 }

void G4ScreeningMottCrossSection::SetScreeningCoefficient (  ) 

Definition at line 136 of file G4ScreeningMottCrossSection.cc.

Referenced by GetScreeningAngle().

00136                                                          {
00137 
00138         G4double alpha2=alpha*alpha;
00139         //Bohr radius 
00140         G4double a0=  Bohr_radius  ;//0.529e-8*cm;
00141         //Thomas-Fermi screening length
00142         G4double aU=0.88534*a0/pow(targetZ,1./3.);
00143         G4double twoR2=aU*aU;
00144 
00145         G4double factor= 1.13 + 3.76*targetZ*targetZ*invbeta2*alpha2;
00146         As=0.25*(htc2)/(twoR2*mom2)*factor;
00147 
00148 
00149 
00150 
00151 //cout<<"0k .........................As  "<<As<<endl;
00152 
00153 }

void G4ScreeningMottCrossSection::SetupKinematic ( G4double  kinEnergy,
G4double  Z 
)

Definition at line 177 of file G4ScreeningMottCrossSection.cc.

References G4NistManager::GetAtomicMassAmu(), G4NucleiProperties::GetNuclearMass(), GetScreeningAngle(), and G4MottCoefficients::SetMottCoeff().

Referenced by G4eSingleCoulombScatteringModel::ComputeCrossSectionPerAtom().

00178 {
00179 
00180                 //...Target
00181                 G4int iz = G4int(Z);
00182                 G4double A = fNistManager->GetAtomicMassAmu(iz);
00183                 G4int ia = G4int(A);
00184                 G4double mass2 = G4NucleiProperties::GetNuclearMass(ia, iz);
00185 
00186                 targetZ = Z;
00187                 targetA = fNistManager->GetAtomicMassAmu(iz);
00188                 targetMass= mass2;
00189 
00190                 mottcoeff->SetMottCoeff(targetZ, coeffb);
00191 
00192         //cout<<"......... targetA "<< targetA <<endl;
00193         //cout<<"......... targetMass "<< targetMass/MeV <<endl;
00194 
00195                 // incident particle lab
00196                 tkinLab = ekin;
00197                 momLab2 = tkinLab*(tkinLab + 2.0*mass);
00198                 invbetaLab2 = 1.0 +  mass*mass/momLab2;
00199 
00200                 G4double etot = tkinLab + mass;
00201                 G4double ptot = sqrt(momLab2);
00202                 G4double m12  = mass*mass;
00203                 
00204 
00205         // relativistic reduced mass from publucation
00206         // A.P. Martynenko, R.N. Faustov, Teoret. mat. Fiz. 64 (1985) 179
00207         
00208                 //incident particle & target nucleus
00209                 G4double Ecm=sqrt(m12 + mass2*mass2 + 2.0*etot*mass2);
00210                  mu_rel=mass*mass2/Ecm;
00211                 G4double momCM= ptot*mass2/Ecm;
00212                 // relative system
00213                 mom2 = momCM*momCM;
00214                 invbeta2 = 1.0 +  mu_rel*mu_rel/mom2;
00215                 tkin = momCM*sqrt(invbeta2) - mu_rel;//Ekin of mu_rel
00216                         G4double  beta2=1./invbeta2;
00217                 beta=std::sqrt(beta2) ;
00218                         G4double gamma2= 1./(1.-beta2);
00219                 gamma=std::sqrt(gamma2);
00220 
00221 
00222 
00223                //.........................................................
00224 
00225 
00226                 G4double screenangle=GetScreeningAngle()/10.;
00227                 //cout<<" screenangle [rad] "<<screenangle/rad <<endl;
00228 
00229                 cosTetMinNuc =min( cosThetaMin ,cos(screenangle));
00230                 cosTetMaxNuc =cosThetaMax;
00231         
00232 //cout<<"ok..............mu_rel "<<mu_rel<<endl;
00233 
00234 }

void G4ScreeningMottCrossSection::SetupParticle ( const G4ParticleDefinition *   )  [inline]

Definition at line 184 of file G4ScreeningMottCrossSection.hh.

References G4ParticleDefinition::GetPDGMass(), and G4ParticleDefinition::GetPDGSpin().

Referenced by Initialise(), and G4eSingleCoulombScatteringModel::SetupParticle().

00185 {
00186         particle = p;
00187         mass = particle->GetPDGMass();
00188         spin = particle->GetPDGSpin();
00189                 if(0.0 != spin) { spin = 0.5; }
00190         tkin = 0.0;
00191 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:21 2013 for Geant4 by  doxygen 1.4.7