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 #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;
00062 b0 = 11.92;
00063 b2 = 0.3036;
00064 SqrtS0 = 20.74;
00065 S0 = 33.0625;
00066
00067 }
00068
00070
00071
00072
00073 G4ComponentAntiNuclNuclearXS::~G4ComponentAntiNuclNuclearXS()
00074 {
00075 }
00076
00078
00079
00080
00082
00083
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
00097 fRadiusNN2=sigmaTotal*sigmaTotal*0.1/(8.*sigmaElastic*pi) ;
00098
00099
00100
00101
00102
00103
00104
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.);
00110
00111 if( (Z==1) && (A==2) ) fRadiusEff = 3.800;
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
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;
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
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;
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
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;
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.)));
00155 xsection =xsection *millibarn;
00156 fTotalXsc = xsection;
00157
00158 return fTotalXsc;
00159 }
00160
00161
00163
00164
00166 G4double G4ComponentAntiNuclNuclearXS::GetTotalIsotopeCrossSection
00167 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4int A )
00168 { return GetTotalElementCrossSection(aParticle, kinEnergy, Z, (G4double) A); }
00169
00171
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
00185 fRadiusNN2=sigmaTotal*sigmaTotal*0.1/(8.*sigmaElastic*pi);
00186
00187
00188
00189
00190
00191
00192
00193
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.);
00199
00200 if( (Z==1) && (A==2) ) fRadiusEff = 3.582;
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
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;
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
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;
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
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;
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.)));
00247 inelxsection = inelxsection * millibarn;
00248 fInelasticXsc = inelxsection;
00249 return fInelasticXsc;
00250 }
00251
00253
00254
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
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
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
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);
00303 S = 2.*Mn*Mn + 2. *Mn*Elab;
00304 SqrtS = std::sqrt(S);
00305
00306 B = b0+b2*std::log(SqrtS/SqrtS0)*std::log(SqrtS/SqrtS0);
00307 SigAss = 36.04 +0.304*std::log(S/S0)*std::log(S/S0);
00308 R0 = std::sqrt(0.40874044*SigAss - B);
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
00318
00319 fAntiHadronNucleonTotXsc = xsection;
00320 return fAntiHadronNucleonTotXsc;
00321 }
00322
00323
00324
00325
00326
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);
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
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 }