| 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 // author: Vladimir.Grichine@cern.ch
27 //
28 // Implements data from: Barashenkov V.S., Nucleon-Nucleus Cross Section,
29 // Preprint JINR P2-89-770, p. 12, Dubna 1989 (scanned version from KEK)
30 // Based on G4NucleonNuclearCrossSection class
31 //
32 //
33
34 #include "G4ComponentBarNucleonNucleusXsc.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4DynamicParticle.hh"
37 #include "G4Neutron.hh"
38 #include "G4Proton.hh"
39 #include "G4Pow.hh"
40 #include "G4BarashenkovData.hh"
41 #include "G4IsotopeList.hh"
42 #include "G4HadronXSDataTable.hh"
43
44 ///////////////////////////////////////////////////////////////////////////////
45
46 G4double G4ComponentBarNucleonNucleusXsc::A75[93] = {0.0};
47 G4int G4ComponentBarNucleonNucleusXsc::theZ[] =
48 {2,4,6,7,8,11,13,14,20,26,29,42,48,50,74,82,92};
49 std::vector<G4PiData*>* G4ComponentBarNucleonNucleusXsc::thePData = nullptr;
50 std::vector<G4PiData*>* G4ComponentBarNucleonNucleusXsc::theNData = nullptr;
51
52 G4ComponentBarNucleonNucleusXsc::G4ComponentBarNucleonNucleusXsc()
53 : G4VComponentCrossSection("BarashenkovNucleonNucleusXsc")
54 {
55 theNeutron = G4Neutron::Neutron();
56 theProton = G4Proton::Proton();
57 if (nullptr == thePData) {
58 LoadData();
59 }
60 }
61
62 ////////////////////////////////////////////////////////////////////
63
64 G4double G4ComponentBarNucleonNucleusXsc::GetTotalIsotopeCrossSection(
65 const G4ParticleDefinition* aParticle,
66 G4double kinEnergy, G4int Z, G4int)
67 {
68 ComputeCrossSections(aParticle, kinEnergy, Z);
69 return fTotalXsc;
70 }
71
72 //////////////////////////////////////////////////////////////////////
73
74 G4double G4ComponentBarNucleonNucleusXsc::GetTotalElementCrossSection(
75 const G4ParticleDefinition* aParticle,
76 G4double kinEnergy, G4int Z, G4double)
77 {
78 ComputeCrossSections(aParticle, kinEnergy, Z);
79 return fTotalXsc;
80 }
81
82 ////////////////////////////////////////////////////////////////////
83
84 G4double G4ComponentBarNucleonNucleusXsc::GetInelasticIsotopeCrossSection(
85 const G4ParticleDefinition* aParticle,
86 G4double kinEnergy, G4int Z, G4int)
87 {
88 ComputeCrossSections(aParticle, kinEnergy, Z);
89 return fInelasticXsc;
90 }
91
92 /////////////////////////////////////////////////////////////////////
93
94 G4double G4ComponentBarNucleonNucleusXsc::GetInelasticElementCrossSection(
95 const G4ParticleDefinition* aParticle,
96 G4double kinEnergy, G4int Z, G4double)
97 {
98 ComputeCrossSections(aParticle, kinEnergy, Z);
99 return fInelasticXsc;
100 }
101
102 //////////////////////////////////////////////////////////////////
103
104 G4double G4ComponentBarNucleonNucleusXsc::GetElasticElementCrossSection(
105 const G4ParticleDefinition* aParticle,
106 G4double kinEnergy, G4int Z, G4double)
107 {
108 ComputeCrossSections(aParticle, kinEnergy, Z);
109 return fElasticXsc;
110 }
111
112 ///////////////////////////////////////////////////////////////////
113
114 G4double G4ComponentBarNucleonNucleusXsc::GetElasticIsotopeCrossSection(
115 const G4ParticleDefinition* aParticle,
116 G4double kinEnergy, G4int Z, G4int)
117 {
118 ComputeCrossSections(aParticle, kinEnergy, Z);
119 return fElasticXsc;
120 }
121
122 ////////////////////////////////////////////////////////////////////////////
123
124 void G4ComponentBarNucleonNucleusXsc::ComputeCrossSections(
125 const G4ParticleDefinition* aParticle, G4double kineticEnergy, G4int ZZ)
126 {
127 G4int Z = std::min(ZZ, 92);
128 G4int it = 0;
129 for(; it<NZ; ++it) { if(Z <= theZ[it]) { break; } }
130 if( it >= NZ ) { it = NZ-1; }
131
132 std::vector<G4PiData*>* theData = (aParticle == theNeutron) ? theNData : thePData;
133
134 if( theZ[it] == Z ) {
135 fInelasticXsc = (*theData)[it]->ReactionXSection(kineticEnergy);
136 fTotalXsc = (*theData)[it]->TotalXSection(kineticEnergy);
137 } else {
138 if(0 == it) { it = 1; }
139 G4double x1 = (*theData)[it-1]->ReactionXSection(kineticEnergy);
140 G4double xt1 = (*theData)[it-1]->TotalXSection(kineticEnergy);
141 G4double x2 = (*theData)[it]->ReactionXSection(kineticEnergy);
142 G4double xt2 = (*theData)[it]->TotalXSection(kineticEnergy);
143 G4int Z1 = theZ[it-1];
144 G4int Z2 = theZ[it];
145
146 fInelasticXsc = Interpolate(Z1, Z2, Z, x1, x2);
147 fTotalXsc = Interpolate(Z1, Z2, Z, xt1, xt2);
148 }
149
150 fElasticXsc = std::max(fTotalXsc - fInelasticXsc, 0.0);
151 }
152
153 /////////////////////////////////////////////////////////////////////////////
154
155 G4double G4ComponentBarNucleonNucleusXsc::
156 Interpolate(G4int Z1, G4int Z2, G4int Z, G4double x1, G4double x2) const
157 {
158 // for tabulated data, cross section scales with A^(2/3)
159 G4double r1 = x1* A75[Z] / A75[Z1];
160 G4double r2 = x2* A75[Z] / A75[Z2];
161 G4double alp1 = (aeff[Z] - aeff[Z1]);
162 G4double alp2 = (aeff[Z2] - aeff[Z]);
163 G4double result = (r1*alp2 + r2*alp1)/(alp1 + alp2);
164 // G4cout << "x1/2, z1/2 z" <<x1<<" "<<x2<<" "<<Z1<<" "<<Z2<<" "<<Z<<G4endl;
165 // G4cout << "res1/2 " << r1 <<" " << r2 <<" " << result<< G4endl;
166 return result;
167 }
168
169 /////////////////////////////////////////////////////////////////////////////
170
171 void G4ComponentBarNucleonNucleusXsc::Description(std::ostream& outFile) const
172 {
173 outFile << "G4ComponentBarNucleonNucleusXsc is a variant of the Barashenkov\n"
174 << "cross section parameterization to be used of protons and\n"
175 << "neutrons on targets heavier than hydrogen. It is intended for\n"
176 << "use as a cross section component and is currently used by\n"
177 << "G4BGGNucleonInelasticXS. It is valid for incident energies up\n"
178 << "to 1 TeV.\n";
179 }
180
181 /////////////////////////////////////////////////////////////////////////////
182
183 void G4ComponentBarNucleonNucleusXsc::LoadData()
184 {
185 theNData = new std::vector<G4PiData*>;
186 thePData = new std::vector<G4PiData*>;
187 theNData->resize(NZ, nullptr);
188 thePData->resize(NZ, nullptr);
189 auto ptr = G4HadronXSDataTable::Instance();
190 ptr->AddPiData(theNData);
191 ptr->AddPiData(thePData);
192
193 // He, Be, C
194 (*theNData)[0] = new G4PiData(he_m_t, he_m_in, e1, 44);
195 (*thePData)[0] = new G4PiData(he_m_t, he_p_in, e1, 44);
196
197 (*theNData)[1] = new G4PiData(be_m_t, be_m_in, e1, 44);
198 (*thePData)[1] = new G4PiData(be_m_t, be_p_in, e1, 44);
199
200 (*theNData)[2] = new G4PiData(c_m_t, c_m_in, e1, 44);
201 (*thePData)[2] = new G4PiData(c_m_t, c_p_in, e1, 44);
202
203 // N, O, Na
204 (*theNData)[3] = new G4PiData(n_m_t, n_m_in, e2, 44);
205 (*thePData)[3] = new G4PiData(n_m_t, n_p_in, e2, 44);
206
207 (*theNData)[4] = new G4PiData(o_m_t, o_m_in, e2, 44);
208 (*thePData)[4] = new G4PiData(o_m_t, o_p_in, e2, 44);
209
210 (*theNData)[5] = new G4PiData(na_m_t, na_m_in, e2, 44);
211 (*thePData)[5] = new G4PiData(na_m_t, na_p_in, e2, 44);
212
213 // Al, Si, Ca
214 (*theNData)[6] = new G4PiData(al_m_t, al_m_in, e3, 45);
215 (*thePData)[6] = new G4PiData(al_m_t, al_p_in, e3, 45);
216
217 (*theNData)[7] = new G4PiData(si_m_t, si_m_in, e3, 45);
218 (*thePData)[7] = new G4PiData(si_m_t, si_p_in, e3, 45);
219
220 (*theNData)[8] = new G4PiData(ca_m_t, ca_m_in, e3, 45);
221 (*thePData)[8] = new G4PiData(ca_m_t, ca_p_in, e3, 45);
222
223 // Fe, Cu, Mo
224 (*theNData)[9] = new G4PiData(fe_m_t, fe_m_in, e4, 47);
225 (*thePData)[9] = new G4PiData(fe_m_t, fe_p_in, e4, 47);
226
227 (*theNData)[10] = new G4PiData(cu_m_t, cu_m_in, e4, 47);
228 (*thePData)[10] = new G4PiData(cu_m_t, cu_p_in, e4, 47);
229
230 (*theNData)[11] = new G4PiData(mo_m_t, mo_m_in, e4, 47);
231 (*thePData)[11] = new G4PiData(mo_m_t, mo_p_in, e4, 47);
232
233 // Cd, Sn, W
234 (*theNData)[12] = new G4PiData(cd_m_t, cd_m_in, e5, 48);
235 (*thePData)[12] = new G4PiData(cd_m_t, cd_p_in, e5, 48);
236
237 (*theNData)[13] = new G4PiData(sn_m_t, sn_m_in, e5, 48);
238 (*thePData)[13] = new G4PiData(sn_m_t, sn_p_in, e5, 48);
239
240 (*theNData)[14] = new G4PiData(w_m_t, w_m_in, e5, 48);
241 (*thePData)[14] = new G4PiData(w_m_t, w_p_in, e5, 48);
242
243 // Pb, U
244 (*theNData)[15] = new G4PiData(pb_m_t, pb_m_in, e6, 46);
245 (*thePData)[15] = new G4PiData(pb_m_t, pb_p_in, e6, 46);
246
247 (*theNData)[16] = new G4PiData(u_m_t, u_m_in, e6, 46);
248 (*thePData)[16] = new G4PiData(u_m_t, u_p_in, e6, 46);
249
250 A75[0] = 1.0;
251 G4Pow* g4pow = G4Pow::GetInstance();
252 for(G4int i=1; i<93; ++i) {
253 A75[i] = g4pow->A23(aeff[i]); // interpolate by square ~ A^(2/3)
254 }
255 }
256
257 /////////////////////////////////////////////////////////////////////////////
258