ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BGGNucleonInelasticXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BGGNucleonInelasticXS.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 // -------------------------------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4BGGNucleonInelasticXS
32 //
33 // Author: Vladimir Ivanchenko
34 //
35 // Creation date: 13.03.2007
36 // Modifications:
37 //
38 //
39 // -------------------------------------------------------------------
40 //
41 
42 #include "G4BGGNucleonInelasticXS.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4ComponentGGHadronNucleusXsc.hh"
45 #include "G4NucleonNuclearCrossSection.hh"
46 #include "G4HadronNucleonXsc.hh"
47 #include "G4ComponentSAIDTotalXS.hh"
48 #include "G4Proton.hh"
49 #include "G4Neutron.hh"
50 #include "G4NistManager.hh"
51 #include "G4Material.hh"
52 #include "G4Element.hh"
53 #include "G4Isotope.hh"
54 #include "G4Log.hh"
55 #include "G4Exp.hh"
56 #include "G4NuclearRadii.hh"
57 
58 #include "G4CrossSectionDataSetRegistry.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 const G4double llog10 = G4Log(10.);
63 
64 G4double G4BGGNucleonInelasticXS::theGlauberFacP[93] = {0.0};
65 G4double G4BGGNucleonInelasticXS::theCoulombFacP[93] = {0.0};
66 G4double G4BGGNucleonInelasticXS::theGlauberFacN[93] = {0.0};
67 G4double G4BGGNucleonInelasticXS::theCoulombFacN[93] = {0.0};
68 G4int G4BGGNucleonInelasticXS::theA[93] = {0};
69 
70 #ifdef G4MULTITHREADED
71 G4Mutex G4BGGNucleonInelasticXS::nucleonInelasticXSMutex = G4MUTEX_INITIALIZER;
72 #endif
73 
74 G4BGGNucleonInelasticXS::G4BGGNucleonInelasticXS(const G4ParticleDefinition* p)
75  : G4VCrossSectionDataSet("BarashenkovGlauberGribov")
76 {
77  verboseLevel = 0;
78  fGlauberEnergy = 91.*GeV;
79  fLowEnergy = 14.*MeV;
80 
81  fNucleon = nullptr;
82  fGlauber = nullptr;
83  fHadron = nullptr;
84 
85  theProton= G4Proton::Proton();
86  isProton = (theProton == p);
87  isMaster = false;
88  SetForAllAtomsAndEnergies(true);
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 
93 G4BGGNucleonInelasticXS::~G4BGGNucleonInelasticXS()
94 {
95  delete fHadron;
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 
100 G4bool G4BGGNucleonInelasticXS::IsElementApplicable(const G4DynamicParticle*,
101  G4int, const G4Material*)
102 {
103  return true;
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 
108 G4bool G4BGGNucleonInelasticXS::IsIsoApplicable(const G4DynamicParticle*,
109  G4int Z, G4int,
110  const G4Element*,
111  const G4Material*)
112 {
113  return (1 == Z);
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
118 G4double
119 G4BGGNucleonInelasticXS::GetElementCrossSection(const G4DynamicParticle* dp,
120  G4int ZZ, const G4Material*)
121 {
122  G4double cross = 0.0;
123  G4double ekin = dp->GetKineticEnergy();
124  G4int Z = std::min(ZZ, 92);
125  if(1 == Z) {
126  cross = 1.0115*GetIsoCrossSection(dp,1,1);
127  } else {
128  if(ekin <= fLowEnergy) {
129  cross = (isProton) ? theCoulombFacP[Z] : theCoulombFacN[Z];
130  cross *= CoulombFactor(ekin, Z);
131  } else if(ekin > fGlauberEnergy) {
132  cross = (isProton) ? theGlauberFacP[Z] : theGlauberFacN[Z];
133  cross *= fGlauber->GetInelasticGlauberGribov(dp, Z, theA[Z]);
134  } else {
135  cross = fNucleon->GetElementCrossSection(dp, Z);
136  }
137  }
138 
139  if(verboseLevel > 1) {
140  G4cout << "G4BGGNucleonInelasticXS::GetCrossSection for "
141  << dp->GetDefinition()->GetParticleName()
142  << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV
143  << " in nucleus Z= " << Z << " A= " << theA[Z]
144  << " XS(b)= " << cross/barn
145  << G4endl;
146  }
147  return cross;
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
151 
152 G4double
153 G4BGGNucleonInelasticXS::GetIsoCrossSection(const G4DynamicParticle* dp,
154  G4int Z, G4int A,
155  const G4Isotope*,
156  const G4Element*,
157  const G4Material*)
158 {
159  // this method should be called only for Z = 1
160  fHadron->HadronNucleonXscNS(dp->GetDefinition(), theProton,
161  dp->GetKineticEnergy());
162  G4double cross = A*fHadron->GetInelasticHadronNucleonXsc();
163 
164  if(verboseLevel > 1) {
165  G4cout << "G4BGGNucleonInelasticXS::GetIsoCrossSection for "
166  << dp->GetDefinition()->GetParticleName()
167  << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV
168  << " in nucleus Z= " << Z << " A= " << theA[Z]
169  << " XS(b)= " << cross/barn
170  << G4endl;
171  }
172  return cross;
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176 
177 void G4BGGNucleonInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
178 {
179  if(fNucleon) { return; }
180  if(&p == theProton || &p == G4Neutron::Neutron()) {
181  isProton = (theProton == &p);
182  } else {
183  G4ExceptionDescription ed;
184  ed << "This BGG cross section is applicable only to nucleons and not to "
185  << p.GetParticleName() << G4endl;
186  G4Exception("G4BGGNucleonInelasticXS::BuildPhysicsTable", "had001",
187  FatalException, ed);
188  return;
189  }
190 
191  fNucleon = new G4NucleonNuclearCrossSection();
192  fGlauber = new G4ComponentGGHadronNucleusXsc();
193  fHadron = new G4HadronNucleonXsc();
194 
195  fNucleon->BuildPhysicsTable(p);
196 
197  if(0 == theA[0]) {
198 #ifdef G4MULTITHREADED
199  G4MUTEXLOCK(&nucleonInelasticXSMutex);
200  if(0 == theA[0]) {
201 #endif
202  isMaster = true;
203 #ifdef G4MULTITHREADED
204  }
205  G4MUTEXUNLOCK(&nucleonInelasticXSMutex);
206 #endif
207  } else {
208  return;
209  }
210 
211  if(isMaster && 0 == theA[0]) {
212 
213  theA[0] = theA[1] = 1;
214  G4ThreeVector mom(0.0,0.0,1.0);
215  G4DynamicParticle dp(theProton, mom, fGlauberEnergy);
216 
217  G4NistManager* nist = G4NistManager::Instance();
218  G4double csup, csdn;
219 
220  if(verboseLevel > 0) {
221  G4cout << "### G4BGGNucleonInelasticXS::Initialise for "
222  << p.GetParticleName() << G4endl;
223  }
224  for(G4int iz=2; iz<93; ++iz) {
225 
226  G4int A = G4lrint(nist->GetAtomicMassAmu(iz));
227  theA[iz] = A;
228 
229  csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, A);
230  csdn = fNucleon->GetElementCrossSection(&dp, iz);
231  theGlauberFacP[iz] = csdn/csup;
232  }
233 
234  dp.SetDefinition(G4Neutron::Neutron());
235  for(G4int iz=2; iz<93; ++iz) {
236  csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, theA[iz]);
237  csdn = fNucleon->GetElementCrossSection(&dp, iz);
238  theGlauberFacN[iz] = csdn/csup;
239 
240  if(verboseLevel > 0) {
241  G4cout << "Z= " << iz << " A= " << theA[iz]
242  << " GFactorP= " << theGlauberFacP[iz]
243  << " GFactorN= " << theGlauberFacN[iz] << G4endl;
244  }
245  }
246  theCoulombFacP[1] = theCoulombFacN[1] = 1.0;
247  dp.SetDefinition(theProton);
248  dp.SetKineticEnergy(fLowEnergy);
249  for(G4int iz=2; iz<93; ++iz) {
250  theCoulombFacP[iz] = fNucleon->GetElementCrossSection(&dp, iz)
251  /CoulombFactor(fLowEnergy, iz);
252  }
253  dp.SetDefinition(G4Neutron::Neutron());
254  for(G4int iz=2; iz<93; ++iz) {
255  theCoulombFacN[iz] = fNucleon->GetElementCrossSection(&dp, iz)
256  /CoulombFactor(fLowEnergy, iz);
257 
258  if(verboseLevel > 0) {
259  G4cout << "Z= " << iz << " A= " << theA[iz]
260  << " CFactorP= " << theCoulombFacP[iz]
261  << " CFactorN= " << theCoulombFacN[iz] << G4endl;
262  }
263  }
264  }
265 }
266 
267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
268 
269 G4double G4BGGNucleonInelasticXS::CoulombFactor(G4double kinEnergy, G4int Z)
270 {
271  G4double res = 0.0;
272 
273  if(kinEnergy <= 0.0) { return res; }
274 
275  G4double elog = G4Log(kinEnergy/GeV)/llog10;
276  G4double aa = theA[Z];
277 
278  if(isProton) {
279 
280  res = G4NuclearRadii::CoulombFactor(Z, theA[Z], theProton, kinEnergy);
281 
282  // from G4ProtonInelasticCrossSection
283  if(res > 0.0) {
284  G4double ff1 = 5.6 - 0.016*aa; // slope of the drop at medium energies.
285  G4double ff2 = 1.37 + 1.37/aa; // start of the slope.
286  G4double ff3 = 0.8 + 18./aa - 0.002*aa; // stephight
287  res *= (1.0 + ff3*(1.0 - (1.0/(1+G4Exp(-ff1*(elog + ff2))))));
288  ff1 = 8. - 8./aa - 0.008*aa; // slope of the rise
289  ff2 = 2.34 - 5.4/aa - 0.0028*aa; // start of the rise
290  res /= (1.0 + G4Exp(-ff1*(elog + ff2)));
291  }
292  } else {
293  // from G4NeutronInelasticCrossSection
294  G4double p3 = 0.6 + 13./aa - 0.0005*aa;
295  G4double p4 = 7.2449 - 0.018242*aa;
296  G4double p5 = 1.36 + 1.8/aa + 0.0005*aa;
297  G4double p6 = 1. + 200./aa + 0.02*aa;
298  G4double p7 = 3.0 - (aa-70.)*(aa-200.)/11000.;
299 
300  G4double firstexp = G4Exp(-p4*(elog + p5));
301  G4double secondexp = G4Exp(-p6*(elog + p7));
302 
303  res = (1.+p3*firstexp/(1. + firstexp))/(1. + secondexp);
304  }
305  return res;
306 }
307 
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309 
310 void G4BGGNucleonInelasticXS::CrossSectionDescription(std::ostream& outFile) const
311 {
312  outFile << "The Barashenkov-Glauber-Gribov cross section calculates inelastic\n"
313  << "scattering of protons and neutrons from nuclei using the\n"
314  << "Barashenkov parameterization below 91 GeV and the Glauber-Gribov\n"
315  << "parameterization above 91 GeV. It uses the G4HadronNucleonXsc\n"
316  << "cross section component for hydrogen targets, and the\n"
317  << "G4ComponentGGHadronNucleusXsc component for other targets.\n";
318 }
319 
320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......