| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // ------------------------------------------- 26 // -------------------------------------------------------------------
27 // 27 //
28 // GEANT4 Class file 28 // GEANT4 Class file
29 // 29 //
30 // 30 //
31 // File name: G4ChargeExchangeXS 31 // File name: G4ChargeExchangeXS
32 // 32 //
33 33
34 #include "G4ChargeExchangeXS.hh" 34 #include "G4ChargeExchangeXS.hh"
35 #include "G4DynamicParticle.hh" 35 #include "G4DynamicParticle.hh"
36 #include "G4ElementTable.hh" 36 #include "G4ElementTable.hh"
37 #include "G4Material.hh" 37 #include "G4Material.hh"
38 #include "G4Element.hh" 38 #include "G4Element.hh"
39 #include "G4IsotopeList.hh" << 39 #include "G4Isotope.hh"
40 #include "G4HadronicParameters.hh" 40 #include "G4HadronicParameters.hh"
41 #include "Randomize.hh" 41 #include "Randomize.hh"
42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
43 #include "G4NucleiProperties.hh" 43 #include "G4NucleiProperties.hh"
44 #include "G4Pow.hh" 44 #include "G4Pow.hh"
45 45
46 #include "G4PionZero.hh" 46 #include "G4PionZero.hh"
47 #include "G4Eta.hh" 47 #include "G4Eta.hh"
48 #include "G4KaonZeroLong.hh" 48 #include "G4KaonZeroLong.hh"
49 #include "G4KaonZeroShort.hh" 49 #include "G4KaonZeroShort.hh"
50 #include "G4KaonPlus.hh" 50 #include "G4KaonPlus.hh"
51 #include "G4KaonMinus.hh" 51 #include "G4KaonMinus.hh"
52 #include "G4ParticleTable.hh" 52 #include "G4ParticleTable.hh"
53 53
54 namespace { 54 namespace {
55 // V. Lyubovitsky parameterisation 55 // V. Lyubovitsky parameterisation
56 const G4double piA[5] = {430., 36., 1.37, 56 const G4double piA[5] = {430., 36., 1.37, 2.0, 60.}; // A
57 const G4double pAP[5] = {1.04, 1.26, 1.35, 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, 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, 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, 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., 61 const G4double pG1[5] = {-0.23, -0.5, 0., 0., -2.}; // g1
62 62
>> 63 // beta_prime value for calculation of cross section of pi0 and eta
>> 64 // absorption inside different nuclei
63 const G4double beta_prime_pi = 0.0410; 65 const G4double beta_prime_pi = 0.0410;
>> 66 const G4double beta_prime_eta = 0.0402;
64 } 67 }
65 68
>> 69
66 G4ChargeExchangeXS::G4ChargeExchangeXS() 70 G4ChargeExchangeXS::G4ChargeExchangeXS()
67 { 71 {
68 if (verboseLevel > 1) { 72 if (verboseLevel > 1) {
69 G4cout << "G4ChargeExchangeXS::G4ChargeEx 73 G4cout << "G4ChargeExchangeXS::G4ChargeExchangeXS" << G4endl;
70 } 74 }
71 g4calc = G4Pow::GetInstance(); 75 g4calc = G4Pow::GetInstance();
72 auto table = G4ParticleTable::GetParticleTab 76 auto table = G4ParticleTable::GetParticleTable();
73 const G4String nam[5] = {"pi0", "eta", "eta_ 77 const G4String nam[5] = {"pi0", "eta", "eta_prime", "omega", "f2(1270)"};
74 for (G4int i=0; i<5; ++i) { 78 for (G4int i=0; i<5; ++i) {
75 fPionSecPD[i] = table->FindParticle(nam[i] 79 fPionSecPD[i] = table->FindParticle(nam[i]);
76 if (nullptr == fPionSecPD[i]) { 80 if (nullptr == fPionSecPD[i]) {
77 G4ExceptionDescription ed; 81 G4ExceptionDescription ed;
78 ed << "### meson " << nam[i] << " is not 82 ed << "### meson " << nam[i] << " is not found out in the particle table";
79 G4Exception("G4ChargeExchangeXS::G4Charg 83 G4Exception("G4ChargeExchangeXS::G4ChargeExchangeXS()","had044",
80 FatalException, ed,""); 84 FatalException, ed,"");
81 } 85 }
82 } 86 }
83 } 87 }
84 88
85 // Print the information of this .cc file 89 // Print the information of this .cc file
86 void G4ChargeExchangeXS::CrossSectionDescripti 90 void G4ChargeExchangeXS::CrossSectionDescription(std::ostream& outFile) const
87 { 91 {
88 outFile << "G4ChargeExchangeXS calculates ch 92 outFile << "G4ChargeExchangeXS calculates charge exchange cross section for "
89 << "pi+, pi-, K+, K-, KL\n"; 93 << "pi+, pi-, K+, K-, KL\n";
90 } 94 }
91 95
92 G4bool G4ChargeExchangeXS::IsElementApplicable << 96 G4bool G4ChargeExchangeXS::IsIsoApplicable(const G4DynamicParticle*,
93 << 97 G4int, G4int,
>> 98 const G4Element*, const G4Material*)
94 { 99 {
95 return true; 100 return true;
96 } 101 }
97 102
98 G4double 103 G4double
99 G4ChargeExchangeXS::GetElementCrossSection(con << 104 G4ChargeExchangeXS::GetIsoCrossSection(const G4DynamicParticle* aParticle,
100 G4int ZZ, const G4Material* << 105 G4int Z, G4int A,
>> 106 const G4Isotope*, const G4Element*,
>> 107 const G4Material*)
101 { 108 {
102 G4double result = 0.0; 109 G4double result = 0.0;
103 const G4double pE = aParticle->GetTotalEnerg 110 const G4double pE = aParticle->GetTotalEnergy();
104 if (pE <= fEnergyLimit) { return result; } << 111 if (pE <= fEnergyLimit) { return result; }
105 <<
106 auto part = aParticle->GetDefinition(); 112 auto part = aParticle->GetDefinition();
107 G4int pdg = part->GetPDGEncoding(); 113 G4int pdg = part->GetPDGEncoding();
108 114
109 // Get or calculate the proton mass, particl << 115 // Get or calculate the nucleus mass, particle mass,particle kinetic energy
110 G4double tM = CLHEP::proton_mass_c2; << 116 // and particle total energy
111 G4double pM = part->GetPDGMass(); << 117 G4double tM = G4NucleiProperties::GetNuclearMass(A, Z);
>> 118 G4double pM = part->GetPDGMass();
>> 119
>> 120 // Calculate s(lorentz invariant)
112 G4double lorentz_s = tM*tM + 2*tM*pE + pM*p 121 G4double lorentz_s = tM*tM + 2*tM*pE + pM*pM;
113 if (lorentz_s <= (tM + pM)*(tM + pM)) { retu 122 if (lorentz_s <= (tM + pM)*(tM + pM)) { return result; }
114 123
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 124 // For unit conversion
125 const G4double inv1e7 = 0.1/(CLHEP::GeV*CLHE << 125 const G4double inv1e7 = 1e-7;
126 const G4double fact = 1e-30*CLHEP::cm2; 126 const G4double fact = 1e-30*CLHEP::cm2;
127 const G4double pfact = 0.1/CLHEP::GeV; 127 const G4double pfact = 0.1/CLHEP::GeV;
128 const G4double kfact = 56.3*fact; 128 const G4double kfact = 56.3*fact;
129 const G4double csmax = 1e-16; << 129
>> 130 G4double logA = g4calc->logZ(A);
130 131
131 // The approximation of Glauber-Gribov formu 132 // The approximation of Glauber-Gribov formula -> extend it from interaction with
132 // proton to nuclei Z^(2/3). The factor g4ca << 133 // proton to nuclei Z^(2/3). The factor g4calc->powA(A,-beta_prime_pi*G4Log(A))
133 // takes into account absorption of pi0 and 134 // takes into account absorption of pi0 and eta
134 << 135 // pi- + p -> sum of (pi0 + eta) + n
135 // pi- + p -> n + meson (0- pi0, 1- eta, 2- <<
136 if (pdg == -211) { 136 if (pdg == -211) {
137 G4double z23 = g4calc->Z23(Z); << 137 const G4double z23 = g4calc->Z23(Z);
138 G4double x = lorentz_s*inv1e7; << 138 const G4int z = A/2;
139 G4double logX = G4Log(x); << 139 const G4double a23 = g4calc->Z23(z);
140 G4double logA = g4calc->logZ(A); << 140 const G4double x = lorentz_s*inv1e7;
141 G4double xf = g4calc->powZ(A, -beta_prime_ << 141 G4double sum = 122.*z23*g4calc->powA(x, -1.23)*g4calc->powZ(A,-beta_prime_pi*logA);
142 G4double sum = 0.0; << 142 fXSecPion[0] = sum;
143 for (G4int i=0; i<5; ++i) { << 143 sum += 31.*z23*g4calc->powA(x, -1.53)*g4calc->powZ(A,-beta_prime_eta*logA);
144 G4double xg = std::max(1.0 + pG0[i] + pG << 144 fXSecPion[1] = sum;
145 G4double xc = std::max(pC0[i] + pC1[i]*l << 145 const G4double logX = G4Log(x);
146 G4double xs = z23*piA[i]*g4calc->powA(x, << 146 for (G4int i=2; i<5; ++i) {
147 sum += xs; << 147 sum += piA[i]*z23*g4calc->powA(x, -pAP[i])*(1.0 + pG0[i] + pG1[i]*logX)
148 fXSecPion[i] = sum; << 148 *g4calc->powA(z23, -0.15*a23)/(pC0[i] + pC1[i]*logX);
>> 149 fXSecPion[i] = sum;
149 } 150 }
150 result = sum*fact; 151 result = sum*fact;
151 } 152 }
152 153
153 // pi+ + n -> p + meson (0- pi0, 1- eta, 2- << 154 // pi+ + n -> sum of (pi0 + eta) + p
154 else if (pdg == 211) { 155 else if (pdg == 211) {
155 G4double n23 = g4calc->Z23(A - Z); << 156 const G4double n23 = g4calc->Z23(A - Z);
156 G4double x = lorentz_s*inv1e7; << 157 const G4int z = A/2;
157 G4double logX = G4Log(x); << 158 const G4double a23 = g4calc->Z23(z);
158 G4double logA = g4calc->logZ(A); << 159 const G4double x = lorentz_s*inv1e7;
159 G4double xf = g4calc->powZ(A, -beta_prime_ << 160 G4double sum = 122.*n23*g4calc->powA(x, -1.23)*g4calc->powZ(A,-beta_prime_pi*logA);
160 << 161 fXSecPion[0] = sum;
161 // hydrogen target case Z = A = 1 << 162 sum += 31.*n23*g4calc->powA(x, -1.53)*g4calc->powZ(A,-beta_prime_eta*logA);
162 // the cross section is defined by fractio << 163 fXSecPion[1] = sum;
163 if (1 == Z) { n23 = ComputeDeuteronFractio << 164 const G4double logX = G4Log(x);
164 G4double sum = 0.0; << 165 for (G4int i=2; i<5; ++i) {
165 for (G4int i=0; i<5; ++i) { << 166 sum += piA[i]*n23*g4calc->powA(x, -pAP[i])*(1.0 + pG0[i] + pG1[i]*logX)
166 G4double xg = std::max(1.0 + pG0[i] + pG << 167 *g4calc->powA(n23, -0.15*a23)/(pC0[i] + pC1[i]*logX);
167 G4double xc = std::max(pC0[i] + pC1[i]*l << 168 fXSecPion[i] = sum;
168 G4double xs = n23*piA[i]*g4calc->powA(x, <<
169 sum += xs; <<
170 fXSecPion[i] = sum; <<
171 } 169 }
172 result = sum*fact; 170 result = sum*fact;
173 } 171 }
174 172
175 // Kaon x-sections depend on the primary par <<
176 // K- + p -> Kbar + n 173 // K- + p -> Kbar + n
177 else if (pdg == -321) { << 174 else if (pdg == -321){
178 G4double p_momentum = std::sqrt(pE*pE - pM << 175 // Calculate the momentum of the bombarding particles and convert
>> 176 // it to GeV/c^2 unit
>> 177 const G4double p_momentum = std::sqrt(pE*pE - pM*pM)*pfact;
179 result = g4calc->Z23(Z)*g4calc->powA(p_mom 178 result = g4calc->Z23(Z)*g4calc->powA(p_momentum, -1.60)*kfact;
180 } 179 }
181 180
182 // K+ + n -> Kbar + p 181 // K+ + n -> Kbar + p
183 else if (pdg == 321) { 182 else if (pdg == 321) {
184 G4double p_momentum = std::sqrt(pE*pE - pM << 183 const G4double p_momentum = std::sqrt(pE*pE - pM*pM)*pfact;
185 G4double n23 = g4calc->Z23(A-Z); << 184 result = g4calc->Z23(A-Z)*g4calc->powA(p_momentum, -1.60)*kfact;
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 } 185 }
191 186
192 // KL 187 // KL
193 else if (pdg == 130) { 188 else if (pdg == 130) {
194 // Cross section of KL = 0.5*(Cross sectio << 189 // Cross section of K-long = 0.5*(Cross section of K+ + Cross section of K-)
195 const G4double p_momentum = std::sqrt(pE*p 190 const G4double p_momentum = std::sqrt(pE*pE - pM*pM)*pfact;
196 result = 0.5*(g4calc->Z23(Z) + g4calc->Z23 191 result = 0.5*(g4calc->Z23(Z) + g4calc->Z23(A-Z))*
197 g4calc->powA(p_momentum, -1.60)*kfact; 192 g4calc->powA(p_momentum, -1.60)*kfact;
198 } 193 }
199 result *= fFactor; << 194
200 if (verboseLevel > 1) { << 195 return result*fFactor;
201 G4cout << " Done for " << part->GetPart <<
202 << " res(mb)=" << result/CLHEP::millibar <<
203 } <<
204 return result; <<
205 } 196 }
206 197
207 const G4ParticleDefinition* 198 const G4ParticleDefinition*
208 G4ChargeExchangeXS::SampleSecondaryType(const 199 G4ChargeExchangeXS::SampleSecondaryType(const G4ParticleDefinition* part,
209 const 200 const G4int Z, const G4int A)
210 { 201 {
211 const G4ParticleDefinition* pd = nullptr; 202 const G4ParticleDefinition* pd = nullptr;
212 G4int pdg = std::abs(part->GetPDGEncoding()) << 203 G4int pdg = part->GetPDGEncoding();
213 204
214 // pi- + p / pi+ + n 205 // pi- + p / pi+ + n
215 if (pdg == 211) { << 206 if (std::abs(pdg) == 211) {
216 pd = fPionSecPD[0]; << 207 const G4double x = fXSecPion[4]*G4UniformRand();
217 G4double x = fXSecPion[4]*G4UniformRand(); <<
218 for (G4int i=0; i<5; ++i) { 208 for (G4int i=0; i<5; ++i) {
219 if (x <= fXSecPion[i]) { 209 if (x <= fXSecPion[i]) {
220 pd = fPionSecPD[i]; << 210 return fPionSecPD[i];
221 break; <<
222 } 211 }
223 } 212 }
224 } 213 }
225 214
226 // K- + p / K+ + n 215 // K- + p / K+ + n
227 // Equal opportunity of producing k-short an 216 // Equal opportunity of producing k-short and k-long
228 else if (pdg == 321) { << 217 else if (std::abs(pdg) == 321) {
229 if (G4UniformRand() >= 0.5) { << 218 if (G4UniformRand() > 0.5) {
230 pd = G4KaonZeroLong::KaonZeroLong(); << 219 pd = G4KaonZeroLong::KaonZeroLong();
231 } 220 }
232 else { 221 else {
233 pd = G4KaonZeroShort::KaonZeroShort(); 222 pd = G4KaonZeroShort::KaonZeroShort();
234 } 223 }
235 } 224 }
236 225
237 // KL + nucleus << 226 // KL + atom
238 else if (pdg == 130) { << 227 else if (std::abs(pdg) == 130) {
239 G4double prob = (G4double)Z/(G4double)A; 228 G4double prob = (G4double)Z/(G4double)A;
240 if (G4UniformRand() >= prob) { << 229 if (G4UniformRand() > prob) {
241 pd = G4KaonMinus::KaonMinus(); 230 pd = G4KaonMinus::KaonMinus();
242 } 231 }
243 else { 232 else {
244 pd = G4KaonPlus::KaonPlus(); 233 pd = G4KaonPlus::KaonPlus();
245 } 234 }
246 } 235 }
247 236
248 return pd; 237 return pd;
249 } 238 }
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 239