00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067 #include "G4ScreeningMottCrossSection.hh"
00068 #include "G4PhysicalConstants.hh"
00069 #include "G4SystemOfUnits.hh"
00070 #include "G4MottCoefficients.hh"
00071 #include "Randomize.hh"
00072 #include "G4Proton.hh"
00073 #include "G4LossTableManager.hh"
00074 #include "G4NucleiProperties.hh"
00075 #include "G4Element.hh"
00076 #include "G4UnitsTable.hh"
00077
00078
00079
00080
00081 using namespace std;
00082
00083 G4ScreeningMottCrossSection::G4ScreeningMottCrossSection():
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 }
00118
00119
00120 G4ScreeningMottCrossSection::~G4ScreeningMottCrossSection()
00121 {}
00122
00123
00124 void G4ScreeningMottCrossSection::Initialise(const G4ParticleDefinition* p,
00125 G4double CosThetaLim)
00126 {
00127
00128 SetupParticle(p);
00129 tkin = targetZ = mom2 = DBL_MIN;
00130 ecut = etag = DBL_MAX;
00131 particle = p;
00132 cosThetaMin = CosThetaLim;
00133
00134 }
00135
00136 void G4ScreeningMottCrossSection::SetScreeningCoefficient(){
00137
00138 G4double alpha2=alpha*alpha;
00139
00140 G4double a0= Bohr_radius ;
00141
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
00152
00153 }
00154
00155
00156
00157 G4double G4ScreeningMottCrossSection::GetScreeningAngle(){
00158
00159 SetScreeningCoefficient();
00160
00161
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
00168
00169 if(screenangle>=pi) screenangle=pi;
00170
00171
00172 return screenangle;
00173
00174 }
00175
00176
00177 void G4ScreeningMottCrossSection::SetupKinematic(G4double ekin, G4double Z )
00178 {
00179
00180
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
00193
00194
00195
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
00206
00207
00208
00209 G4double Ecm=sqrt(m12 + mass2*mass2 + 2.0*etot*mass2);
00210 mu_rel=mass*mass2/Ecm;
00211 G4double momCM= ptot*mass2/Ecm;
00212
00213 mom2 = momCM*momCM;
00214 invbeta2 = 1.0 + mu_rel*mu_rel/mom2;
00215 tkin = momCM*sqrt(invbeta2) - 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
00228
00229 cosTetMinNuc =min( cosThetaMin ,cos(screenangle));
00230 cosTetMaxNuc =cosThetaMax;
00231
00232
00233
00234 }
00235
00236
00237 G4double G4ScreeningMottCrossSection::FormFactor2ExpHof(G4double angle){
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;
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
00257 }
00258
00259
00260
00261 G4double G4ScreeningMottCrossSection::McFcorrection(G4double angle ){
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 }
00272
00273 G4double G4ScreeningMottCrossSection::RatioMottRutherford(G4double angle){
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 }
00301
00302
00303 G4double G4ScreeningMottCrossSection::NuclearCrossSection()
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
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
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
00359 R=RatioMottRutherford(tet[i]);
00360 }
00361
00362 else if (coeffb[0][0]==0){
00363
00364 R=McFcorrection(tet[i]);
00365 }
00366
00367
00368
00369
00370
00371
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 }
00384
00385
00386
00387
00388
00389 return TotalCross;
00390
00391 }
00392
00393 G4double G4ScreeningMottCrossSection::AngleDistribution(G4double anglein){
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
00404 R=RatioMottRutherford(anglein);
00405 }
00406
00407 else if (coeffb[0][0]==0){
00408
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 }
00419
00420
00421
00422 G4double G4ScreeningMottCrossSection::GetScatteringAngle()
00423 {
00424
00425
00426
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
00440
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
00482 break;
00483
00484 }
00485
00486 }
00487
00488
00489 return scattangle;
00490
00491
00492 }
00493
00494
00495
00496 G4ThreeVector G4ScreeningMottCrossSection::GetNewDirection(){
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
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 }
00527
00528