G4AntiNuclElastic.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4AntiNuclElastic.cc  - A.Galoyan 02.05.2011
00027 // GEANT4 tag $Name: not supported by cvs2svn $
00028 //
00029 // Geant4 Header : G4AntiNuclElastic
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   //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 }
00090 
00092 G4AntiNuclElastic::~G4AntiNuclElastic()
00093 {
00094   delete cs;  
00095 }
00096 
00098 // sample momentum transfer in the CMS system 
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   //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 }
00378 
00380 //  Sample of Theta in CMS
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    // 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 }  
00413 
00414 
00416 //  Sample of Theta in Lab System
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  // 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 }
00470 
00472 //   Calculation of Damp factor
00473  G4double G4AntiNuclElastic::DampFactor(G4double x)
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 }
00488 
00489 
00491 //  Calculation of particle velocity Beta
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 //   Calculation of parameter Zommerfeld
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 // Bessel J0 function based on rational approximation from
00530 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
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 // Bessel J1 function based on rational approximation from
00581 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
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 // return J1(x)/x with special case for small x 
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 // return  angle from which Coulomb scattering is calculated 
00649 G4double G4AntiNuclElastic::GetcosTeta1(G4double plab, G4int A)
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 }
00660 
00661 
00662 
00663 
00664 
00665 
00666 

Generated on Mon May 27 17:47:39 2013 for Geant4 by  doxygen 1.4.7