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 #include "G4AntiNuclElastic.hh"
00034
00035 #include "G4PhysicalConstants.hh"
00036 #include "G4SystemOfUnits.hh"
00037 #include "G4ParticleTable.hh"
00038 #include "G4ParticleDefinition.hh"
00039 #include "G4IonTable.hh"
00040 #include "Randomize.hh"
00041 #include "G4AntiProton.hh"
00042 #include "G4AntiNeutron.hh"
00043 #include "G4AntiDeuteron.hh"
00044 #include "G4AntiAlpha.hh"
00045 #include "G4AntiTriton.hh"
00046 #include "G4AntiHe3.hh"
00047 #include "G4Proton.hh"
00048 #include "G4Neutron.hh"
00049 #include "G4Deuteron.hh"
00050 #include "G4Alpha.hh"
00051 #include "G4Pow.hh"
00052
00053 #include "G4NucleiProperties.hh"
00054
00055 G4AntiNuclElastic::G4AntiNuclElastic()
00056 : G4HadronElastic("AntiAElastic")
00057 {
00058
00059
00060
00061
00062
00063 theAProton = G4AntiProton::AntiProton();
00064 theANeutron = G4AntiNeutron::AntiNeutron();
00065 theADeuteron = G4AntiDeuteron::AntiDeuteron();
00066 theATriton = G4AntiTriton::AntiTriton();
00067 theAAlpha = G4AntiAlpha::AntiAlpha();
00068 theAHe3 = G4AntiHe3::AntiHe3();
00069
00070 theProton = G4Proton::Proton();
00071 theNeutron = G4Neutron::Neutron();
00072 theDeuteron = G4Deuteron::Deuteron();
00073 theAlpha = G4Alpha::Alpha();
00074
00075
00076 cs = new G4ComponentAntiNuclNuclearXS();
00077 fParticle = 0;
00078 fWaveVector = 0.;
00079 fBeta = 0.;
00080 fZommerfeld = 0.;
00081 fAm = 0.;
00082 fTetaCMS = 0.;
00083 fRa = 0.;
00084 fRef = 0.;
00085 fceff = 0.;
00086 fptot = 0.;
00087 fTmax = 0.;
00088 fThetaLab = 0.;
00089 }
00090
00092 G4AntiNuclElastic::~G4AntiNuclElastic()
00093 {
00094 delete cs;
00095 }
00096
00098
00099 G4double G4AntiNuclElastic::SampleInvariantT(const G4ParticleDefinition* particle,
00100 G4double Plab, G4int Z, G4int A)
00101 {
00102 G4double T;
00103 G4double Mproj = particle->GetPDGMass();
00104 G4LorentzVector Pproj(0.,0.,Plab,std::sqrt(Plab*Plab+Mproj*Mproj));
00105 G4double ctet1 = GetcosTeta1(Plab, A);
00106
00107 G4double energy=Pproj.e()-Mproj;
00108
00109 const G4ParticleDefinition* theParticle = particle;
00110
00111 G4ParticleDefinition * theDef = 0;
00112
00113 if(Z == 1 && A == 1) theDef = theProton;
00114 else if (Z == 1 && A == 2) theDef = theDeuteron;
00115 else if (Z == 1 && A == 3) theDef = G4Triton::Triton();
00116 else if (Z == 2 && A == 3) theDef = G4He3::He3();
00117 else if (Z == 2 && A == 4) theDef = theAlpha;
00118
00119
00120 G4double TargMass =G4NucleiProperties::GetNuclearMass(A,Z);
00121
00122
00123
00124 G4LorentzVector lv(0.0,0.0,0.0,TargMass);
00125 lv += Pproj;
00126 G4double S = lv.mag2()/GeV/GeV;
00127
00128 G4ThreeVector bst = lv.boostVector();
00129 Pproj.boost(-bst);
00130
00131 G4ThreeVector p1 = Pproj.vect();
00132 G4double ptot = p1.mag();
00133
00134 fbst = bst;
00135 fptot= ptot;
00136 fTmax = 4.0*ptot*ptot;
00137
00138 if(Plab/std::abs(particle->GetBaryonNumber()) < 100.*MeV)
00139 {return fTmax*G4UniformRand();}
00140
00141 G4double Z1 = particle->GetPDGCharge();
00142 G4double Z2 = Z;
00143
00144 G4double beta = CalculateParticleBeta(particle, ptot);
00145 G4double n = CalculateZommerfeld( beta, Z1, Z2 );
00146 G4double Am = CalculateAm( ptot, n, Z2 );
00147 fWaveVector = ptot;
00148
00149 G4LorentzVector Fproj(0.,0.,0.,0.);
00150 G4double XsCoulomb = sqr(n/fWaveVector)*pi*(1+ctet1)/(1.+Am)/(1.+2.*Am-ctet1);
00151 XsCoulomb=XsCoulomb*0.38938e+6;
00152
00153
00154 G4double XsElastHad =cs->GetElasticElementCrossSection(particle, energy, Z, (G4double)A);
00155 G4double XstotalHad =cs->GetTotalElementCrossSection(particle, energy, Z, (G4double)A);
00156
00157
00158 XsElastHad/=millibarn; XstotalHad/=millibarn;
00159
00160
00161 G4double CoulombProb = XsCoulomb/(XsCoulomb+XsElastHad);
00162
00163
00164
00165
00166
00167 if(G4UniformRand() < CoulombProb)
00168 {
00169
00170 G4double phi = twopi * G4UniformRand();
00171 G4double Ksi = G4UniformRand();
00172
00173 G4double par1 = 2.*(1.+Am)/(1.+ctet1);
00174
00175
00176
00177 G4double cosThetaCMS = (par1*ctet1- Ksi*(1.+2.*Am))/(par1-Ksi);
00178
00179 G4double PtZ=ptot*cosThetaCMS;
00180 Fproj.setPz(PtZ);
00181 G4double PtProjCMS = ptot*std::sqrt(1.0 - cosThetaCMS*cosThetaCMS);
00182 G4double PtX= PtProjCMS * std::cos(phi);
00183 G4double PtY= PtProjCMS * std::sin(phi);
00184 Fproj.setPx(PtX);
00185 Fproj.setPy(PtY);
00186 Fproj.setE(std::sqrt(PtX*PtX+PtY*PtY+PtZ*PtZ+Mproj*Mproj));
00187 T = -(Pproj-Fproj).mag2();
00188 } else
00189
00190 {
00192
00193
00194 G4double Qmax = 2.*3.0*197.33;
00195 G4double Amag = 70*70;
00196 G4double SlopeMag = 2.*3.0;
00197
00198 G4double sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy);
00199
00200 fRa = 1.113*G4Pow::GetInstance()->Z13(A) -
00201 0.227/G4Pow::GetInstance()->Z13(A);
00202 if(A == 3) fRa=1.81;
00203 if(A == 4) fRa=1.37;
00204
00205 if((A>=12.) && (A<27) ) fRa=fRa*0.85;
00206 if((A>=27.) && (A<48) ) fRa=fRa*0.90;
00207 if((A>=48.) && (A<65) ) fRa=fRa*0.95;
00208
00209 G4double Ref2 = 0;
00210 G4double ceff2 =0;
00211 G4double rho = 0;
00212 if ((theParticle == theAProton) || (theParticle == theANeutron))
00213 {
00214 if(theDef == theProton)
00215 {
00216
00217
00218
00219
00220 if(Plab < 610.)
00221 { rho = 1.3347-10.342*Plab/1000.+22.277*Plab/1000.*Plab/1000.-
00222 13.634*Plab/1000.*Plab/1000.*Plab/1000. ;}
00223 if((Plab < 5500.)&&(Plab >= 610.) )
00224 { rho = 0.22; }
00225 if((Plab >= 5500.)&&(Plab < 12300.) )
00226 { rho = -0.32; }
00227 if( Plab >= 12300.)
00228 { rho = 0.135-2.26/(std::sqrt(S)) ;}
00229
00230 Ref2 = 0.35 + 0.9/std::sqrt(std::sqrt(S-4.*0.88))+0.04*std::log(S) ;
00231 ceff2 = 0.375 - 2./S + 0.44/(sqr(S-4.)+1.5) ;
00232
00233
00234
00235
00236
00237
00238
00239
00240 Ref2=Ref2*Ref2;
00241 ceff2 = ceff2*ceff2;
00242
00243 SlopeMag = 0.5;
00244 Amag= 1.;
00245 }
00246
00247 if(Z>2)
00248 { Ref2 = fRa*fRa +2.48*0.01*sig_pbarp*fRa - 2.23e-6*sig_pbarp*sig_pbarp*fRa*fRa;
00249 ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*std::exp(-0.03*sig_pbarp);
00250 }
00251 if( (Z==2)&&(A==4) )
00252 { Ref2 = fRa*fRa -0.46 +0.03*sig_pbarp - 2.98e-6*sig_pbarp*sig_pbarp;
00253 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*std::exp(-0.03*sig_pbarp);
00254 }
00255 if( (Z==1)&&(A==3) )
00256 { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
00257 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp);
00258 }
00259 if( (Z==2)&&(A==3) )
00260 { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
00261 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp);
00262 }
00263 if( (Z==1)&&(A==2) )
00264 {
00265 Ref2 = fRa*fRa - 0.28 + 0.019 * sig_pbarp + 2.06e-6 * sig_pbarp*sig_pbarp;
00266 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*std::exp(-0.03*sig_pbarp);
00267 }
00268 }
00269
00270 if (theParticle == theADeuteron)
00271 {
00272 sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy/2.);
00273 Ref2 = XstotalHad/10./2./pi ;
00274 if(Z>2)
00275 {
00276 ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5 * std::exp(-0.03*sig_pbarp);
00277 }
00278 if(theDef == theProton)
00279 {
00280 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*std::exp(-0.03*sig_pbarp);
00281 }
00282 if(theDef == theDeuteron)
00283 {
00284 ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55 * std::exp(-0.03*sig_pbarp);
00285 }
00286 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
00287 {
00288 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * std::exp(-0.02*sig_pbarp);
00289 }
00290 if(theDef == theAlpha)
00291 {
00292 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * std::exp(-0.02*sig_pbarp);
00293 }
00294 }
00295
00296 if( (theParticle ==theAHe3) || (theParticle ==theATriton) )
00297 {
00298 sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.);
00299 Ref2 = XstotalHad/10./2./pi ;
00300 if(Z>2)
00301 {
00302 ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33*std::exp(-0.03*sig_pbarp);
00303 }
00304 if(theDef == theProton)
00305 {
00306 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp);
00307 }
00308 if(theDef == theDeuteron)
00309 {
00310 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * std::exp(-0.02*sig_pbarp);
00311 }
00312 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
00313 {
00314 ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7 * std::exp(-0.02*sig_pbarp);
00315 }
00316 if(theDef == theAlpha)
00317 {
00318 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * std::exp(-0.03*sig_pbarp);
00319 }
00320 }
00321
00322
00323 if (theParticle == theAAlpha)
00324 {
00325 sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.);
00326 Ref2 = XstotalHad/10./2./pi ;
00327 if(Z>2)
00328 {
00329 ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2 * std::exp(-0.03*sig_pbarp);
00330 }
00331 if(theDef == theProton)
00332 {
00333 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*std::exp(-0.03*sig_pbarp);
00334 }
00335 if(theDef == theDeuteron)
00336 {
00337 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * std::exp(-0.02*sig_pbarp);
00338 }
00339 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
00340 {
00341 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * std::exp(-0.03*sig_pbarp);
00342 }
00343 if(theDef == theAlpha)
00344 {
00345 ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45 * std::exp(-0.03*sig_pbarp);
00346 }
00347 }
00348
00349 fRef=std::sqrt(Ref2);
00350 fceff = std::sqrt(ceff2);
00351
00352
00353
00354 G4double Q = 0.0 ;
00355 G4double BracFunct;
00356 do
00357 {
00358 Q = -std::log(1.-(1.- std::exp(-SlopeMag * Qmax))* G4UniformRand() )/SlopeMag;
00359 G4double x = fRef * Q;
00360 BracFunct = ( ( sqr(BesselOneByArg(x))+sqr(rho/2. * BesselJzero(x)) )
00361 * sqr(DampFactor(pi*fceff*Q))) /(Amag*std::exp(-SlopeMag*Q));
00362
00363 BracFunct = BracFunct * Q * sqr(sqr(fRef));
00364 }
00365 while (G4UniformRand()>BracFunct);
00366
00367 T= sqr(Q);
00368 T*=3.893913e+4;
00369 }
00370
00371 G4double cosTet=1.0-T/(2.*ptot*ptot);
00372 if(cosTet > 1.0 ) cosTet= 1.;
00373 if(cosTet < -1.0 ) cosTet=-1.;
00374 fTetaCMS=std::acos(cosTet);
00375
00376 return T;
00377 }
00378
00380
00381 G4double G4AntiNuclElastic::SampleThetaCMS(const G4ParticleDefinition* p, G4double plab,
00382 G4int Z, G4int A)
00383 {
00384 G4double T;
00385 T = SampleInvariantT( p, plab, Z, A);
00386
00387
00388 if(!(T < 0.0 || T >= 0.0))
00389 {
00390 if (verboseLevel > 0)
00391 {
00392 G4cout << "G4DiffuseElastic:WARNING: A = " << A
00393 << " mom(GeV)= " << plab/GeV
00394 << " S-wave will be sampled"
00395 << G4endl;
00396 }
00397 T = G4UniformRand()*fTmax;
00398
00399 }
00400
00401 if(fptot > 0.)
00402 {
00403 G4double cosTet=1.0-T/(2.*fptot*fptot);
00404 if(cosTet > 1.0 ) cosTet= 1.;
00405 if(cosTet < -1.0 ) cosTet=-1.;
00406 fTetaCMS=std::acos(cosTet);
00407 return fTetaCMS;
00408 } else
00409 {
00410 return 2.*G4UniformRand()-1.;
00411 }
00412 }
00413
00414
00416
00417 G4double G4AntiNuclElastic::SampleThetaLab(const G4ParticleDefinition* p, G4double plab,
00418 G4int Z, G4int A)
00419 {
00420 G4double T;
00421 T = SampleInvariantT( p, plab, Z, A);
00422
00423
00424 if(!(T < 0.0 || T >= 0.0))
00425 {
00426 if (verboseLevel > 0)
00427 {
00428 G4cout << "G4DiffuseElastic:WARNING: A = " << A
00429 << " mom(GeV)= " << plab/GeV
00430 << " S-wave will be sampled"
00431 << G4endl;
00432 }
00433 T = G4UniformRand()*fTmax;
00434 }
00435
00436 G4double phi = G4UniformRand()*twopi;
00437
00438 G4double cost(1.);
00439 if(fTmax > 0.) {cost = 1. - 2.0*T/fTmax;}
00440
00441 G4double sint;
00442 if( cost >= 1.0 )
00443 {
00444 cost = 1.0;
00445 sint = 0.0;
00446 }
00447 else if( cost <= -1.0)
00448 {
00449 cost = -1.0;
00450 sint = 0.0;
00451 }
00452 else
00453 {
00454 sint = std::sqrt((1.0-cost)*(1.0+cost));
00455 }
00456
00457 G4double m1 = p->GetPDGMass();
00458 G4ThreeVector v(sint*std::cos(phi),sint*std::sin(phi),cost);
00459 v *= fptot;
00460 G4LorentzVector nlv(v.x(),v.y(),v.z(),std::sqrt(fptot*fptot + m1*m1));
00461
00462 nlv.boost(fbst);
00463
00464 G4ThreeVector np = nlv.vect();
00465 G4double theta = np.theta();
00466 fThetaLab = theta;
00467
00468 return theta;
00469 }
00470
00472
00473 G4double G4AntiNuclElastic::DampFactor(G4double x)
00474 {
00475 G4double df;
00476 G4double f3 = 6.;
00477
00478 if( std::fabs(x) < 0.01 )
00479 {
00480 df=1./(1.+x*x/f3);
00481 }
00482 else
00483 {
00484 df = x/std::sinh(x);
00485 }
00486 return df;
00487 }
00488
00489
00491
00492
00493 G4double G4AntiNuclElastic::CalculateParticleBeta( const G4ParticleDefinition* particle,
00494 G4double momentum )
00495 {
00496 G4double mass = particle->GetPDGMass();
00497 G4double a = momentum/mass;
00498 fBeta = a/std::sqrt(1+a*a);
00499
00500 return fBeta;
00501 }
00502
00503
00505
00506
00507 G4double G4AntiNuclElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
00508 {
00509 fZommerfeld = fine_structure_const*Z1*Z2/beta;
00510
00511 return fZommerfeld;
00512 }
00513
00515
00516 G4double G4AntiNuclElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
00517 {
00518 G4double k = momentum/hbarc;
00519 G4double ch = 1.13 + 3.76*n*n;
00520 G4double zn = 1.77*k/G4Pow::GetInstance()->A13(Z)*Bohr_radius;
00521 G4double zn2 = zn*zn;
00522 fAm = ch/zn2;
00523
00524 return fAm;
00525 }
00526
00528
00529
00530
00531
00532 G4double G4AntiNuclElastic::BesselJzero(G4double value)
00533 {
00534 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00535
00536 modvalue = std::fabs(value);
00537
00538 if ( value < 8.0 && value > -8.0 )
00539 {
00540 value2 = value*value;
00541
00542 fact1 = 57568490574.0 + value2*(-13362590354.0
00543 + value2*( 651619640.7
00544 + value2*(-11214424.18
00545 + value2*( 77392.33017
00546 + value2*(-184.9052456 ) ) ) ) );
00547
00548 fact2 = 57568490411.0 + value2*( 1029532985.0
00549 + value2*( 9494680.718
00550 + value2*(59272.64853
00551 + value2*(267.8532712
00552 + value2*1.0 ) ) ) );
00553
00554 bessel = fact1/fact2;
00555 }
00556 else
00557 {
00558 arg = 8.0/modvalue;
00559
00560 value2 = arg*arg;
00561
00562 shift = modvalue-0.785398164;
00563
00564 fact1 = 1.0 + value2*(-0.1098628627e-2
00565 + value2*(0.2734510407e-4
00566 + value2*(-0.2073370639e-5
00567 + value2*0.2093887211e-6 ) ) );
00568 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
00569 + value2*(-0.6911147651e-5
00570 + value2*(0.7621095161e-6
00571 - value2*0.934945152e-7 ) ) );
00572
00573 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
00574 }
00575 return bessel;
00576 }
00577
00578
00580
00581
00582
00583 G4double G4AntiNuclElastic::BesselJone(G4double value)
00584 {
00585 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00586
00587 modvalue = std::fabs(value);
00588
00589 if ( modvalue < 8.0 )
00590 {
00591 value2 = value*value;
00592 fact1 = value*(72362614232.0 + value2*(-7895059235.0
00593 + value2*( 242396853.1
00594 + value2*(-2972611.439
00595 + value2*( 15704.48260
00596 + value2*(-30.16036606 ) ) ) ) ) );
00597
00598 fact2 = 144725228442.0 + value2*(2300535178.0
00599 + value2*(18583304.74
00600 + value2*(99447.43394
00601 + value2*(376.9991397
00602 + value2*1.0 ) ) ) );
00603 bessel = fact1/fact2;
00604 }
00605 else
00606 {
00607 arg = 8.0/modvalue;
00608 value2 = arg*arg;
00609
00610 shift = modvalue - 2.356194491;
00611
00612 fact1 = 1.0 + value2*( 0.183105e-2
00613 + value2*(-0.3516396496e-4
00614 + value2*(0.2457520174e-5
00615 + value2*(-0.240337019e-6 ) ) ) );
00616
00617 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
00618 + value2*( 0.8449199096e-5
00619 + value2*(-0.88228987e-6
00620 + value2*0.105787412e-6 ) ) );
00621
00622 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
00623 if (value < 0.0) bessel = -bessel;
00624 }
00625 return bessel;
00626 }
00627
00629
00630 G4double G4AntiNuclElastic::BesselOneByArg(G4double x)
00631 {
00632 G4double x2, result;
00633
00634 if( std::fabs(x) < 0.01 )
00635 {
00636 x *= 0.5;
00637 x2 = x*x;
00638 result = (2.- x2 + x2*x2/6.)/4.;
00639 }
00640 else
00641 {
00642 result = BesselJone(x)/x;
00643 }
00644 return result;
00645 }
00646
00648
00649 G4double G4AntiNuclElastic::GetcosTeta1(G4double plab, G4int A)
00650 {
00651
00652
00653 G4double p0 = 1.*hbarc/fermi;
00654
00655 G4double cteta1 = 1.0 - p0*p0/2.0 * G4Pow::GetInstance()->Z23(A)/(plab*plab);
00657 if(cteta1 < -1.) cteta1 = -1.0;
00658 return cteta1;
00659 }
00660
00661
00662
00663
00664
00665
00666