| Geant4 Cross Reference |
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 // -------------------------------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4ChargeExchangeXS
32 //
33
34 #include "G4ChargeExchangeXS.hh"
35 #include "G4DynamicParticle.hh"
36 #include "G4ElementTable.hh"
37 #include "G4Material.hh"
38 #include "G4Element.hh"
39 #include "G4IsotopeList.hh"
40 #include "G4HadronicParameters.hh"
41 #include "Randomize.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4NucleiProperties.hh"
44 #include "G4Pow.hh"
45
46 #include "G4PionZero.hh"
47 #include "G4Eta.hh"
48 #include "G4KaonZeroLong.hh"
49 #include "G4KaonZeroShort.hh"
50 #include "G4KaonPlus.hh"
51 #include "G4KaonMinus.hh"
52 #include "G4ParticleTable.hh"
53
54 namespace {
55 // V. Lyubovitsky parameterisation
56 const G4double piA[5] = {430., 36., 1.37, 2.0, 60.}; // A
57 const G4double pAP[5] = {1.04, 1.26, 1.35, 0.94, 0.94}; // 2 - 2alphaP
58 const G4double pC0[5] = {12.7, 6.0, 6.84, 6.5, 8.0}; // c0
59 const G4double pC1[5] = {1.57, 1.6, 1.7, 1.23, 2.6}; // c1
60 const G4double pG0[5] = {2.55, 4.6, 3.7, 5.5, 4.6}; // g0
61 const G4double pG1[5] = {-0.23, -0.5, 0., 0., -2.}; // g1
62
63 const G4double beta_prime_pi = 0.0410;
64 }
65
66 G4ChargeExchangeXS::G4ChargeExchangeXS()
67 {
68 if (verboseLevel > 1) {
69 G4cout << "G4ChargeExchangeXS::G4ChargeExchangeXS" << G4endl;
70 }
71 g4calc = G4Pow::GetInstance();
72 auto table = G4ParticleTable::GetParticleTable();
73 const G4String nam[5] = {"pi0", "eta", "eta_prime", "omega", "f2(1270)"};
74 for (G4int i=0; i<5; ++i) {
75 fPionSecPD[i] = table->FindParticle(nam[i]);
76 if (nullptr == fPionSecPD[i]) {
77 G4ExceptionDescription ed;
78 ed << "### meson " << nam[i] << " is not found out in the particle table";
79 G4Exception("G4ChargeExchangeXS::G4ChargeExchangeXS()","had044",
80 FatalException, ed,"");
81 }
82 }
83 }
84
85 // Print the information of this .cc file
86 void G4ChargeExchangeXS::CrossSectionDescription(std::ostream& outFile) const
87 {
88 outFile << "G4ChargeExchangeXS calculates charge exchange cross section for "
89 << "pi+, pi-, K+, K-, KL\n";
90 }
91
92 G4bool G4ChargeExchangeXS::IsElementApplicable(const G4DynamicParticle*,
93 G4int, const G4Material*)
94 {
95 return true;
96 }
97
98 G4double
99 G4ChargeExchangeXS::GetElementCrossSection(const G4DynamicParticle* aParticle,
100 G4int ZZ, const G4Material* mat)
101 {
102 G4double result = 0.0;
103 const G4double pE = aParticle->GetTotalEnergy();
104 if (pE <= fEnergyLimit) { return result; }
105
106 auto part = aParticle->GetDefinition();
107 G4int pdg = part->GetPDGEncoding();
108
109 // Get or calculate the proton mass, particle mass, and s(Lorentz invariant)
110 G4double tM = CLHEP::proton_mass_c2;
111 G4double pM = part->GetPDGMass();
112 G4double lorentz_s = tM*tM + 2*tM*pE + pM*pM;
113 if (lorentz_s <= (tM + pM)*(tM + pM)) { return result; }
114
115 const G4int Z = std::min(ZZ, ZMAXNUCLEARDATA);
116 const G4int A = G4lrint(aeff[Z]);
117
118 if (verboseLevel > 1) {
119 G4cout << "### G4ChargeExchangeXS: " << part->GetParticleName()
120 << " Z=" << Z << " A=" << A << " Etot(GeV)=" << pE/CLHEP::GeV
121 << " s(GeV^2)=" << lorentz_s/(CLHEP::GeV*CLHEP::GeV) << G4endl;
122 }
123
124 // For unit conversion
125 const G4double inv1e7 = 0.1/(CLHEP::GeV*CLHEP::GeV);
126 const G4double fact = 1e-30*CLHEP::cm2;
127 const G4double pfact = 0.1/CLHEP::GeV;
128 const G4double kfact = 56.3*fact;
129 const G4double csmax = 1e-16;
130
131 // The approximation of Glauber-Gribov formula -> extend it from interaction with
132 // proton to nuclei Z^(2/3). The factor g4calc->powA(A,-beta_prime_pi*G4Log(A))
133 // takes into account absorption of pi0 and eta
134
135 // pi- + p -> n + meson (0- pi0, 1- eta, 2- eta', 3- omega, 4- f2(1270))
136 if (pdg == -211) {
137 G4double z23 = g4calc->Z23(Z);
138 G4double x = lorentz_s*inv1e7;
139 G4double logX = G4Log(x);
140 G4double logA = g4calc->logZ(A);
141 G4double xf = g4calc->powZ(A, -beta_prime_pi*logA);
142 G4double sum = 0.0;
143 for (G4int i=0; i<5; ++i) {
144 G4double xg = std::max(1.0 + pG0[i] + pG1[i]*logX, 0.0);
145 G4double xc = std::max(pC0[i] + pC1[i]*logX, csmax);
146 G4double xs = z23*piA[i]*g4calc->powA(x, -pAP[i])*xf*xg/xc;
147 sum += xs;
148 fXSecPion[i] = sum;
149 }
150 result = sum*fact;
151 }
152
153 // pi+ + n -> p + meson (0- pi0, 1- eta, 2- eta', 3- omega, 4- f2(1270))
154 else if (pdg == 211) {
155 G4double n23 = g4calc->Z23(A - Z);
156 G4double x = lorentz_s*inv1e7;
157 G4double logX = G4Log(x);
158 G4double logA = g4calc->logZ(A);
159 G4double xf = g4calc->powZ(A, -beta_prime_pi*logA);
160
161 // hydrogen target case Z = A = 1
162 // the cross section is defined by fraction of deuteron and tritium
163 if (1 == Z) { n23 = ComputeDeuteronFraction(mat); }
164 G4double sum = 0.0;
165 for (G4int i=0; i<5; ++i) {
166 G4double xg = std::max(1.0 + pG0[i] + pG1[i]*logX, 0.0);
167 G4double xc = std::max(pC0[i] + pC1[i]*logX, csmax);
168 G4double xs = n23*piA[i]*g4calc->powA(x, -pAP[i])*xf*xg/xc;
169 sum += xs;
170 fXSecPion[i] = sum;
171 }
172 result = sum*fact;
173 }
174
175 // Kaon x-sections depend on the primary particles momentum
176 // K- + p -> Kbar + n
177 else if (pdg == -321) {
178 G4double p_momentum = std::sqrt(pE*pE - pM*pM)*pfact;
179 result = g4calc->Z23(Z)*g4calc->powA(p_momentum, -1.60)*kfact;
180 }
181
182 // K+ + n -> Kbar + p
183 else if (pdg == 321) {
184 G4double p_momentum = std::sqrt(pE*pE - pM*pM)*pfact;
185 G4double n23 = g4calc->Z23(A-Z);
186 // hydrogen target case Z = A = 1
187 // the cross section is defined by fraction of deuteron and tritium
188 if (1 == Z) { n23 = ComputeDeuteronFraction(mat); }
189 result = n23*g4calc->powA(p_momentum, -1.60)*kfact;
190 }
191
192 // KL
193 else if (pdg == 130) {
194 // Cross section of KL = 0.5*(Cross section of K+ + Cross section of K-)
195 const G4double p_momentum = std::sqrt(pE*pE - pM*pM)*pfact;
196 result = 0.5*(g4calc->Z23(Z) + g4calc->Z23(A-Z))*
197 g4calc->powA(p_momentum, -1.60)*kfact;
198 }
199 result *= fFactor;
200 if (verboseLevel > 1) {
201 G4cout << " Done for " << part->GetParticleName() << " Etot(GeV)=" << pE/CLHEP::GeV
202 << " res(mb)=" << result/CLHEP::millibarn << G4endl;
203 }
204 return result;
205 }
206
207 const G4ParticleDefinition*
208 G4ChargeExchangeXS::SampleSecondaryType(const G4ParticleDefinition* part,
209 const G4int Z, const G4int A)
210 {
211 const G4ParticleDefinition* pd = nullptr;
212 G4int pdg = std::abs(part->GetPDGEncoding());
213
214 // pi- + p / pi+ + n
215 if (pdg == 211) {
216 pd = fPionSecPD[0];
217 G4double x = fXSecPion[4]*G4UniformRand();
218 for (G4int i=0; i<5; ++i) {
219 if (x <= fXSecPion[i]) {
220 pd = fPionSecPD[i];
221 break;
222 }
223 }
224 }
225
226 // K- + p / K+ + n
227 // Equal opportunity of producing k-short and k-long
228 else if (pdg == 321) {
229 if (G4UniformRand() >= 0.5) {
230 pd = G4KaonZeroLong::KaonZeroLong();
231 }
232 else {
233 pd = G4KaonZeroShort::KaonZeroShort();
234 }
235 }
236
237 // KL + nucleus
238 else if (pdg == 130) {
239 G4double prob = (G4double)Z/(G4double)A;
240 if (G4UniformRand() >= prob) {
241 pd = G4KaonMinus::KaonMinus();
242 }
243 else {
244 pd = G4KaonPlus::KaonPlus();
245 }
246 }
247
248 return pd;
249 }
250
251 G4double
252 G4ChargeExchangeXS::ComputeDeuteronFraction(const G4Material* mat)
253 {
254 for (auto const & elm : *mat->GetElementVector()) {
255 if (1 == elm->GetZasInt()) {
256 G4double ab = 0.0;
257 const G4int nIso = (G4int)elm->GetNumberOfIsotopes();
258 const G4double* abu = elm->GetRelativeAbundanceVector();
259 for (G4int j = 0; j < nIso; ++j) {
260 auto const iso = elm->GetIsotope(j);
261 ab += (iso->GetN() - iso->GetZ())*abu[j];
262 }
263 return ab;
264 }
265 }
266 return 0.0;
267 }
268