Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BGGNucleonInelasticXS.cc
Go to the documentation of this file.
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 // $Id$
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4BGGNucleonInelasticXS
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 13.03.2007
38 // Modifications:
39 //
40 //
41 // -------------------------------------------------------------------
42 //
43 
44 #include "G4BGGNucleonInelasticXS.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4GlauberGribovCrossSection.hh"
47 #include "G4NucleonNuclearCrossSection.hh"
48 #include "G4HadronNucleonXsc.hh"
49 //#include "G4HadronInelasticDataSet.hh"
50 #include "G4ComponentSAIDTotalXS.hh"
51 #include "G4Proton.hh"
52 #include "G4Neutron.hh"
53 #include "G4NistManager.hh"
54 #include "G4Material.hh"
55 #include "G4Element.hh"
56 #include "G4Isotope.hh"
57 
58 #include "G4CrossSectionDataSetRegistry.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 G4BGGNucleonInelasticXS::G4BGGNucleonInelasticXS(const G4ParticleDefinition* p)
63  : G4VCrossSectionDataSet("Barashenkov-Glauber")
64 {
65  verboseLevel = 0;
66  fGlauberEnergy = 91.*GeV;
67  fLowEnergy = 14.*MeV;
68  fHighEnergy = 5.*GeV;
69  fSAIDHighEnergyLimit = 1.3*GeV;
70  for (G4int i = 0; i < 93; ++i) {
71  theGlauberFac[i] = 0.0;
72  theCoulombFac[i] = 0.0;
73  theA[i] = 1;
74  }
75  fNucleon = 0;
76  fGlauber = 0;
77  fHadron = 0;
78  // fGHEISHA = 0;
79  fSAID = 0;
80  particle = p;
81  theProton= G4Proton::Proton();
82  isProton = false;
83  isInitialized = false;
84 }
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 
88 G4BGGNucleonInelasticXS::~G4BGGNucleonInelasticXS()
89 {
90  delete fHadron;
91  delete fSAID;
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 
96 G4bool G4BGGNucleonInelasticXS::IsElementApplicable(const G4DynamicParticle*,
97  G4int, const G4Material*)
98 {
99  return true;
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 
104 G4bool G4BGGNucleonInelasticXS::IsIsoApplicable(const G4DynamicParticle*,
105  G4int Z, G4int A,
106  const G4Element*,
107  const G4Material*)
108 {
109  return (1 == Z && 2 >= A);
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113 
114 G4double
115 G4BGGNucleonInelasticXS::GetElementCrossSection(const G4DynamicParticle* dp,
116  G4int ZZ, const G4Material*)
117 {
118  // this method should be called only for Z > 1
119 
120  G4double cross = 0.0;
121  G4double ekin = dp->GetKineticEnergy();
122  G4int Z = ZZ;
123  if(1 == Z) {
124  cross = 1.0115*GetIsoCrossSection(dp,1,1);
125  } else if(2 == Z) {
126  if(ekin > fGlauberEnergy) {
127  cross = theGlauberFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, theA[Z]);
128  } else {
129  cross = fNucleon->GetElementCrossSection(dp, Z);
130  }
131 
132  } else {
133  if(Z > 92) { Z = 92; }
134 
135  if(ekin <= fLowEnergy) {
136  cross = theCoulombFac[Z]*CoulombFactor(ekin, Z);
137  } else if(ekin > fGlauberEnergy) {
138  cross = theGlauberFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, theA[Z]);
139  } else {
140  cross = fNucleon->GetElementCrossSection(dp, Z);
141  }
142  }
143 
144  if(verboseLevel > 1) {
145  G4cout << "G4BGGNucleonInelasticXS::GetCrossSection for "
146  << dp->GetDefinition()->GetParticleName()
147  << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV
148  << " in nucleus Z= " << Z << " A= " << theA[Z]
149  << " XS(b)= " << cross/barn
150  << G4endl;
151  }
152  return cross;
153 }
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156 
157 G4double
158 G4BGGNucleonInelasticXS::GetIsoCrossSection(const G4DynamicParticle* dp,
159  G4int Z, G4int A,
160  const G4Isotope*,
161  const G4Element*,
162  const G4Material*)
163 {
164  // this method should be called only for Z = 1
165 
166  G4double cross = 0.0;
167  G4double ekin = dp->GetKineticEnergy();
168 
169  if(ekin <= fSAIDHighEnergyLimit) {
170  cross = fSAID->GetInelasticIsotopeCrossSection(particle, ekin, 1, 1);
171  } else if(ekin < fHighEnergy) {
172  fHadron->GetHadronNucleonXscNS(dp, theProton);
173  cross = (theCoulombFac[0]/ekin + 1)*fHadron->GetInelasticHadronNucleonXsc();
174  } else {
175  fHadron->GetHadronNucleonXscPDG(dp, theProton);
176  cross = (theCoulombFac[1]/ekin + 1)*fHadron->GetInelasticHadronNucleonXsc();
177  }
178  cross *= A;
179 
180  if(verboseLevel > 1) {
181  G4cout << "G4BGGNucleonInelasticXS::GetCrossSection for "
182  << dp->GetDefinition()->GetParticleName()
183  << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV
184  << " in nucleus Z= " << Z << " A= " << A
185  << " XS(b)= " << cross/barn
186  << G4endl;
187  }
188  return cross;
189 }
190 
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
192 
193 void G4BGGNucleonInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
194 {
195  if(&p == theProton || &p == G4Neutron::Neutron()) {
196  particle = &p;
197  } else {
198  G4cout << "### G4BGGNucleonInelasticXS WARNING: is not applicable to "
199  << p.GetParticleName()
200  << G4endl;
201  throw G4HadronicException(__FILE__, __LINE__,
202  "G4BGGNucleonElasticXS::BuildPhysicsTable is used for wrong particle");
203  return;
204  }
205 
206  if(isInitialized) { return; }
207  isInitialized = true;
208 
209  fNucleon = (G4NucleonNuclearCrossSection*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4NucleonNuclearCrossSection::Default_Name());
210  fGlauber = (G4GlauberGribovCrossSection*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4GlauberGribovCrossSection::Default_Name());
211 
212  fHadron = new G4HadronNucleonXsc();
213  //fGHEISHA = new G4HadronInelasticDataSet();
214  fSAID = new G4ComponentSAIDTotalXS();
215 
216  fNucleon->BuildPhysicsTable(*particle);
217  fGlauber->BuildPhysicsTable(*particle);
218 
219  if(particle == theProton) {
220  isProton = true;
221  fSAIDHighEnergyLimit = 2*GeV;
222  fHighEnergy = 2*GeV;
223  }
224 
225  G4ParticleDefinition* part = const_cast<G4ParticleDefinition*>(particle);
226  G4ThreeVector mom(0.0,0.0,1.0);
227  G4DynamicParticle dp(part, mom, fGlauberEnergy);
228 
229  G4NistManager* nist = G4NistManager::Instance();
230  G4int A;
231 
232  G4double csup, csdn;
233 
234  if(verboseLevel > 0) {
235  G4cout << "### G4BGGNucleonInelasticXS::Initialise for "
236  << particle->GetParticleName() << G4endl;
237  }
238  for(G4int iz=2; iz<93; iz++) {
239 
240  A = G4lrint(nist->GetAtomicMassAmu(iz));
241  theA[iz] = A;
242 
243  csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, A);
244  csdn = fNucleon->GetElementCrossSection(&dp, iz);
245 
246  theGlauberFac[iz] = csdn/csup;
247  if(verboseLevel > 0) {
248  G4cout << "Z= " << iz << " A= " << A
249  << " GlauberFactor= " << theGlauberFac[iz] << G4endl;
250  }
251  }
252  //const G4Material* mat = 0;
253 
254  dp.SetKineticEnergy(fSAIDHighEnergyLimit);
255  fHadron->GetHadronNucleonXscNS(&dp, theProton);
256  theCoulombFac[0] = fSAIDHighEnergyLimit*
257  (fSAID->GetInelasticIsotopeCrossSection(particle,fSAIDHighEnergyLimit,1,1)
258  /fHadron->GetInelasticHadronNucleonXsc() - 1);
259 
260  //G4cout << "Z=1 E(GeV)= " << fSAIDHighEnergyLimit/GeV
261  // << " xsNS(b)= " << fHadron->GetInelasticHadronNucleonXsc()/barn;
262  fHadron->GetHadronNucleonXscPDG(&dp, theProton);
263  //G4cout << " xsPDG(b)= " << fHadron->GetInelasticHadronNucleonXsc()/barn;
264  //G4cout << " xsSAID(b)= " << fSAID->GetInelasticIsotopeCrossSection(particle,fSAIDHighEnergyLimit,1,1)/barn << G4endl;
265 
266  dp.SetKineticEnergy(fHighEnergy);
267  fHadron->GetHadronNucleonXscPDG(&dp, theProton);
268  G4double x = fHadron->GetInelasticHadronNucleonXsc();
269 
270  //G4cout << "Z=1 E(GeV)= " << fHighEnergy/GeV
271  // << " xsPDG(b)= " << fHadron->GetInelasticHadronNucleonXsc()/barn;
272 
273  fHadron->GetHadronNucleonXscNS(&dp, theProton);
274  theCoulombFac[1] = fHighEnergy*((theCoulombFac[0]/fHighEnergy + 1)
275  *fHadron->GetInelasticHadronNucleonXsc()/x - 1);
276 
277  fHadron->GetHadronNucleonXscNS(&dp, theProton);
278  //G4cout << " xsNS(b)= " << fHadron->GetInelasticHadronNucleonXsc()/barn << G4endl;
279 
280  if(verboseLevel > 0) {
281  G4cout << "Z=1 A=1" << " CoulombFactor[0]= " << theCoulombFac[0]
282  << " CoulombFactor[1]= " << theCoulombFac[1] << G4endl;
283  }
284  theCoulombFac[2] = 1.0;
285 
286  dp.SetKineticEnergy(fLowEnergy);
287  for(G4int iz=3; iz<93; iz++) {
288  theCoulombFac[iz] =
289  fNucleon->GetElementCrossSection(&dp, iz)/CoulombFactor(fLowEnergy, iz);
290 
291  if(verboseLevel > 0) {
292  G4cout << "Z= " << iz << " A= " << theA[iz]
293  << " CoulombFactor= " << theCoulombFac[iz] << G4endl;
294  }
295  }
296 }
297 
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
299 
300 G4double G4BGGNucleonInelasticXS::CoulombFactor(G4double kinEnergy, G4int Z)
301 {
302  G4double res= 0.0;
303  if(kinEnergy <= 0.0) { return res; }
304  else if (Z <= 1) { return kinEnergy*kinEnergy; }
305 
306  G4double elog = std::log10(kinEnergy/GeV);
307  G4double aa = theA[Z];
308 
309  // from G4ProtonInelasticCrossSection
310  if(isProton) {
311 
312  G4double ff1 = 0.70 - 0.002*aa; // slope of the drop at medium energies.
313  G4double ff2 = 1.00 + 1/aa; // start of the slope.
314  G4double ff3 = 0.8 + 18/aa - 0.002*aa; // stephight
315  res = 1.0 + ff3*(1.0 - (1.0/(1+std::exp(-8*ff1*(elog + 1.37*ff2)))));
316 
317  ff1 = 1. - 1./aa - 0.001*aa; // slope of the rise
318  ff2 = 1.17 - 2.7/aa-0.0014*aa; // start of the rise
319  res /= (1 + std::exp(-8.*ff1*(elog + 2*ff2)));
320 
321  } else {
322 
323  // from G4NeutronInelasticCrossSection
324  G4double p3 = 0.6 + 13./aa - 0.0005*aa;
325  G4double p4 = 7.2449 - 0.018242*aa;
326  G4double p5 = 1.36 + 1.8/aa + 0.0005*aa;
327  G4double p6 = 1. + 200./aa + 0.02*aa;
328  G4double p7 = 3.0 - (aa-70.)*(aa-200.)/11000.;
329 
330  G4double firstexp = std::exp(-p4*(elog + p5));
331  G4double secondexp = std::exp(-p6*(elog + p7));
332 
333  res = (1.+p3*firstexp/(1. + firstexp))/(1. + secondexp);
334 
335  }
336  return res;
337 }
338 
339 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
340 
341 void G4BGGNucleonInelasticXS::CrossSectionDescription(std::ostream& outFile) const
342 {
343  outFile << "The Barashenkov-Glauber-Gribov cross section calculates inelastic\n"
344  << "scattering of protons and neutrons from nuclei using the\n"
345  << "Barashenkov parameterization below 91 GeV and the Glauber-Gribov\n"
346  << "parameterization above 91 GeV. It uses the G4HadronNucleonXsc\n"
347  << "cross section component for hydrogen targets, and the\n"
348  << "G4GlauberGribovCrossSection component for other targets.\n";
349 }
350 
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Generated on Sat May 25 2013 14:33:40 for Geant4 by   doxygen 1.8.4