Geant4: G4ComponentAntiNuclNuclearXS.cc Source 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 //  Calculation of the total, elastic and inelastic cross-sections
00027 //  of anti-nucleon and anti-nucleus interactions with nuclei
00028 //  based on Glauber approach and V. Grishine formulaes for
00029 //  interpolations (ref. V.M.Grichine, Eur.Phys.J., C62(2009) 399;
00030 //  NIM, B267 (2009) 2460) and our parametrization of hadron-nucleon
00031 //  cross-sections
00032 // 
00033 // 
00034 //   Created by A.Galoyan and V. Uzhinsky, 18.11.2010
00035 
00036 
00037 #include "G4ComponentAntiNuclNuclearXS.hh"
00038 
00039 #include "G4PhysicalConstants.hh"
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4ParticleTable.hh"
00042 #include "G4IonTable.hh"
00043 #include "G4ParticleDefinition.hh"
00044 
00046 
00047 G4ComponentAntiNuclNuclearXS::G4ComponentAntiNuclNuclearXS() 
00048 : G4VComponentCrossSection("AntiAGlauber"), fUpperLimit(10000*GeV),
00049   fLowerLimit(10*MeV), fRadiusEff(0.0), fRadiusNN2(0.0),
00050   fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0),
00051   fAntiHadronNucleonTotXsc(0.0), fAntiHadronNucleonElXsc(0.0),
00052   Elab(0.0), S(0.0), SqrtS(0) 
00053 {
00054   theAProton       = G4AntiProton::AntiProton();
00055   theANeutron      = G4AntiNeutron::AntiNeutron();
00056   theADeuteron     = G4AntiDeuteron::AntiDeuteron();
00057   theATriton       = G4AntiTriton::AntiTriton();
00058   theAAlpha        = G4AntiAlpha::AntiAlpha();
00059   theAHe3          = G4AntiHe3::AntiHe3();
00060 
00061  Mn       = 0.93827231;           // GeV
00062  b0       = 11.92;                // GeV^(-2)
00063  b2       = 0.3036;               // GeV^(-2)
00064  SqrtS0   = 20.74;                // GeV
00065  S0       = 33.0625;              // GeV^2
00066 
00067 }
00068 
00070 //
00071 //
00072 
00073 G4ComponentAntiNuclNuclearXS::~G4ComponentAntiNuclNuclearXS()
00074 {
00075 }
00076 
00078 
00079 
00080 
00082 //
00083 // Calculation of total CrossSection of Anti-Nucleus - Nucleus 
00084 
00085 
00086 G4double G4ComponentAntiNuclNuclearXS::GetTotalElementCrossSection
00087 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4double A)
00088 {
00089   G4double xsection,   sigmaTotal, sigmaElastic;
00090 
00091  const G4ParticleDefinition* theParticle = aParticle;
00092 
00093     sigmaTotal        = GetAntiHadronNucleonTotCrSc(theParticle,kinEnergy);
00094     sigmaElastic      = GetAntiHadronNucleonElCrSc(theParticle,kinEnergy);
00095 
00096 // calculation of squared radius of  NN-collision
00097    fRadiusNN2=sigmaTotal*sigmaTotal*0.1/(8.*sigmaElastic*pi) ;  //fm^2   
00098 
00099 // calculation of effective nuclear radius for Pbar and Nbar interactions (can be changed)
00100 
00101   //A.R. 29-Jan-2013 : use antiprotons/antineutrons as the default case,
00102   //                   to be used for instance, as first approximation
00103   //                   without validation, for anti-hyperons. 
00104   //if ( (theParticle == theAProton) || (theParticle == theANeutron) ) {   
00105      if(A==1)
00106      { fTotalXsc = sigmaTotal * millibarn;
00107         return fTotalXsc;  }
00108  
00109    fRadiusEff = 1.34*std::pow(A,0.23)+1.35/std::pow(A,1./3.);   //fm
00110   
00111    if( (Z==1) && (A==2) ) fRadiusEff = 3.800;     //fm
00112    if( (Z==1) && (A==3) ) fRadiusEff = 3.300;  
00113    if( (Z==2) && (A==3) ) fRadiusEff = 3.300;  
00114    if( (Z==2) && (A==4) ) fRadiusEff = 2.376;     
00115  //}
00116       
00117 //calculation of effective nuclear radius for AntiDeuteron interaction (can be changed)
00118   if (theParticle == theADeuteron) 
00119  { fRadiusEff = 1.46 * std::pow(A,0.21) + 1.45 / std::pow(A,1./3.);
00120 
00121     if( (Z==1) && (A==2) ) fRadiusEff = 3.238;     //fm
00122     if( (Z==1) && (A==3) ) fRadiusEff = 3.144;     
00123     if( (Z==2) && (A==3) ) fRadiusEff = 3.144;      
00124     if( (Z==2) && (A==4) ) fRadiusEff = 2.544;     
00125  }
00126 // calculation of effective nuclear radius for AntiHe3 interaction (can be changed)
00127 
00128   if( (theParticle ==theAHe3) || (theParticle ==theATriton) )
00129  { fRadiusEff = 1.40* std::pow(A,0.21)+1.63/std::pow(A,1./3.);
00130 
00131     if( (Z==1) && (A==2) ) fRadiusEff = 3.144;     //fm
00132     if( (Z==1) && (A==3) ) fRadiusEff = 3.075;  
00133     if( (Z==2) && (A==3) ) fRadiusEff = 3.075;  
00134     if( (Z==2) && (A==4) ) fRadiusEff = 2.589;  
00135   }
00136 
00137 //calculation of effective nuclear radius for AntiAlpha interaction (can be changed)
00138 
00139   if (theParticle == theAAlpha) 
00140  {
00141   fRadiusEff = 1.35* std::pow(A,0.21)+1.1/std::pow(A,1./3.);
00142   
00143     if( (Z==1) && (A==2) ) fRadiusEff = 2.544;     //fm
00144     if( (Z==1) && (A==3) ) fRadiusEff = 2.589;   
00145     if( (Z==2) && (A==3) ) fRadiusEff = 2.589;   
00146     if( (Z==2) && (A==4) ) fRadiusEff = 2.241;    
00147   
00148  }
00149 
00150    G4double R2 = fRadiusEff*fRadiusEff;
00151    G4double REf2  = R2+fRadiusNN2;
00152    G4double ApAt = std::abs(theParticle->GetBaryonNumber())  *  A;
00153 
00154  xsection = 2*pi*REf2*10.*std::log(1+(ApAt*sigmaTotal/(2*pi*REf2*10.)));  //mb
00155  xsection =xsection *millibarn; 
00156  fTotalXsc   = xsection;
00157 
00158   return fTotalXsc; 
00159 }
00160 
00161 
00163 // 
00164 // Calculation of total CrossSection of Anti-Nucleus - Nucleus 
00166 G4double G4ComponentAntiNuclNuclearXS::GetTotalIsotopeCrossSection
00167 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4int A )
00168 { return GetTotalElementCrossSection(aParticle, kinEnergy, Z, (G4double) A);  }
00169 
00171 // Calculation of inelastic CrossSection of Anti-Nucleus - Nucleus
00173 
00174 G4double G4ComponentAntiNuclNuclearXS::GetInelasticElementCrossSection
00175 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4double A)
00176 {
00177   G4double  inelxsection,  sigmaTotal, sigmaElastic;
00178 
00179   const G4ParticleDefinition* theParticle = aParticle;
00180 
00181     sigmaTotal        = GetAntiHadronNucleonTotCrSc(theParticle,kinEnergy);
00182     sigmaElastic      = GetAntiHadronNucleonElCrSc(theParticle,kinEnergy);
00183   
00184 // calculation of sqr of radius NN-collision
00185    fRadiusNN2=sigmaTotal*sigmaTotal*0.1/(8.*sigmaElastic*pi);   // fm^2   
00186 
00187 
00188 // calculation of effective nuclear radius for Pbar and Nbar interaction (can be changed)
00189 
00190   //A.R. 29-Jan-2013 : use antiprotons/antineutrons as the default case,
00191   //                   to be used for instance, as first approximation
00192   //                   without validation, for anti-hyperons. 
00193   //if ( (theParticle == theAProton) || (theParticle == theANeutron) ) {   
00194   if (A==1)
00195       { fInelasticXsc = (sigmaTotal - sigmaElastic) * millibarn;
00196         return fInelasticXsc;  
00197       } 
00198  fRadiusEff = 1.31*std::pow(A, 0.22)+0.9/std::pow(A, 1./3.);  //fm
00199     
00200     if( (Z==1) && (A==2) ) fRadiusEff = 3.582;               //fm
00201     if( (Z==1) && (A==3) ) fRadiusEff = 3.105;               
00202     if( (Z==2) && (A==3) ) fRadiusEff = 3.105;
00203     if( (Z==2) && (A==4) ) fRadiusEff = 2.209;
00204  //}
00205 
00206 //calculation of effective nuclear radius for AntiDeuteron interaction (can be changed)
00207 
00208   if (theParticle ==theADeuteron) 
00209 { 
00210  fRadiusEff = 1.38*std::pow(A, 0.21)+1.55/std::pow(A, 1./3.);
00211   
00212     if( (Z==1) && (A==2) ) fRadiusEff = 3.169;            //fm
00213     if( (Z==1) && (A==3) ) fRadiusEff = 3.066;
00214     if( (Z==2) && (A==3) ) fRadiusEff = 3.066;
00215     if( (Z==2) && (A==4) ) fRadiusEff = 2.498;
00216  }
00217 
00218 //calculation of effective nuclear radius for AntiHe3 interaction (can be changed)
00219 
00220   if( (theParticle ==theAHe3) || (theParticle ==theATriton) )
00221  {
00222   fRadiusEff = 1.34 * std::pow(A, 0.21)+1.51/std::pow(A, 1./3.);
00223   
00224     if( (Z==1) && (A==2) ) fRadiusEff = 3.066;           //fm
00225     if( (Z==1) && (A==3) ) fRadiusEff = 2.973;
00226     if( (Z==2) && (A==3) ) fRadiusEff = 2.973;
00227     if( (Z==2) && (A==4) ) fRadiusEff = 2.508;
00228   
00229  }
00230 
00231 //calculation of effective nuclear radius for AntiAlpha interaction (can be changed)
00232 
00233   if (theParticle == theAAlpha) 
00234  {
00235   fRadiusEff = 1.3*std::pow(A, 0.21)+1.05/std::pow(A, 1./3.);
00236     
00237     if( (Z==1) && (A==2) ) fRadiusEff = 2.498;            //fm
00238     if( (Z==1) && (A==3) ) fRadiusEff = 2.508;
00239     if( (Z==2) && (A==3) ) fRadiusEff = 2.508;
00240     if( (Z==2) && (A==4) ) fRadiusEff = 2.158;
00241  }
00242   G4double R2 = fRadiusEff*fRadiusEff;
00243   G4double REf2  = R2+fRadiusNN2;
00244   G4double  ApAt= std::abs(theParticle->GetBaryonNumber())  *  A;
00245 
00246  inelxsection  = pi*REf2 *10* std::log(1+(ApAt*sigmaTotal/(pi*REf2*10.))); //mb
00247  inelxsection  = inelxsection * millibarn;  
00248    fInelasticXsc =  inelxsection; 
00249    return fInelasticXsc;
00250 }
00251 
00253 //
00254 // Calculates Inelastic Anti-nucleus-Nucleus cross-section   
00255 //
00256 G4double G4ComponentAntiNuclNuclearXS::GetInelasticIsotopeCrossSection
00257 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4int A)
00258 {return GetInelasticElementCrossSection(aParticle, kinEnergy, Z, (G4double) A); }
00259  
00260 
00261 
00263 //
00264 // Calculates elastic Anti-nucleus-Nucleus cross-section  as Total - Inelastic 
00265 //
00266 G4double G4ComponentAntiNuclNuclearXS::GetElasticElementCrossSection
00267 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4double A)
00268 {
00269  fElasticXsc = GetTotalElementCrossSection(aParticle, kinEnergy, Z, A)-
00270    GetInelasticElementCrossSection(aParticle, kinEnergy, Z, A);
00271 
00272  if (fElasticXsc < 0.) fElasticXsc = 0.;
00273 
00274  return fElasticXsc;
00275 }
00276  
00278 //
00279 // Calculates elastic Anti-nucleus-Nucleus cross-section   
00280 //
00281 G4double G4ComponentAntiNuclNuclearXS::GetElasticIsotopeCrossSection
00282 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4int A)
00283 { return GetElasticElementCrossSection(aParticle, kinEnergy, Z, (G4double) A); }
00284 
00285 
00287 // Calculation of  Antihadron - hadron Total Cross-section  
00288 
00289 G4double G4ComponentAntiNuclNuclearXS::GetAntiHadronNucleonTotCrSc
00290 (const G4ParticleDefinition* aParticle, G4double kinEnergy)
00291 {
00292   G4double xsection, Pmass, Energy, momentum;
00293   const G4ParticleDefinition* theParticle = aParticle;
00294   Pmass=theParticle->GetPDGMass();
00295   Energy=Pmass+kinEnergy;
00296   momentum=std::sqrt(Energy*Energy-Pmass*Pmass)/std::abs(theParticle->GetBaryonNumber());
00297   G4double Plab = momentum / GeV;
00298 
00299  G4double   B, SigAss;
00300  G4double   C, d1, d2, d3  ;
00301 
00302  Elab     = std::sqrt(Mn*Mn + Plab*Plab);   // GeV
00303  S        = 2.*Mn*Mn + 2. *Mn*Elab;         // GeV^2
00304  SqrtS    = std::sqrt(S);                   // GeV 
00305 
00306  B        = b0+b2*std::log(SqrtS/SqrtS0)*std::log(SqrtS/SqrtS0); //GeV^(-2)
00307  SigAss   = 36.04 +0.304*std::log(S/S0)*std::log(S/S0);          //mb 
00308  R0       = std::sqrt(0.40874044*SigAss - B);                   //GeV^(-2)
00309  
00310  C        = 13.55;
00311  d1       = -4.47;
00312  d2       = 12.38;
00313  d3       = -12.43;
00314  xsection = SigAss*(1 + 1./(std::sqrt(S-4.*Mn*Mn)) / (std::pow(R0, 3.))
00315   *C* (1+d1/SqrtS+d2/(std::pow(SqrtS,2.))+d3/(std::pow(SqrtS,3.)) ));  
00316 
00317 //  xsection *= millibarn;
00318 
00319   fAntiHadronNucleonTotXsc = xsection;
00320   return fAntiHadronNucleonTotXsc;
00321 }
00322 
00323 
00324 //
00325 // /////////////////////////////////////////////////////////////////////////////////
00326 // Calculation of  Antihadron - hadron Elastic Cross-section  
00327 
00328 G4double G4ComponentAntiNuclNuclearXS :: 
00329 GetAntiHadronNucleonElCrSc(const G4ParticleDefinition* aParticle, G4double kinEnergy)
00330 {
00331  G4double xsection;
00332 
00333  G4double   SigAss;
00334  G4double   C, d1, d2, d3  ;
00335 
00336  GetAntiHadronNucleonTotCrSc(aParticle,kinEnergy);
00337 
00338  SigAss   = 4.5 + 0.101*std::log(S/S0)*std::log(S/S0);            //mb
00339   
00340  C        = 59.27;
00341  d1       = -6.95;
00342  d2       = 23.54;
00343  d3       = -25.34;
00344 
00345  xsection = SigAss* (1 + 1. / (std::sqrt(S-4.*Mn*Mn)) / (std::pow(R0, 3.))
00346   *C* ( 1+d1/SqrtS+d2/(std::pow(SqrtS,2.))+d3/(std::pow(SqrtS,3.)) ));  
00347 
00348 //  xsection *= millibarn;
00349 
00350   fAntiHadronNucleonElXsc = xsection;
00351   return fAntiHadronNucleonElXsc;
00352 }
00353 
00354 void G4ComponentAntiNuclNuclearXS::CrossSectionDescription(std::ostream& outFile) const
00355 {
00356   outFile << "The G4ComponentAntiNuclNuclearXS calculates total,\n"
00357           << "inelastic, elastic cross sections  of anti-nucleons and light \n"
00358           << "anti-nucleus interactions with nuclei using Glauber's approach.\n"  
00359           << "It uses parametrizations of antiproton-proton total and elastic \n"
00360           << "cross sections and Wood-Saxon distribution of nuclear density.\n"   
00361           << "The lower limit is 10 MeV, the upper limit is 10 TeV.   \n"
00362           << "See details in Phys.Lett. B705 (2011) 235. \n";
00363 }