ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ComponentBarNucleonNucleusXsc.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ComponentBarNucleonNucleusXsc.cc
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 "G4NistManager.hh"
42 
44 
45 G4double G4ComponentBarNucleonNucleusXsc::theA[93] = {0.0};
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 #ifdef G4MULTITHREADED
53 G4Mutex G4ComponentBarNucleonNucleusXsc::barNNXSMutex = G4MUTEX_INITIALIZER;
54 #endif
55 
56 G4ComponentBarNucleonNucleusXsc::G4ComponentBarNucleonNucleusXsc()
57  : G4VComponentCrossSection("BarashenkovNucleonNucleusXsc"),
58  fTotalXsc(0.0), fInelasticXsc(0.0), fElasticXsc(0.0), isMaster(false)
59 {
60  theNeutron = G4Neutron::Neutron();
61  theProton = G4Proton::Proton();
62 }
63 
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;
77  }
78 }
79 
81 
82 G4double G4ComponentBarNucleonNucleusXsc::GetTotalIsotopeCrossSection(
83  const G4ParticleDefinition* aParticle,
84  G4double kinEnergy, G4int Z, G4int)
85 {
86  ComputeCrossSections(aParticle, kinEnergy, Z);
87  return fTotalXsc;
88 }
89 
91 
92 G4double G4ComponentBarNucleonNucleusXsc::GetTotalElementCrossSection(
93  const G4ParticleDefinition* aParticle,
94  G4double kinEnergy, G4int Z, G4double)
95 {
96  ComputeCrossSections(aParticle, kinEnergy, Z);
97  return fTotalXsc;
98 }
99 
101 
102 G4double G4ComponentBarNucleonNucleusXsc::GetInelasticIsotopeCrossSection(
103  const G4ParticleDefinition* aParticle,
104  G4double kinEnergy, G4int Z, G4int)
105 {
106  ComputeCrossSections(aParticle, kinEnergy, Z);
107  return fInelasticXsc;
108 }
109 
111 
112 G4double G4ComponentBarNucleonNucleusXsc::GetInelasticElementCrossSection(
113  const G4ParticleDefinition* aParticle,
114  G4double kinEnergy, G4int Z, G4double)
115 {
116  ComputeCrossSections(aParticle, kinEnergy, Z);
117  return fInelasticXsc;
118 }
119 
121 
122 G4double G4ComponentBarNucleonNucleusXsc::GetElasticElementCrossSection(
123  const G4ParticleDefinition* aParticle,
124  G4double kinEnergy, G4int Z, G4double)
125 {
126  ComputeCrossSections(aParticle, kinEnergy, Z);
127  return fElasticXsc;
128 }
129 
131 
132 G4double G4ComponentBarNucleonNucleusXsc::GetElasticIsotopeCrossSection(
133  const G4ParticleDefinition* aParticle,
134  G4double kinEnergy, G4int Z, G4int)
135 {
136  ComputeCrossSections(aParticle, kinEnergy, Z);
137  return fElasticXsc;
138 }
139 
141 
142 void G4ComponentBarNucleonNucleusXsc::ComputeCrossSections(
143  const G4ParticleDefinition* aParticle, G4double kineticEnergy, G4int ZZ)
144 {
145  G4int Z = std::min(ZZ, 92);
146  G4int it = 0;
147  for(; it<NZ; ++it) { if(Z <= theZ[it]) { break; } }
148  if( it >= NZ ) { it = NZ-1; }
149 
150  std::vector<G4PiData*>* theData = (aParticle == theNeutron) ? theNData : thePData;
151 
152  if( theZ[it] == Z ) {
153  fInelasticXsc = (*theData)[it]->ReactionXSection(kineticEnergy);
154  fTotalXsc = (*theData)[it]->TotalXSection(kineticEnergy);
155  } else {
156  if(0 == it) { it = 1; }
157  G4double x1 = (*theData)[it-1]->ReactionXSection(kineticEnergy);
158  G4double xt1 = (*theData)[it-1]->TotalXSection(kineticEnergy);
159  G4double x2 = (*theData)[it]->ReactionXSection(kineticEnergy);
160  G4double xt2 = (*theData)[it]->TotalXSection(kineticEnergy);
161  G4int Z1 = theZ[it-1];
162  G4int Z2 = theZ[it];
163 
164  fInelasticXsc = Interpolate(Z1, Z2, Z, x1, x2);
165  fTotalXsc = Interpolate(Z1, Z2, Z, xt1, xt2);
166  }
167 
168  fElasticXsc = std::max(fTotalXsc - fInelasticXsc, 0.0);
169 }
170 
172 
173 G4double G4ComponentBarNucleonNucleusXsc::
174 Interpolate(G4int Z1, G4int Z2, G4int Z, G4double x1, G4double x2) const
175 {
176  // for tabulated data, cross section scales with A^(2/3)
177  G4double r1 = x1* A75[Z] / A75[Z1];
178  G4double r2 = x2* A75[Z] / A75[Z2];
179  G4double alp1 = (theA[Z] - theA[Z1]);
180  G4double alp2 = (theA[Z2] - theA[Z]);
181  G4double result = (r1*alp2 + r2*alp1)/(alp1 + alp2);
182  // G4cout << "x1/2, z1/2 z" <<x1<<" "<<x2<<" "<<Z1<<" "<<Z2<<" "<<Z<<G4endl;
183  // G4cout << "res1/2 " << r1 <<" " << r2 <<" " << result<< G4endl;
184  return result;
185 }
186 
188 
189 void G4ComponentBarNucleonNucleusXsc::Description(std::ostream& outFile) const
190 {
191  outFile << "G4ComponentBarNucleonNucleusXsc is a variant of the Barashenkov\n"
192  << "cross section parameterization to be used of protons and\n"
193  << "neutrons on targets heavier than hydrogen. It is intended for\n"
194  << "use as a cross section component and is currently used by\n"
195  << "G4BGGNucleonInelasticXS. It is valid for incident energies up\n"
196  << "to 1 TeV.\n";
197 }
198 
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 
219 
220 void G4ComponentBarNucleonNucleusXsc::LoadData()
221 {
222  theNData = new std::vector<G4PiData*>;
223  thePData = new std::vector<G4PiData*>;
224  theNData->resize(NZ, nullptr);
225  thePData->resize(NZ, nullptr);
226 
227  // He, Be, C
228  (*theNData)[0] = new G4PiData(he_m_t, he_m_in, e1, 44);
229  (*thePData)[0] = new G4PiData(he_m_t, he_p_in, e1, 44);
230 
231  (*theNData)[1] = new G4PiData(be_m_t, be_m_in, e1, 44);
232  (*thePData)[1] = new G4PiData(be_m_t, be_p_in, e1, 44);
233 
234  (*theNData)[2] = new G4PiData(c_m_t, c_m_in, e1, 44);
235  (*thePData)[2] = new G4PiData(c_m_t, c_p_in, e1, 44);
236 
237  // N, O, Na
238  (*theNData)[3] = new G4PiData(n_m_t, n_m_in, e2, 44);
239  (*thePData)[3] = new G4PiData(n_m_t, n_p_in, e2, 44);
240 
241  (*theNData)[4] = new G4PiData(o_m_t, o_m_in, e2, 44);
242  (*thePData)[4] = new G4PiData(o_m_t, o_p_in, e2, 44);
243 
244  (*theNData)[5] = new G4PiData(na_m_t, na_m_in, e2, 44);
245  (*thePData)[5] = new G4PiData(na_m_t, na_p_in, e2, 44);
246 
247  // Al, Si, Ca
248  (*theNData)[6] = new G4PiData(al_m_t, al_m_in, e3, 45);
249  (*thePData)[6] = new G4PiData(al_m_t, al_p_in, e3, 45);
250 
251  (*theNData)[7] = new G4PiData(si_m_t, si_m_in, e3, 45);
252  (*thePData)[7] = new G4PiData(si_m_t, si_p_in, e3, 45);
253 
254  (*theNData)[8] = new G4PiData(ca_m_t, ca_m_in, e3, 45);
255  (*thePData)[8] = new G4PiData(ca_m_t, ca_p_in, e3, 45);
256 
257  // Fe, Cu, Mo
258  (*theNData)[9] = new G4PiData(fe_m_t, fe_m_in, e4, 47);
259  (*thePData)[9] = new G4PiData(fe_m_t, fe_p_in, e4, 47);
260 
261  (*theNData)[10] = new G4PiData(cu_m_t, cu_m_in, e4, 47);
262  (*thePData)[10] = new G4PiData(cu_m_t, cu_p_in, e4, 47);
263 
264  (*theNData)[11] = new G4PiData(mo_m_t, mo_m_in, e4, 47);
265  (*thePData)[11] = new G4PiData(mo_m_t, mo_p_in, e4, 47);
266 
267  // Cd, Sn, W
268  (*theNData)[12] = new G4PiData(cd_m_t, cd_m_in, e5, 48);
269  (*thePData)[12] = new G4PiData(cd_m_t, cd_p_in, e5, 48);
270 
271  (*theNData)[13] = new G4PiData(sn_m_t, sn_m_in, e5, 48);
272  (*thePData)[13] = new G4PiData(sn_m_t, sn_p_in, e5, 48);
273 
274  (*theNData)[14] = new G4PiData(w_m_t, w_m_in, e5, 48);
275  (*thePData)[14] = new G4PiData(w_m_t, w_p_in, e5, 48);
276 
277  // Pb, U
278  (*theNData)[15] = new G4PiData(pb_m_t, pb_m_in, e6, 46);
279  (*thePData)[15] = new G4PiData(pb_m_t, pb_p_in, e6, 46);
280 
281  (*theNData)[16] = new G4PiData(u_m_t, u_m_in, e6, 46);
282  (*thePData)[16] = new G4PiData(u_m_t, u_p_in, e6, 46);
283 
284  G4NistManager* nist = G4NistManager::Instance();
285  A75[0] = theA[0] = 1.0;
286  G4Pow* g4pow = G4Pow::GetInstance();
287  for(G4int i=1; i<93; ++i) {
288  theA[i] = nist->GetAtomicMassAmu(i);
289  A75[i] = g4pow->A23(theA[i]); // interpolate by square ~ A^(2/3)
290  }
291 }
292