ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ComponentAntiNuclNuclearXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ComponentAntiNuclNuclearXS.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Calculation of the total, elastic and inelastic cross-sections
27 // of anti-nucleon and anti-nucleus interactions with nuclei
28 // based on Glauber approach and V. Grishine formulaes for
29 // interpolations (ref. V.M.Grichine, Eur.Phys.J., C62(2009) 399;
30 // NIM, B267 (2009) 2460) and our parametrization of hadron-nucleon
31 // cross-sections
32 //
33 //
34 // Created by A.Galoyan and V. Uzhinsky, 18.11.2010
35 
36 
37 #include "G4ComponentAntiNuclNuclearXS.hh"
38 
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4ParticleTable.hh"
42 #include "G4IonTable.hh"
43 #include "G4ParticleDefinition.hh"
44 #include "G4Pow.hh"
45 
46 
48 
49 G4ComponentAntiNuclNuclearXS::G4ComponentAntiNuclNuclearXS()
50 : G4VComponentCrossSection("AntiAGlauber"),
51  fRadiusEff(0.0), fRadiusNN2(0.0),
52  fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0),
53  fAntiHadronNucleonTotXsc(0.0), fAntiHadronNucleonElXsc(0.0),
54  Elab(0.0), S(0.0), SqrtS(0)
55 {
56  theAProton = G4AntiProton::AntiProton();
57  theANeutron = G4AntiNeutron::AntiNeutron();
58  theADeuteron = G4AntiDeuteron::AntiDeuteron();
59  theATriton = G4AntiTriton::AntiTriton();
60  theAAlpha = G4AntiAlpha::AntiAlpha();
61  theAHe3 = G4AntiHe3::AntiHe3();
62  Mn = 0.93827231; // GeV
63  b0 = 11.92; // GeV^(-2)
64  b2 = 0.3036; // GeV^(-2)
65  SqrtS0 = 20.74; // GeV
66  S0 = 33.0625; // GeV^2
67  R0 = 1.0; // default value (V.Ivanchenko)
68 }
69 
70 
72 
73 G4ComponentAntiNuclNuclearXS::~G4ComponentAntiNuclNuclearXS()
74 {
75 }
76 
77 
79 //
80 // Calculation of total CrossSection of Anti-Nucleus - Nucleus
81 
82 G4double G4ComponentAntiNuclNuclearXS::GetTotalElementCrossSection
83 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4double A)
84 {
85  G4double xsection, sigmaTotal, sigmaElastic;
86  const G4ParticleDefinition* theParticle = aParticle;
87  sigmaTotal = GetAntiHadronNucleonTotCrSc(theParticle,kinEnergy);
88  sigmaElastic = GetAntiHadronNucleonElCrSc(theParticle,kinEnergy);
89 
90  // calculation of squared radius of NN-collision
91  fRadiusNN2=sigmaTotal*sigmaTotal*0.1/(8.*sigmaElastic*pi) ; //fm^2
92 
93  // calculation of effective nuclear radius for Pbar and Nbar interactions (can be changed)
94  //A.R. 29-Jan-2013 : use antiprotons/antineutrons as the default case,
95  // to be used for instance, as first approximation
96  // without validation, for anti-hyperons.
97  if (A==1) {
98  fTotalXsc = sigmaTotal * millibarn;
99  return fTotalXsc;
100  }
101  fRadiusEff = 1.34*G4Pow::GetInstance()->powA(A,0.23)+1.35/G4Pow::GetInstance()->powA(A,1./3.); //fm
102  if ( (Z==1) && (A==2) ) fRadiusEff = 3.800; //fm
103  if ( (Z==1) && (A==3) ) fRadiusEff = 3.300;
104  if ( (Z==2) && (A==3) ) fRadiusEff = 3.300;
105  if ( (Z==2) && (A==4) ) fRadiusEff = 2.376;
106 
107  // calculation of effective nuclear radius for AntiDeuteron interaction (can be changed)
108  if (theParticle == theADeuteron) {
109  fRadiusEff = 1.46 * G4Pow::GetInstance()->powA(A,0.21) + 1.45 / G4Pow::GetInstance()->powA(A,1./3.);
110  if ( (Z==1) && (A==2) ) fRadiusEff = 3.238; //fm
111  if ( (Z==1) && (A==3) ) fRadiusEff = 3.144;
112  if ( (Z==2) && (A==3) ) fRadiusEff = 3.144;
113  if ( (Z==2) && (A==4) ) fRadiusEff = 2.544;
114  }
115 
116  // calculation of effective nuclear radius for AntiHe3 interaction (can be changed)
117  if ( (theParticle ==theAHe3) || (theParticle ==theATriton) ) {
118  fRadiusEff = 1.40* G4Pow::GetInstance()->powA(A,0.21)+1.63/G4Pow::GetInstance()->powA(A,1./3.);
119  if ( (Z==1) && (A==2) ) fRadiusEff = 3.144; //fm
120  if ( (Z==1) && (A==3) ) fRadiusEff = 3.075;
121  if ( (Z==2) && (A==3) ) fRadiusEff = 3.075;
122  if ( (Z==2) && (A==4) ) fRadiusEff = 2.589;
123  }
124 
125  // calculation of effective nuclear radius for AntiAlpha interaction (can be changed)
126  if (theParticle == theAAlpha) {
127  fRadiusEff = 1.35* G4Pow::GetInstance()->powA(A,0.21)+1.1/G4Pow::GetInstance()->powA(A,1./3.);
128  if ( (Z==1) && (A==2) ) fRadiusEff = 2.544; //fm
129  if ( (Z==1) && (A==3) ) fRadiusEff = 2.589;
130  if ( (Z==2) && (A==3) ) fRadiusEff = 2.589;
131  if ( (Z==2) && (A==4) ) fRadiusEff = 2.241;
132  }
133 
134  G4double R2 = fRadiusEff*fRadiusEff;
135  G4double REf2 = R2+fRadiusNN2;
136  G4double ApAt = std::abs(theParticle->GetBaryonNumber()) * A;
137 
138  xsection = 2*pi*REf2*10.*G4Log(1+(ApAt*sigmaTotal/(2*pi*REf2*10.))); //mb
139  xsection = xsection *millibarn;
140  fTotalXsc = xsection;
141 
142  return fTotalXsc;
143 }
144 
145 
147 //
148 // Calculation of total CrossSection of Anti-Nucleus - Nucleus
149 
150 G4double G4ComponentAntiNuclNuclearXS::GetTotalIsotopeCrossSection
151 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4int A )
152 {
153  return GetTotalElementCrossSection(aParticle, kinEnergy, Z, (G4double) A);
154 }
155 
156 
158 // Calculation of inelastic CrossSection of Anti-Nucleus - Nucleus
159 
160 G4double G4ComponentAntiNuclNuclearXS::GetInelasticElementCrossSection
161 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4double A)
162 {
163  G4double inelxsection, sigmaTotal, sigmaElastic;
164  const G4ParticleDefinition* theParticle = aParticle;
165  sigmaTotal = GetAntiHadronNucleonTotCrSc(theParticle,kinEnergy);
166  sigmaElastic = GetAntiHadronNucleonElCrSc(theParticle,kinEnergy);
167 
168  // calculation of sqr of radius NN-collision
169  fRadiusNN2=sigmaTotal*sigmaTotal*0.1/(8.*sigmaElastic*pi); // fm^2
170 
171  // calculation of effective nuclear radius for Pbar and Nbar interaction (can be changed)
172  //A.R. 29-Jan-2013 : use antiprotons/antineutrons as the default case,
173  // to be used for instance, as first approximation
174  // without validation, for anti-hyperons.
175  if (A==1) {
176  fInelasticXsc = (sigmaTotal - sigmaElastic) * millibarn;
177  return fInelasticXsc;
178  }
179  fRadiusEff = 1.31*G4Pow::GetInstance()->powA(A, 0.22)+0.9/G4Pow::GetInstance()->powA(A, 1./3.); //fm
180  if ( (Z==1) && (A==2) ) fRadiusEff = 3.582; //fm
181  if ( (Z==1) && (A==3) ) fRadiusEff = 3.105;
182  if ( (Z==2) && (A==3) ) fRadiusEff = 3.105;
183  if ( (Z==2) && (A==4) ) fRadiusEff = 2.209;
184 
185  // calculation of effective nuclear radius for AntiDeuteron interaction (can be changed)
186  if (theParticle ==theADeuteron) {
187  fRadiusEff = 1.38*G4Pow::GetInstance()->powA(A, 0.21)+1.55/G4Pow::GetInstance()->powA(A, 1./3.);
188  if ( (Z==1) && (A==2) ) fRadiusEff = 3.169; //fm
189  if ( (Z==1) && (A==3) ) fRadiusEff = 3.066;
190  if ( (Z==2) && (A==3) ) fRadiusEff = 3.066;
191  if ( (Z==2) && (A==4) ) fRadiusEff = 2.498;
192  }
193 
194  // calculation of effective nuclear radius for AntiHe3 interaction (can be changed)
195  if ( (theParticle ==theAHe3) || (theParticle ==theATriton) ) {
196  fRadiusEff = 1.34 * G4Pow::GetInstance()->powA(A, 0.21)+1.51/G4Pow::GetInstance()->powA(A, 1./3.);
197  if ( (Z==1) && (A==2) ) fRadiusEff = 3.066; //fm
198  if ( (Z==1) && (A==3) ) fRadiusEff = 2.973;
199  if ( (Z==2) && (A==3) ) fRadiusEff = 2.973;
200  if ( (Z==2) && (A==4) ) fRadiusEff = 2.508;
201  }
202 
203  // calculation of effective nuclear radius for AntiAlpha interaction (can be changed)
204  if (theParticle == theAAlpha) {
205  fRadiusEff = 1.3*G4Pow::GetInstance()->powA(A, 0.21)+1.05/G4Pow::GetInstance()->powA(A, 1./3.);
206  if ( (Z==1) && (A==2) ) fRadiusEff = 2.498; //fm
207  if ( (Z==1) && (A==3) ) fRadiusEff = 2.508;
208  if ( (Z==2) && (A==3) ) fRadiusEff = 2.508;
209  if ( (Z==2) && (A==4) ) fRadiusEff = 2.158;
210  }
211 
212  G4double R2 = fRadiusEff*fRadiusEff;
213  G4double REf2 = R2+fRadiusNN2;
214  G4double ApAt = std::abs(theParticle->GetBaryonNumber()) * A;
215 
216  inelxsection = pi*REf2 *10* G4Log(1+(ApAt*sigmaTotal/(pi*REf2*10.))); //mb
217  inelxsection = inelxsection * millibarn;
218  fInelasticXsc = inelxsection;
219 
220  return fInelasticXsc;
221 }
222 
223 
225 //
226 // Calculates Inelastic Anti-nucleus-Nucleus cross-section
227 
228 G4double G4ComponentAntiNuclNuclearXS::GetInelasticIsotopeCrossSection
229 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4int A)
230 {
231  return GetInelasticElementCrossSection(aParticle, kinEnergy, Z, (G4double) A);
232 }
233 
234 
236 //
237 // Calculates elastic Anti-nucleus-Nucleus cross-section as Total - Inelastic
238 
239 G4double G4ComponentAntiNuclNuclearXS::GetElasticElementCrossSection
240 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4double A)
241 {
242  fElasticXsc = GetTotalElementCrossSection(aParticle, kinEnergy, Z, A)-
243  GetInelasticElementCrossSection(aParticle, kinEnergy, Z, A);
244  if (fElasticXsc < 0.) fElasticXsc = 0.;
245  return fElasticXsc;
246 }
247 
248 
250 //
251 // Calculates elastic Anti-nucleus-Nucleus cross-section
252 
253 G4double G4ComponentAntiNuclNuclearXS::GetElasticIsotopeCrossSection
254 (const G4ParticleDefinition* aParticle, G4double kinEnergy, G4int Z, G4int A)
255 {
256  return GetElasticElementCrossSection(aParticle, kinEnergy, Z, (G4double) A);
257 }
258 
259 
261 // Calculation of Antihadron - hadron Total Cross-section
262 
263 G4double G4ComponentAntiNuclNuclearXS::GetAntiHadronNucleonTotCrSc
264 (const G4ParticleDefinition* aParticle, G4double kinEnergy)
265 {
266  G4double xsection, Pmass, Energy, momentum;
267  const G4ParticleDefinition* theParticle = aParticle;
268  Pmass=theParticle->GetPDGMass();
269  Energy=Pmass+kinEnergy;
270  momentum=std::sqrt(Energy*Energy-Pmass*Pmass)/std::abs(theParticle->GetBaryonNumber());
271  G4double Plab = momentum / GeV;
272 
273  G4double B, SigAss;
274  G4double C, d1, d2, d3;
275  Elab = std::sqrt(Mn*Mn + Plab*Plab); // GeV
276  S = 2.*Mn*Mn + 2. *Mn*Elab; // GeV^2
277  SqrtS = std::sqrt(S); // GeV
278  B = b0+b2*G4Log(SqrtS/SqrtS0)*G4Log(SqrtS/SqrtS0); //GeV^(-2)
279  SigAss = 36.04 +0.304*G4Log(S/S0)*G4Log(S/S0); //mb
280  R0 = std::sqrt(0.40874044*SigAss - B); //GeV^(-2)
281  C = 13.55;
282  d1 = -4.47;
283  d2 = 12.38;
284  d3 = -12.43;
285 
286  xsection = SigAss * ( 1 + 1./(std::sqrt(S-4.*Mn*Mn)) / (G4Pow::GetInstance()->powA(R0, 3.))
287  * C * ( 1 + d1/SqrtS + d2/(G4Pow::GetInstance()->powA(SqrtS,2.))
288  + d3/(G4Pow::GetInstance()->powA(SqrtS,3.)) ) );
289 
290  //xsection *= millibarn;
291  fAntiHadronNucleonTotXsc = xsection;
292 
293  return fAntiHadronNucleonTotXsc;
294 }
295 
296 
297 // //////////////////////////////////////////////////////////////////////////
298 // Calculation of Antihadron - hadron Elastic Cross-section
299 
300 G4double G4ComponentAntiNuclNuclearXS ::
301 GetAntiHadronNucleonElCrSc(const G4ParticleDefinition* aParticle, G4double kinEnergy)
302 {
303  G4double xsection;
304  G4double SigAss;
305  G4double C, d1, d2, d3;
306  GetAntiHadronNucleonTotCrSc(aParticle,kinEnergy);
307  SigAss = 4.5 + 0.101*G4Log(S/S0)*G4Log(S/S0); //mb
308  C = 59.27;
309  d1 = -6.95;
310  d2 = 23.54;
311  d3 = -25.34;
312 
313  xsection = SigAss * ( 1 + 1. / (std::sqrt(S-4.*Mn*Mn)) / (G4Pow::GetInstance()->powA(R0, 3.))
314  * C * ( 1 + d1/SqrtS + d2/(G4Pow::GetInstance()->powA(SqrtS,2.))
315  + d3/(G4Pow::GetInstance()->powA(SqrtS,3.)) ) );
316 
317  //xsection *= millibarn;
318  fAntiHadronNucleonElXsc = xsection;
319 
320  return fAntiHadronNucleonElXsc;
321 }
322 
323 
325 
326 void G4ComponentAntiNuclNuclearXS::CrossSectionDescription(std::ostream& outFile) const
327 {
328  outFile << "The G4ComponentAntiNuclNuclearXS calculates total,\n"
329  << "inelastic, elastic cross sections of anti-nucleons and light \n"
330  << "anti-nucleus interactions with nuclei using Glauber's approach.\n"
331  << "It uses parametrizations of antiproton-proton total and elastic \n"
332  << "cross sections and Wood-Saxon distribution of nuclear density.\n"
333  << "See details in Phys.Lett. B705 (2011) 235. \n";
334 }
335