Skip to content
Snippets Groups Projects
Commit 0d3dfa0b authored by Ivan Vorobyev's avatar Ivan Vorobyev
Browse files

fix merging issues 2

parent 70bd9b76
No related branches found
No related tags found
No related merge requests found
......@@ -52,84 +52,52 @@
#include "G4VComponentCrossSection.hh"
class G4ParticleDefinition;
class G4ComponentAntiNuclNuclearXS : public G4VComponentCrossSection
{
public:
G4ComponentAntiNuclNuclearXS ();
virtual ~G4ComponentAntiNuclNuclearXS ();
virtual
G4double GetTotalIsotopeCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy,
G4int Z, G4int A);
virtual
G4double GetTotalElementCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy,
G4int Z, G4double A);
virtual
G4double GetInelasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy,
G4int Z, G4int A);
virtual
G4double GetInelasticElementCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy,
G4int Z, G4double A);
virtual
G4double GetElasticElementCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy,
G4int Z, G4double A);
virtual
G4double GetElasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy,
G4int Z, G4int A);
virtual
void BuildPhysicsTable(const G4ParticleDefinition&)
{}
virtual
void DumpPhysicsTable(const G4ParticleDefinition&)
{}
virtual void CrossSectionDescription(std::ostream&) const;
// Method for calculation of Anti-Hadron Nucleon Total Cross-section
G4double GetAntiHadronNucleonTotCrSc(const G4ParticleDefinition* aParticle, G4double kinEnergy);
// Method for calculation of Anti-Hadron Nucleon Elastic Cross-section
G4double GetAntiHadronNucleonElCrSc(const G4ParticleDefinition* aParticle, G4double kinEnergy);
// Method for calculation of scaling factor for cross-section
G4double GetScalingFactorCrSc(const G4ParticleDefinition* aParticle, G4double kinEnergy);
private:
class G4ParticleDefinition;
// const G4double fUpperLimit;
// const G4double fLowerLimit;
G4double fRadiusEff; // Effective Radius for AntiNucleus
G4double fRadiusNN2; // Sqr of radius of NN collision
G4double fTotalXsc, fElasticXsc, fInelasticXsc;
G4double fAntiHadronNucleonTotXsc, fAntiHadronNucleonElXsc;
G4double Elab, S, SqrtS ;
G4double Mn, b0, b2, SqrtS0, S0, R0; //parameters for AntiHadron-Nucleon Xsc
G4double fScalingFactorXsc; // scaling factor for cross-sections
G4ParticleDefinition* theAProton;
G4ParticleDefinition* theANeutron;
G4ParticleDefinition* theADeuteron;
G4ParticleDefinition* theATriton;
G4ParticleDefinition* theAAlpha;
G4ParticleDefinition* theAHe3;
class G4ComponentAntiNuclNuclearXS : public G4VComponentCrossSection {
public:
G4ComponentAntiNuclNuclearXS ();
virtual ~G4ComponentAntiNuclNuclearXS ();
virtual G4double GetTotalIsotopeCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy, G4int Z, G4int A);
virtual G4double GetTotalElementCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy, G4int Z, G4double A);
virtual G4double GetInelasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy, G4int Z, G4int A);
virtual G4double GetInelasticElementCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy, G4int Z, G4double A);
virtual G4double GetElasticElementCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy, G4int Z, G4double A);
virtual G4double GetElasticIsotopeCrossSection(const G4ParticleDefinition* aParticle,
G4double kinEnergy, G4int Z, G4int A);
virtual void BuildPhysicsTable(const G4ParticleDefinition&) {}
virtual void DumpPhysicsTable(const G4ParticleDefinition&) {}
virtual void CrossSectionDescription(std::ostream&) const;
// Method for calculation of Anti-Hadron Nucleon Total Cross-section
G4double GetAntiHadronNucleonTotCrSc(const G4ParticleDefinition* aParticle, G4double kinEnergy);
// Method for calculation of Anti-Hadron Nucleon Elastic Cross-section
G4double GetAntiHadronNucleonElCrSc(const G4ParticleDefinition* aParticle, G4double kinEnergy);
// Method for calculation of scaling factor for cross-section
G4double GetScalingFactorCrSc(const G4ParticleDefinition* aParticle, G4double kinEnergy);
private:
G4double fRadiusEff; // Effective Radius for AntiNucleus
G4double fRadiusNN2; // Sqr of radius of NN collision
G4double fTotalXsc, fElasticXsc, fInelasticXsc;
G4double fAntiHadronNucleonTotXsc, fAntiHadronNucleonElXsc;
G4double Elab, S, SqrtS ;
G4double Mn, b0, b2, SqrtS0, S0, R0; // Parameters for AntiHadron-Nucleon Xsc
G4double fScalingFactorXsc; // scaling factor for cross-sections
G4ParticleDefinition* theAProton;
G4ParticleDefinition* theANeutron;
G4ParticleDefinition* theADeuteron;
G4ParticleDefinition* theATriton;
G4ParticleDefinition* theAAlpha;
G4ParticleDefinition* theAHe3;
};
......
......@@ -45,266 +45,228 @@
#include <dlfcn.h>
///////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
G4ComponentAntiNuclNuclearXS::G4ComponentAntiNuclNuclearXS()
: G4VComponentCrossSection("AntiAGlauber"),
// fUpperLimit(10000*GeV), fLowerLimit(10*MeV),
fRadiusEff(0.0), fRadiusNN2(0.0),
fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0),
fAntiHadronNucleonTotXsc(0.0), fAntiHadronNucleonElXsc(0.0),
Elab(0.0), S(0.0), SqrtS(0), fScalingFactorXsc(1.0)
{
theAProton = G4AntiProton::AntiProton();
theANeutron = G4AntiNeutron::AntiNeutron();
theADeuteron = G4AntiDeuteron::AntiDeuteron();
theATriton = G4AntiTriton::AntiTriton();
theAAlpha = G4AntiAlpha::AntiAlpha();
theAHe3 = G4AntiHe3::AntiHe3();
Mn = 0.93827231; // GeV
b0 = 11.92; // GeV^(-2)
b2 = 0.3036; // GeV^(-2)
SqrtS0 = 20.74; // GeV
S0 = 33.0625; // GeV^2
R0 = 1.0; // default value (V.Ivanchenko)
theAProton = G4AntiProton::AntiProton();
theANeutron = G4AntiNeutron::AntiNeutron();
theADeuteron = G4AntiDeuteron::AntiDeuteron();
theATriton = G4AntiTriton::AntiTriton();
theAAlpha = G4AntiAlpha::AntiAlpha();
theAHe3 = G4AntiHe3::AntiHe3();
Mn = 0.93827231; // GeV
b0 = 11.92; // GeV^(-2)
b2 = 0.3036; // GeV^(-2)
SqrtS0 = 20.74; // GeV
S0 = 33.0625; // GeV^2
R0 = 1.0; // default value (V.Ivanchenko)
}
///////////////////////////////////////////////////////////////////////////////////////
//
//
/////////////////////////////////////////////////////////////////////////////
G4ComponentAntiNuclNuclearXS::~G4ComponentAntiNuclNuclearXS()
{
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// Calculation of total CrossSection of Anti-Nucleus - Nucleus
G4double G4ComponentAntiNuclNuclearXS::GetTotalElementCrossSection
(const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4double A)
{
G4double xsection, sigmaTotal, sigmaElastic;
const G4ParticleDefinition* theParticle = aParticle;
sigmaTotal = GetAntiHadronNucleonTotCrSc(theParticle,kinEnergy);
sigmaElastic = GetAntiHadronNucleonElCrSc(theParticle,kinEnergy);
// calculation of squared radius of NN-collision
fRadiusNN2=sigmaTotal*sigmaTotal*0.1/(8.*sigmaElastic*pi) ; //fm^2
G4double xsection, sigmaTotal, sigmaElastic;
const G4ParticleDefinition* theParticle = aParticle;
sigmaTotal = GetAntiHadronNucleonTotCrSc(theParticle,kinEnergy);
sigmaElastic = GetAntiHadronNucleonElCrSc(theParticle,kinEnergy);
// calculation of effective nuclear radius for Pbar and Nbar interactions (can be changed)
// calculation of squared radius of NN-collision
fRadiusNN2=sigmaTotal*sigmaTotal*0.1/(8.*sigmaElastic*pi) ; //fm^2
// calculation of effective nuclear radius for Pbar and Nbar interactions (can be changed)
//A.R. 29-Jan-2013 : use antiprotons/antineutrons as the default case,
// to be used for instance, as first approximation
// without validation, for anti-hyperons.
//if ( (theParticle == theAProton) || (theParticle == theANeutron) ) {
if(A==1)
{ fTotalXsc = sigmaTotal * millibarn;
return fTotalXsc; }
fRadiusEff = 1.34*G4Pow::GetInstance()->powA(A,0.23)+1.35/G4Pow::GetInstance()->powA(A,1./3.); //fm
if( (Z==1) && (A==2) ) fRadiusEff = 3.800; //fm
if( (Z==1) && (A==3) ) fRadiusEff = 3.300;
if( (Z==2) && (A==3) ) fRadiusEff = 3.300;
if( (Z==2) && (A==4) ) fRadiusEff = 2.376;
//}
//calculation of effective nuclear radius for AntiDeuteron interaction (can be changed)
if (theParticle == theADeuteron)
{ fRadiusEff = 1.46 * G4Pow::GetInstance()->powA(A,0.21) + 1.45 / G4Pow::GetInstance()->powA(A,1./3.);
if( (Z==1) && (A==2) ) fRadiusEff = 3.238; //fm
if( (Z==1) && (A==3) ) fRadiusEff = 3.144;
if( (Z==2) && (A==3) ) fRadiusEff = 3.144;
if( (Z==2) && (A==4) ) fRadiusEff = 2.544;
}
// calculation of effective nuclear radius for AntiHe3 interaction (can be changed)
if( (theParticle ==theAHe3) || (theParticle ==theATriton) )
{ fRadiusEff = 1.40* G4Pow::GetInstance()->powA(A,0.21)+1.63/G4Pow::GetInstance()->powA(A,1./3.);
if( (Z==1) && (A==2) ) fRadiusEff = 3.144; //fm
if( (Z==1) && (A==3) ) fRadiusEff = 3.075;
if( (Z==2) && (A==3) ) fRadiusEff = 3.075;
if( (Z==2) && (A==4) ) fRadiusEff = 2.589;
if (A==1) {
fTotalXsc = sigmaTotal * millibarn;
return fTotalXsc;
}
fRadiusEff = 1.34*G4Pow::GetInstance()->powA(A,0.23)+1.35/G4Pow::GetInstance()->powA(A,1./3.); //fm
if ( (Z==1) && (A==2) ) fRadiusEff = 3.800; //fm
if ( (Z==1) && (A==3) ) fRadiusEff = 3.300;
if ( (Z==2) && (A==3) ) fRadiusEff = 3.300;
if ( (Z==2) && (A==4) ) fRadiusEff = 2.376;
// calculation of effective nuclear radius for AntiDeuteron interaction (can be changed)
if (theParticle == theADeuteron) {
fRadiusEff = 1.46 * G4Pow::GetInstance()->powA(A,0.21) + 1.45 / G4Pow::GetInstance()->powA(A,1./3.);
if ( (Z==1) && (A==2) ) fRadiusEff = 3.238; //fm
if ( (Z==1) && (A==3) ) fRadiusEff = 3.144;
if ( (Z==2) && (A==3) ) fRadiusEff = 3.144;
if ( (Z==2) && (A==4) ) fRadiusEff = 2.544;
}
//calculation of effective nuclear radius for AntiAlpha interaction (can be changed)
// calculation of effective nuclear radius for AntiHe3 interaction (can be changed)
if ( (theParticle ==theAHe3) || (theParticle ==theATriton) ) {
fRadiusEff = 1.40* G4Pow::GetInstance()->powA(A,0.21)+1.63/G4Pow::GetInstance()->powA(A,1./3.);
if ( (Z==1) && (A==2) ) fRadiusEff = 3.144; //fm
if ( (Z==1) && (A==3) ) fRadiusEff = 3.075;
if ( (Z==2) && (A==3) ) fRadiusEff = 3.075;
if ( (Z==2) && (A==4) ) fRadiusEff = 2.589;
}
if (theParticle == theAAlpha)
{
fRadiusEff = 1.35* G4Pow::GetInstance()->powA(A,0.21)+1.1/G4Pow::GetInstance()->powA(A,1./3.);
if( (Z==1) && (A==2) ) fRadiusEff = 2.544; //fm
if( (Z==1) && (A==3) ) fRadiusEff = 2.589;
if( (Z==2) && (A==3) ) fRadiusEff = 2.589;
if( (Z==2) && (A==4) ) fRadiusEff = 2.241;
}
// calculation of effective nuclear radius for AntiAlpha interaction (can be changed)
if (theParticle == theAAlpha) {
fRadiusEff = 1.35* G4Pow::GetInstance()->powA(A,0.21)+1.1/G4Pow::GetInstance()->powA(A,1./3.);
if ( (Z==1) && (A==2) ) fRadiusEff = 2.544; //fm
if ( (Z==1) && (A==3) ) fRadiusEff = 2.589;
if ( (Z==2) && (A==3) ) fRadiusEff = 2.589;
if ( (Z==2) && (A==4) ) fRadiusEff = 2.241;
}
G4double R2 = fRadiusEff*fRadiusEff;
G4double REf2 = R2+fRadiusNN2;
G4double ApAt = std::abs(theParticle->GetBaryonNumber()) * A;
G4double R2 = fRadiusEff*fRadiusEff;
G4double REf2 = R2+fRadiusNN2;
G4double ApAt = std::abs(theParticle->GetBaryonNumber()) * A;
xsection = 2*pi*REf2*10.*G4Log(1+(ApAt*sigmaTotal/(2*pi*REf2*10.))); //mb
xsection =xsection *millibarn;
fTotalXsc = xsection;
xsection = 2*pi*REf2*10.*G4Log(1+(ApAt*sigmaTotal/(2*pi*REf2*10.))); //mb
xsection = xsection *millibarn;
fTotalXsc = xsection;
return fTotalXsc;
return fTotalXsc;
}
////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// Calculation of total CrossSection of Anti-Nucleus - Nucleus
//////////////////////////////////////////////////////////////////////////////
G4double G4ComponentAntiNuclNuclearXS::GetTotalIsotopeCrossSection
(const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4int A )
{ return GetTotalElementCrossSection(aParticle, kinEnergy, Z, (G4double) A); }
{
return GetTotalElementCrossSection(aParticle, kinEnergy, Z, (G4double) A);
}
////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
// Calculation of inelastic CrossSection of Anti-Nucleus - Nucleus
////////////////////////////////////////////////////////////////
G4double G4ComponentAntiNuclNuclearXS::GetInelasticElementCrossSection
(const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4double A)
{
G4double inelxsection, sigmaTotal, sigmaElastic;
const G4ParticleDefinition* theParticle = aParticle;
sigmaTotal = GetAntiHadronNucleonTotCrSc(theParticle,kinEnergy);
sigmaElastic = GetAntiHadronNucleonElCrSc(theParticle,kinEnergy);
sigmaTotal = GetAntiHadronNucleonTotCrSc(theParticle,kinEnergy);
sigmaElastic = GetAntiHadronNucleonElCrSc(theParticle,kinEnergy);
// calculation of sqr of radius NN-collision
fRadiusNN2=sigmaTotal*sigmaTotal*0.1/(8.*sigmaElastic*pi); // fm^2
// calculation of effective nuclear radius for Pbar and Nbar interaction (can be changed)
// calculation of sqr of radius NN-collision
fRadiusNN2=sigmaTotal*sigmaTotal*0.1/(8.*sigmaElastic*pi); // fm^2
// calculation of effective nuclear radius for Pbar and Nbar interaction (can be changed)
//A.R. 29-Jan-2013 : use antiprotons/antineutrons as the default case,
// to be used for instance, as first approximation
// without validation, for anti-hyperons.
//if ( (theParticle == theAProton) || (theParticle == theANeutron) ) {
if (A==1)
{ fInelasticXsc = (sigmaTotal - sigmaElastic)*millibarn;
return fInelasticXsc;
}
fRadiusEff = 1.31*G4Pow::GetInstance()->powA(A, 0.22)+0.9/G4Pow::GetInstance()->powA(A, 1./3.); //fm
if( (Z==1) && (A==2) ) fRadiusEff = 3.582; //fm
if( (Z==1) && (A==3) ) fRadiusEff = 3.105;
if( (Z==2) && (A==3) ) fRadiusEff = 3.105;
if( (Z==2) && (A==4) ) fRadiusEff = 2.209;
//}
//calculation of effective nuclear radius for AntiDeuteron interaction (can be changed)
if (theParticle ==theADeuteron)
{
fRadiusEff = 1.38*G4Pow::GetInstance()->powA(A, 0.21)+1.55/G4Pow::GetInstance()->powA(A, 1./3.);
if( (Z==1) && (A==2) ) fRadiusEff = 3.169; //fm
if( (Z==1) && (A==3) ) fRadiusEff = 3.066;
if( (Z==2) && (A==3) ) fRadiusEff = 3.066;
if( (Z==2) && (A==4) ) fRadiusEff = 2.498;
}
if (A==1) {
fInelasticXsc = (sigmaTotal - sigmaElastic) * millibarn;
return fInelasticXsc;
}
fRadiusEff = 1.31*G4Pow::GetInstance()->powA(A, 0.22)+0.9/G4Pow::GetInstance()->powA(A, 1./3.); //fm
if ( (Z==1) && (A==2) ) fRadiusEff = 3.582; //fm
if ( (Z==1) && (A==3) ) fRadiusEff = 3.105;
if ( (Z==2) && (A==3) ) fRadiusEff = 3.105;
if ( (Z==2) && (A==4) ) fRadiusEff = 2.209;
// calculation of effective nuclear radius for AntiDeuteron interaction (can be changed)
if (theParticle ==theADeuteron) {
fRadiusEff = 1.38*G4Pow::GetInstance()->powA(A, 0.21)+1.55/G4Pow::GetInstance()->powA(A, 1./3.);
if ( (Z==1) && (A==2) ) fRadiusEff = 3.169; //fm
if ( (Z==1) && (A==3) ) fRadiusEff = 3.066;
if ( (Z==2) && (A==3) ) fRadiusEff = 3.066;
if ( (Z==2) && (A==4) ) fRadiusEff = 2.498;
}
//calculation of effective nuclear radius for AntiHe3 interaction (can be changed)
// calculation of effective nuclear radius for AntiHe3 interaction (can be changed)
if ( (theParticle ==theAHe3) || (theParticle ==theATriton) ) {
fRadiusEff = 1.34 * G4Pow::GetInstance()->powA(A, 0.21)+1.51/G4Pow::GetInstance()->powA(A, 1./3.);
if ( (Z==1) && (A==2) ) fRadiusEff = 3.066; //fm
if ( (Z==1) && (A==3) ) fRadiusEff = 2.973;
if ( (Z==2) && (A==3) ) fRadiusEff = 2.973;
if ( (Z==2) && (A==4) ) fRadiusEff = 2.508;
}
if( (theParticle ==theAHe3) || (theParticle ==theATriton) )
{
fRadiusEff = 1.34 * G4Pow::GetInstance()->powA(A, 0.21)+1.51/G4Pow::GetInstance()->powA(A, 1./3.);
if( (Z==1) && (A==2) ) fRadiusEff = 3.066; //fm
if( (Z==1) && (A==3) ) fRadiusEff = 2.973;
if( (Z==2) && (A==3) ) fRadiusEff = 2.973;
if( (Z==2) && (A==4) ) fRadiusEff = 2.508;
}
// calculation of effective nuclear radius for AntiAlpha interaction (can be changed)
if (theParticle == theAAlpha) {
fRadiusEff = 1.3*G4Pow::GetInstance()->powA(A, 0.21)+1.05/G4Pow::GetInstance()->powA(A, 1./3.);
if ( (Z==1) && (A==2) ) fRadiusEff = 2.498; //fm
if ( (Z==1) && (A==3) ) fRadiusEff = 2.508;
if ( (Z==2) && (A==3) ) fRadiusEff = 2.508;
if ( (Z==2) && (A==4) ) fRadiusEff = 2.158;
}
//calculation of effective nuclear radius for AntiAlpha interaction (can be changed)
G4double R2 = fRadiusEff*fRadiusEff;
G4double REf2 = R2+fRadiusNN2;
G4double ApAt = std::abs(theParticle->GetBaryonNumber()) * A;
if (theParticle == theAAlpha)
{
fRadiusEff = 1.3*G4Pow::GetInstance()->powA(A, 0.21)+1.05/G4Pow::GetInstance()->powA(A, 1./3.);
if( (Z==1) && (A==2) ) fRadiusEff = 2.498; //fm
if( (Z==1) && (A==3) ) fRadiusEff = 2.508;
if( (Z==2) && (A==3) ) fRadiusEff = 2.508;
if( (Z==2) && (A==4) ) fRadiusEff = 2.158;
}
G4double R2 = fRadiusEff*fRadiusEff;
G4double REf2 = R2+fRadiusNN2;
G4double ApAt= std::abs(theParticle->GetBaryonNumber()) * A;
inelxsection = pi*REf2 *10* G4Log(1+(ApAt*sigmaTotal/(pi*REf2*10.))); //mb
inelxsection = inelxsection * millibarn;
fInelasticXsc = inelxsection;
// scale inelastic cross-section with a momentum-dependent factor
// if (theParticle == theADeuteron || theParticle == theAProton) {
fScalingFactorXsc = GetScalingFactorCrSc(aParticle,kinEnergy);
fInelasticXsc = inelxsection*fScalingFactorXsc;
// }
inelxsection = pi*REf2 *10* G4Log(1+(ApAt*sigmaTotal/(pi*REf2*10.))); //mb
inelxsection = inelxsection * millibarn;
fInelasticXsc = inelxsection;
return fInelasticXsc;
// scale inelastic cross-section with a momentum-dependent factor
fScalingFactorXsc = GetScalingFactorCrSc(aParticle,kinEnergy);
fInelasticXsc = inelxsection*fScalingFactorXsc;
return fInelasticXsc;
}
///////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// Calculates Inelastic Anti-nucleus-Nucleus cross-section
//
G4double G4ComponentAntiNuclNuclearXS::GetInelasticIsotopeCrossSection
(const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4int A)
{return GetInelasticElementCrossSection(aParticle, kinEnergy, Z, (G4double) A); }
{
return GetInelasticElementCrossSection(aParticle, kinEnergy, Z, (G4double) A);
}
///////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// Calculates elastic Anti-nucleus-Nucleus cross-section as Total - Inelastic
//
G4double G4ComponentAntiNuclNuclearXS::GetElasticElementCrossSection
(const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4double A)
{
fElasticXsc = GetTotalElementCrossSection(aParticle, kinEnergy, Z, A)-
GetInelasticElementCrossSection(aParticle, kinEnergy, Z, A);
// keep elastic cross-section unchanged
// if (aParticle == theADeuteron || aParticle == theAProton) {
fScalingFactorXsc = GetScalingFactorCrSc(aParticle,kinEnergy);
fElasticXsc = GetTotalElementCrossSection(aParticle, kinEnergy, Z, A)-
GetInelasticElementCrossSection(aParticle, kinEnergy, Z, A)/fScalingFactorXsc;
// }
// keep elastic cross-section unchanged
fScalingFactorXsc = GetScalingFactorCrSc(aParticle,kinEnergy);
if (fElasticXsc < 0.) fElasticXsc = 0.;
fElasticXsc = GetTotalElementCrossSection(aParticle, kinEnergy, Z, A)-
GetInelasticElementCrossSection(aParticle, kinEnergy, Z, A)/fScalingFactorXsc;
return fElasticXsc;
if (fElasticXsc < 0.) fElasticXsc = 0.;
return fElasticXsc;
}
///////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// Calculates elastic Anti-nucleus-Nucleus cross-section
//
G4double G4ComponentAntiNuclNuclearXS::GetElasticIsotopeCrossSection
(const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4int A)
{ return GetElasticElementCrossSection(aParticle, kinEnergy, Z, (G4double) A); }
{
return GetElasticElementCrossSection(aParticle, kinEnergy, Z, (G4double) A);
}
///////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
// Calculation of Antihadron - hadron Total Cross-section
G4double G4ComponentAntiNuclNuclearXS::GetAntiHadronNucleonTotCrSc
......@@ -317,97 +279,61 @@ G4double G4ComponentAntiNuclNuclearXS::GetAntiHadronNucleonTotCrSc
momentum=std::sqrt(Energy*Energy-Pmass*Pmass)/std::abs(theParticle->GetBaryonNumber());
G4double Plab = momentum / GeV;
G4double B, SigAss;
G4double C, d1, d2, d3 ;
Elab = std::sqrt(Mn*Mn + Plab*Plab); // GeV
S = 2.*Mn*Mn + 2. *Mn*Elab; // GeV^2
SqrtS = std::sqrt(S); // GeV
B = b0+b2*G4Log(SqrtS/SqrtS0)*G4Log(SqrtS/SqrtS0); //GeV^(-2)
SigAss = 36.04 +0.304*G4Log(S/S0)*G4Log(S/S0); //mb
R0 = std::sqrt(0.40874044*SigAss - B); //GeV^(-2)
C = 13.55;
d1 = -4.47;
d2 = 12.38;
d3 = -12.43;
xsection = SigAss*(1 + 1./(std::sqrt(S-4.*Mn*Mn)) / (G4Pow::GetInstance()->powA(R0, 3.))
*C* (1+d1/SqrtS+d2/(G4Pow::GetInstance()->powA(SqrtS,2.))+d3/(G4Pow::GetInstance()->powA(SqrtS,3.)) ));
// xsection *= millibarn;
G4double B, SigAss;
G4double C, d1, d2, d3;
Elab = std::sqrt(Mn*Mn + Plab*Plab); // GeV
S = 2.*Mn*Mn + 2. *Mn*Elab; // GeV^2
SqrtS = std::sqrt(S); // GeV
B = b0+b2*G4Log(SqrtS/SqrtS0)*G4Log(SqrtS/SqrtS0); //GeV^(-2)
SigAss = 36.04 +0.304*G4Log(S/S0)*G4Log(S/S0); //mb
R0 = std::sqrt(0.40874044*SigAss - B); //GeV^(-2)
C = 13.55;
d1 = -4.47;
d2 = 12.38;
d3 = -12.43;
xsection = SigAss * ( 1 + 1./(std::sqrt(S-4.*Mn*Mn)) / (G4Pow::GetInstance()->powA(R0, 3.))
* C * ( 1 + d1/SqrtS + d2/(G4Pow::GetInstance()->powA(SqrtS,2.))
+ d3/(G4Pow::GetInstance()->powA(SqrtS,3.)) ) );
//xsection *= millibarn;
fAntiHadronNucleonTotXsc = xsection;
return fAntiHadronNucleonTotXsc;
}
//
// /////////////////////////////////////////////////////////////////////////////////
// //////////////////////////////////////////////////////////////////////////
// Calculation of Antihadron - hadron Elastic Cross-section
G4double G4ComponentAntiNuclNuclearXS ::
GetAntiHadronNucleonElCrSc(const G4ParticleDefinition* aParticle, G4double kinEnergy)
{
G4double xsection;
G4double SigAss;
G4double C, d1, d2, d3 ;
GetAntiHadronNucleonTotCrSc(aParticle,kinEnergy);
SigAss = 4.5 + 0.101*G4Log(S/S0)*G4Log(S/S0); //mb
C = 59.27;
d1 = -6.95;
d2 = 23.54;
d3 = -25.34;
xsection = SigAss* (1 + 1. / (std::sqrt(S-4.*Mn*Mn)) / (G4Pow::GetInstance()->powA(R0, 3.))
*C* ( 1+d1/SqrtS+d2/(G4Pow::GetInstance()->powA(SqrtS,2.))+d3/(G4Pow::GetInstance()->powA(SqrtS,3.)) ));
// xsection *= millibarn;
G4double xsection;
G4double SigAss;
G4double C, d1, d2, d3;
GetAntiHadronNucleonTotCrSc(aParticle,kinEnergy);
SigAss = 4.5 + 0.101*G4Log(S/S0)*G4Log(S/S0); //mb
C = 59.27;
d1 = -6.95;
d2 = 23.54;
d3 = -25.34;
xsection = SigAss * ( 1 + 1. / (std::sqrt(S-4.*Mn*Mn)) / (G4Pow::GetInstance()->powA(R0, 3.))
* C * ( 1 + d1/SqrtS + d2/(G4Pow::GetInstance()->powA(SqrtS,2.))
+ d3/(G4Pow::GetInstance()->powA(SqrtS,3.)) ) );
//xsection *= millibarn;
fAntiHadronNucleonElXsc = xsection;
return fAntiHadronNucleonElXsc;
}
////////////////////////////////////////////////////////////////////////////
// Calculation of scaling factor for Cross-section
G4double G4ComponentAntiNuclNuclearXS::GetScalingFactorCrSc(const G4ParticleDefinition* aParticle, G4double kinEnergy){
// const G4ParticleDefinition* theParticle = aParticle;
// G4double mass = theParticle->GetPDGMass();
// G4double p = std::sqrt(kinEnergy*kinEnergy + 2.0*kinEnergy*mass); // p in MeV
// p = p/1000.; // make p in GeV
//
// static void* lib = dlopen("./CustomScalingFunction.so", RTLD_NOW);
// static void* fun_ad = dlsym(lib, "CustomScalingFunctionAntiDeuteron");
// static void* fun_ap = dlsym(lib, "CustomScalingFunctionAntiProton");
//
// if (!lib){
// return 1.0;
// }
//
// if ((theParticle == theADeuteron && !fun_ad) ||
// (theParticle == theAProton && !fun_ap)){
// return 1.0;
// }
//
// static double (*scaling_ap)(double) = (double (*)(double)) fun_ap;
// static double (*scaling_ad)(double) = (double (*)(double)) fun_ad;
//
// if (theParticle == theADeuteron){
// return scaling_ad(p);
// }
//
// if (theParticle == theAProton){
// return scaling_ap(p);
// }
//
// return 1.0;
// generic implementation from Max
const G4ParticleDefinition* theParticle = aParticle;
G4double mass = theParticle->GetPDGMass();
G4double p = std::sqrt(kinEnergy*kinEnergy + 2.0*kinEnergy*mass); // p in MeV
......@@ -429,13 +355,15 @@ G4double G4ComponentAntiNuclNuclearXS::GetScalingFactorCrSc(const G4ParticleDefi
}
/////////////////////////////////////////////////////////////////////////////
void G4ComponentAntiNuclNuclearXS::CrossSectionDescription(std::ostream& outFile) const
{
outFile << "The G4ComponentAntiNuclNuclearXS calculates total,\n"
<< "inelastic, elastic cross sections of anti-nucleons and light \n"
<< "anti-nucleus interactions with nuclei using Glauber's approach.\n"
<< "anti-nucleus interactions with nuclei using Glauber's approach.\n"
<< "It uses parametrizations of antiproton-proton total and elastic \n"
<< "cross sections and Wood-Saxon distribution of nuclear density.\n"
<< "The lower limit is 10 MeV, the upper limit is 10 TeV. \n"
<< "cross sections and Wood-Saxon distribution of nuclear density.\n"
<< "See details in Phys.Lett. B705 (2011) 235. \n";
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment