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