G4AntiNuclElastic Class Reference

#include <G4AntiNuclElastic.hh>

Inheritance diagram for G4AntiNuclElastic:

G4HadronElastic G4HadronicInteraction

Public Member Functions

 G4AntiNuclElastic ()
virtual ~G4AntiNuclElastic ()
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4double SampleThetaCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4double SampleThetaLab (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4double CalculateParticleBeta (const G4ParticleDefinition *particle, G4double momentum)
G4double CalculateZommerfeld (G4double beta, G4double Z1, G4double Z2)
G4double CalculateAm (G4double momentum, G4double n, G4double Z)
G4double DampFactor (G4double z)
G4double BesselJzero (G4double z)
G4double BesselJone (G4double z)
G4double BesselOneByArg (G4double z)
G4double GetcosTeta1 (G4double plab, G4int A)
G4ComponentAntiNuclNuclearXS * GetComponentCrossSection ()

Detailed Description

Definition at line 50 of file G4AntiNuclElastic.hh.

Constructor & Destructor Documentation

G4AntiNuclElastic::G4AntiNuclElastic (  ) 

Definition at line 55 of file G4AntiNuclElastic.cc.

References G4Alpha::Alpha(), G4AntiAlpha::AntiAlpha(), G4AntiDeuteron::AntiDeuteron(), G4AntiHe3::AntiHe3(), G4AntiNeutron::AntiNeutron(), G4AntiProton::AntiProton(), G4AntiTriton::AntiTriton(), G4Deuteron::Deuteron(), G4Neutron::Neutron(), and G4Proton::Proton().

00056   : G4HadronElastic("AntiAElastic")
00057 {
00058   //V.Ivanchenko commented out 
00059   //SetMinEnergy( 0.1*GeV );
00060   //SetMaxEnergy( 10.*TeV );
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 }

G4AntiNuclElastic::~G4AntiNuclElastic (  )  [virtual]

Definition at line 92 of file G4AntiNuclElastic.cc.

00093 {
00094   delete cs;  
00095 }

Member Function Documentation

G4double G4AntiNuclElastic::BesselJone ( G4double  z  ) 

Definition at line 583 of file G4AntiNuclElastic.cc.

Referenced by BesselOneByArg().

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 }

G4double G4AntiNuclElastic::BesselJzero ( G4double  z  ) 

Definition at line 532 of file G4AntiNuclElastic.cc.

Referenced by SampleInvariantT().

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 }

G4double G4AntiNuclElastic::BesselOneByArg ( G4double  z  ) 

Definition at line 630 of file G4AntiNuclElastic.cc.

References BesselJone().

Referenced by SampleInvariantT().

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 }

G4double G4AntiNuclElastic::CalculateAm ( G4double  momentum,
G4double  n,
G4double  Z 
)

Definition at line 516 of file G4AntiNuclElastic.cc.

References G4Pow::A13(), and G4Pow::GetInstance().

Referenced by SampleInvariantT().

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 }

G4double G4AntiNuclElastic::CalculateParticleBeta ( const G4ParticleDefinition *  particle,
G4double  momentum 
)

Definition at line 493 of file G4AntiNuclElastic.cc.

References G4ParticleDefinition::GetPDGMass().

Referenced by SampleInvariantT().

00495 {
00496   G4double mass = particle->GetPDGMass();
00497   G4double a    = momentum/mass;
00498   fBeta         = a/std::sqrt(1+a*a);
00499 
00500   return fBeta; 
00501 }

G4double G4AntiNuclElastic::CalculateZommerfeld ( G4double  beta,
G4double  Z1,
G4double  Z2 
)

Definition at line 507 of file G4AntiNuclElastic.cc.

Referenced by SampleInvariantT().

00508 {
00509   fZommerfeld = fine_structure_const*Z1*Z2/beta;
00510 
00511   return fZommerfeld; 
00512 }

G4double G4AntiNuclElastic::DampFactor ( G4double  z  ) 

Definition at line 473 of file G4AntiNuclElastic.cc.

Referenced by SampleInvariantT().

00474 {
00475    G4double df;
00476    G4double  f3 = 6.; // first factorials
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 }

G4ComponentAntiNuclNuclearXS * G4AntiNuclElastic::GetComponentCrossSection (  )  [inline]

Definition at line 126 of file G4AntiNuclElastic.hh.

Referenced by G4HadronQElasticPhysics::ConstructProcess(), G4HadronHElasticPhysics::ConstructProcess(), G4HadronElasticPhysics::ConstructProcess(), and G4HadronDElasticPhysics::ConstructProcess().

00127 {
00128   return cs;
00129 }

G4double G4AntiNuclElastic::GetcosTeta1 ( G4double  plab,
G4int  A 
)

Definition at line 649 of file G4AntiNuclElastic.cc.

References G4Pow::GetInstance(), and G4Pow::Z23().

Referenced by SampleInvariantT().

00650 {
00651 
00652 // G4double p0 =G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
00653   G4double p0 = 1.*hbarc/fermi;
00654 //G4double cteta1 = 1.0 - p0*p0/2.0 * pow(A,2./3.)/(plab*plab);
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 }

G4double G4AntiNuclElastic::SampleInvariantT ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
) [virtual]

Reimplemented from G4HadronElastic.

