| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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,
57 const G4double pAP[5] = {1.04, 1.26, 1.35,
58 const G4double pC0[5] = {12.7, 6.0, 6.84,
59 const G4double pC1[5] = {1.57, 1.6, 1.7,
60 const G4double pG0[5] = {2.55, 4.6, 3.7,
61 const G4double pG1[5] = {-0.23, -0.5, 0.,
62
63 const G4double beta_prime_pi = 0.0410;
64 }
65
66 G4ChargeExchangeXS::G4ChargeExchangeXS()
67 {
68 if (verboseLevel > 1) {
69 G4cout << "G4ChargeExchangeXS::G4ChargeEx
70 }
71 g4calc = G4Pow::GetInstance();
72 auto table = G4ParticleTable::GetParticleTab
73 const G4String nam[5] = {"pi0", "eta", "eta_
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
79 G4Exception("G4ChargeExchangeXS::G4Charg
80 FatalException, ed,"");
81 }
82 }
83 }
84
85 // Print the information of this .cc file
86 void G4ChargeExchangeXS::CrossSectionDescripti
87 {
88 outFile << "G4ChargeExchangeXS calculates ch
89 << "pi+, pi-, K+, K-, KL\n";
90 }
91
92 G4bool G4ChargeExchangeXS::IsElementApplicable
93
94 {
95 return true;
96 }
97
98 G4double
99 G4ChargeExchangeXS::GetElementCrossSection(con
100 G4int ZZ, const G4Material*
101 {
102 G4double result = 0.0;
103 const G4double pE = aParticle->GetTotalEnerg
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, particl
110 G4double tM = CLHEP::proton_mass_c2;
111 G4double pM = part->GetPDGMass();
112 G4double lorentz_s = tM*tM + 2*tM*pE + pM*p
113 if (lorentz_s <= (tM + pM)*(tM + pM)) { retu
114
115 const G4int Z = std::min(ZZ, ZMAXNUCLEARDATA
116 const G4int A = G4lrint(aeff[Z]);
117
118 if (verboseLevel > 1) {
119 G4cout << "### G4ChargeExchangeXS: " << pa
120 << " Z=" << Z << " A=" << A << " Etot(GeV
121 << " s(GeV^2)=" << lorentz_s/(CLHEP::GeV*
122 }
123
124 // For unit conversion
125 const G4double inv1e7 = 0.1/(CLHEP::GeV*CLHE
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 formu
132 // proton to nuclei Z^(2/3). The factor g4ca
133 // takes into account absorption of pi0 and
134
135 // pi- + p -> n + meson (0- pi0, 1- eta, 2-
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_
142 G4double sum = 0.0;
143 for (G4int i=0; i<5; ++i) {
144 G4double xg = std::max(1.0 + pG0[i] + pG
145 G4double xc = std::max(pC0[i] + pC1[i]*l
146 G4double xs = z23*piA[i]*g4calc->powA(x,
147 sum += xs;
148 fXSecPion[i] = sum;
149 }
150 result = sum*fact;
151 }
152
153 // pi+ + n -> p + meson (0- pi0, 1- eta, 2-
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_
160
161 // hydrogen target case Z = A = 1
162 // the cross section is defined by fractio
163 if (1 == Z) { n23 = ComputeDeuteronFractio
164 G4double sum = 0.0;
165 for (G4int i=0; i<5; ++i) {
166 G4double xg = std::max(1.0 + pG0[i] + pG
167 G4double xc = std::max(pC0[i] + pC1[i]*l
168 G4double xs = n23*piA[i]*g4calc->powA(x,
169 sum += xs;
170 fXSecPion[i] = sum;
171 }
172 result = sum*fact;
173 }
174
175 // Kaon x-sections depend on the primary par
176 // K- + p -> Kbar + n
177 else if (pdg == -321) {
178 G4double p_momentum = std::sqrt(pE*pE - pM
179 result = g4calc->Z23(Z)*g4calc->powA(p_mom
180 }
181
182 // K+ + n -> Kbar + p
183 else if (pdg == 321) {
184 G4double p_momentum = std::sqrt(pE*pE - pM
185 G4double n23 = g4calc->Z23(A-Z);
186 // hydrogen target case Z = A = 1
187 // the cross section is defined by fractio
188 if (1 == Z) { n23 = ComputeDeuteronFractio
189 result = n23*g4calc->powA(p_momentum, -1.6
190 }
191
192 // KL
193 else if (pdg == 130) {
194 // Cross section of KL = 0.5*(Cross sectio
195 const G4double p_momentum = std::sqrt(pE*p
196 result = 0.5*(g4calc->Z23(Z) + g4calc->Z23
197 g4calc->powA(p_momentum, -1.60)*kfact;
198 }
199 result *= fFactor;
200 if (verboseLevel > 1) {
201 G4cout << " Done for " << part->GetPart
202 << " res(mb)=" << result/CLHEP::millibar
203 }
204 return result;
205 }
206
207 const G4ParticleDefinition*
208 G4ChargeExchangeXS::SampleSecondaryType(const
209 const
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 an
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(co
253 {
254 for (auto const & elm : *mat->GetElementVect
255 if (1 == elm->GetZasInt()) {
256 G4double ab = 0.0;
257 const G4int nIso = (G4int)elm->GetNumber
258 const G4double* abu = elm->GetRelativeAb
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