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