Definition at line 99 of file G4AntiNuclElastic.cc.

References BesselJzero(), BesselOneByArg(), CalculateAm(), CalculateParticleBeta(), CalculateZommerfeld(), DampFactor(), G4UniformRand, G4ComponentAntiNuclNuclearXS::GetAntiHadronNucleonTotCrSc(), G4ParticleDefinition::GetBaryonNumber(), GetcosTeta1(), G4ComponentAntiNuclNuclearXS::GetElasticElementCrossSection(), G4Pow::GetInstance(), G4NucleiProperties::GetNuclearMass(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4ComponentAntiNuclNuclearXS::GetTotalElementCrossSection(), G4He3::He3(), CLHEP::detail::n, G4INCL::Math::pi, sqr(), G4Triton::Triton(), and G4Pow::Z13().

Referenced by SampleThetaCMS(), and SampleThetaLab().

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   //transform to CMS
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)    // Uzhi 24 Nov. 2011
00139   {return fTmax*G4UniformRand();}                              // Uzhi 24 Nov. 2011
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;     //   /hbarc; 
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 // G4cout<<" XselastHadron " << XsElastHad << "  XsCol "<< XsCoulomb <<G4endl; 
00164 // G4cout <<"  XsTotal" << XstotalHad <<G4endl;  
00165 // G4cout<<"XsInel"<< XstotalHad-XsElastHad<<G4endl;        
00166 
00167   if(G4UniformRand() < CoulombProb)
00168   {  // Simulation of Coulomb scattering
00169 
00170    G4double phi = twopi * G4UniformRand();
00171    G4double Ksi =  G4UniformRand();
00172 
00173    G4double par1 = 2.*(1.+Am)/(1.+ctet1);
00174 
00175 // ////sample ThetaCMS in Coulomb part
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 //   G4double Qmax = 2.*ptot*197.33;    // in fm^-1
00194    G4double Qmax = 2.*3.0*197.33;    // in fm^-1
00195    G4double Amag = 70*70;          //  A1 in Magora funct:A1*exp(-q*A2)
00196    G4double SlopeMag = 2.*3.0;    // A2 in Magora funct:A1*exp(-q*A2)
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 //   G4double Mp2=sqr(theDef->GetPDGMass()/GeV ); 
00217 
00218 // change 30 October
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    Ref2=0.8/std::sqrt(std::sqrt(S-4.*Mp2)) + 0.55; 
00235    if(S>1000.) Ref2=0.62+0.02*std::log(S) ;
00236    ceff2 = 0.035/(sqr(S-4.3)+0.4) + 0.085 * std::log(S) ;
00237    if(S>1000.) ceff2 = 0.005 * std::log(S) + 0.29;
00238 */
00239 
00240      Ref2=Ref2*Ref2;
00241      ceff2 = ceff2*ceff2;
00242 
00243      SlopeMag = 0.5;        // Uzhi
00244      Amag= 1.;              // Uzhi
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 // G4cout<<" Ref  "<<fRef<<" c_eff "<<fceff<< " rho "<< rho<<G4endl; 
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;                // fm -> MeV^2
00369    }
00370 
00371    G4double cosTet=1.0-T/(2.*ptot*ptot);
00372    if(cosTet >  1.0 ) cosTet= 1.;          // Uzhi 30 Nov. 
00373    if(cosTet < -1.0 ) cosTet=-1.;          // Uzhi 30 Nov. 
00374    fTetaCMS=std::acos(cosTet);
00375 
00376    return T;
00377 }

G4double G4AntiNuclElastic::SampleThetaCMS ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)

Definition at line 381 of file G4AntiNuclElastic.cc.

References G4cout, G4endl, G4UniformRand, SampleInvariantT(), and G4HadronicInteraction::verboseLevel.

00383 { 
00384   G4double T;
00385   T =  SampleInvariantT( p, plab,  Z,  A);
00386 
00387    // NaN finder
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.)                             // Uzhi 24 Nov. 2011
00402   {
00403    G4double cosTet=1.0-T/(2.*fptot*fptot);
00404    if(cosTet >  1.0 ) cosTet= 1.;          // Uzhi 30 Nov. 
00405    if(cosTet < -1.0 ) cosTet=-1.;          // Uzhi 30 Nov. 
00406    fTetaCMS=std::acos(cosTet); 
00407    return fTetaCMS;
00408   } else                                    // Uzhi 24 Nov. 2011
00409   {                                         // Uzhi 24 Nov. 2011
00410    return 2.*G4UniformRand()-1.;            // Uzhi 24 Nov. 2011
00411   }                                         // Uzhi 24 Nov. 2011
00412 }

G4double G4AntiNuclElastic::SampleThetaLab ( const G4ParticleDefinition *  p,
G4double  plab,
G4int  Z,
G4int  A 
)

Definition at line 417 of file G4AntiNuclElastic.cc.

References G4cout, G4endl, G4UniformRand, G4ParticleDefinition::GetPDGMass(), SampleInvariantT(), and G4HadronicInteraction::verboseLevel.

00419 { 
00420   G4double T; 
00421   T = SampleInvariantT( p, plab,  Z,  A);
00422 
00423  // NaN finder
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;}             // Uzhi 24 Nov. 2011
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 }

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