ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BGGNucleonElasticXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BGGNucleonElasticXS.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: G4BGGNucleonElasticXS
32 //
33 // Author: Vladimir Ivanchenko
34 //
35 // Creation date: 13.03.2007
36 //
37 // -------------------------------------------------------------------
38 //
39 
40 #include "G4BGGNucleonElasticXS.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4ComponentGGHadronNucleusXsc.hh"
43 #include "G4NucleonNuclearCrossSection.hh"
44 #include "G4HadronNucleonXsc.hh"
45 #include "G4NuclearRadii.hh"
46 #include "G4Proton.hh"
47 #include "G4Neutron.hh"
48 #include "G4NistManager.hh"
49 #include "G4NuclearRadii.hh"
50 
51 #include "G4CrossSectionDataSetRegistry.hh"
52 
53 G4double G4BGGNucleonElasticXS::theGlauberFacP[93] = {0.0};
54 G4double G4BGGNucleonElasticXS::theCoulombFacP[93] = {0.0};
55 G4double G4BGGNucleonElasticXS::theGlauberFacN[93] = {0.0};
56 G4double G4BGGNucleonElasticXS::theCoulombFacN[93] = {0.0};
57 G4int G4BGGNucleonElasticXS::theA[93] = {0};
58 
59 #ifdef G4MULTITHREADED
60 G4Mutex G4BGGNucleonElasticXS::nucleonElasticXSMutex = G4MUTEX_INITIALIZER;
61 #endif
62 
63 G4BGGNucleonElasticXS::G4BGGNucleonElasticXS(const G4ParticleDefinition* p)
64  : G4VCrossSectionDataSet("BarashenkovGlauberGribov")
65 {
66  verboseLevel = 0;
67  fGlauberEnergy = 91.*GeV;
68  fLowEnergy = 14.0*MeV;
69  fNucleon = nullptr;
70  fGlauber = nullptr;
71  fHadron = nullptr;
72 
73  theProton= G4Proton::Proton();
74  isProton = (theProton == p);
75  isMaster = false;
76  SetForAllAtomsAndEnergies(true);
77 }
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 
81 G4BGGNucleonElasticXS::~G4BGGNucleonElasticXS()
82 {
83  delete fHadron;
84 }
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 
88 G4bool
89 G4BGGNucleonElasticXS::IsElementApplicable(const G4DynamicParticle*, G4int,
90  const G4Material*)
91 {
92  return true;
93 }
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96 
97 G4bool G4BGGNucleonElasticXS::IsIsoApplicable(const G4DynamicParticle*,
98  G4int Z, G4int,
99  const G4Element*,
100  const G4Material*)
101 {
102  return (1 == Z);
103 }
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
106 
107 G4double
108 G4BGGNucleonElasticXS::GetElementCrossSection(const G4DynamicParticle* dp,
109  G4int ZZ, const G4Material*)
110 {
111  // this method should be called only for Z > 1
112 
113  G4double cross = 0.0;
114  G4int Z = std::min(ZZ, 92);
115  G4double ekin = dp->GetKineticEnergy();
116  if(1 == Z) {
117  cross = 1.0115*GetIsoCrossSection(dp,1,1);
118  } else {
119  if(ekin <= fLowEnergy) {
120  cross = (isProton) ? theCoulombFacP[Z] : theCoulombFacN[Z];
121  cross *= CoulombFactor(ekin, Z);
122  } else if(ekin > fGlauberEnergy) {
123  cross = (isProton) ? theGlauberFacP[Z] : theGlauberFacN[Z];
124  cross *= fGlauber->GetElasticGlauberGribov(dp, Z, theA[Z]);
125  } else {
126  cross = fNucleon->GetElasticCrossSection(dp, Z);
127  }
128  }
129  if(verboseLevel > 1) {
130  G4cout << "G4BGGNucleonElasticXS::GetElementCrossSection for "
131  << dp->GetDefinition()->GetParticleName()
132  << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV
133  << " in nucleus Z= " << Z << " A= " << theA[Z]
134  << " XS(b)= " << cross/barn
135  << G4endl;
136  }
137  return cross;
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 
142 G4double
143 G4BGGNucleonElasticXS::GetIsoCrossSection(const G4DynamicParticle* dp,
144  G4int Z, G4int A,
145  const G4Isotope*,
146  const G4Element*,
147  const G4Material*)
148 {
149  // this method should be called only for Z = 1
150  fHadron->HadronNucleonXscNS(dp->GetDefinition(), theProton,
151  dp->GetKineticEnergy());
152  G4double cross = A*fHadron->GetElasticHadronNucleonXsc();
153 
154  if(verboseLevel > 1) {
155  G4cout << "G4BGGNucleonElasticXS::GetIsoCrossSection for "
156  << dp->GetDefinition()->GetParticleName()
157  << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV
158  << " in nucleus Z= " << Z << " A= " << A
159  << " XS(b)= " << cross/barn
160  << G4endl;
161  }
162  return cross;
163 }
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
166 
167 void G4BGGNucleonElasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
168 {
169  if(fNucleon) { return; }
170  if(&p == theProton || &p == G4Neutron::Neutron()) {
171  isProton = (theProton == &p);
172 
173  } else {
174  G4ExceptionDescription ed;
175  ed << "This BGG cross section is applicable only to nucleons and not to "
176  << p.GetParticleName() << G4endl;
177  G4Exception("G4BGGNucleonElasticXS::BuildPhysicsTable", "had001",
178  FatalException, ed);
179  return;
180  }
181 
182  fNucleon = new G4NucleonNuclearCrossSection();
183  fGlauber = new G4ComponentGGHadronNucleusXsc();
184  fHadron = new G4HadronNucleonXsc();
185 
186  fNucleon->BuildPhysicsTable(p);
187 
188  if(0 == theA[0]) {
189 #ifdef G4MULTITHREADED
190  G4MUTEXLOCK(&nucleonElasticXSMutex);
191  if(0 == theA[0]) {
192 #endif
193  isMaster = true;
194 #ifdef G4MULTITHREADED
195  }
196  G4MUTEXUNLOCK(&nucleonElasticXSMutex);
197 #endif
198  } else {
199  return;
200  }
201 
202  if(isMaster && 0 == theA[0]) {
203 
204  theA[0] = theA[1] = 1;
205  G4ThreeVector mom(0.0,0.0,1.0);
206  G4DynamicParticle dp(theProton, mom, fGlauberEnergy);
207 
208  G4NistManager* nist = G4NistManager::Instance();
209  G4double csup, csdn;
210 
211  if(verboseLevel > 0) {
212  G4cout << "### G4BGGNucleonElasticXS::Initialise for "
213  << p.GetParticleName() << G4endl;
214  }
215 
216  for(G4int iz=2; iz<93; ++iz) {
217  G4int A = G4lrint(nist->GetAtomicMassAmu(iz));
218  theA[iz] = A;
219 
220  csup = fGlauber->GetElasticGlauberGribov(&dp, iz, A);
221  csdn = fNucleon->GetElasticCrossSection(&dp, iz);
222  theGlauberFacP[iz] = csdn/csup;
223  }
224 
225  dp.SetDefinition(G4Neutron::Neutron());
226  for(G4int iz=2; iz<93; ++iz) {
227  csup = fGlauber->GetElasticGlauberGribov(&dp, iz, theA[iz]);
228  csdn = fNucleon->GetElasticCrossSection(&dp, iz);
229  theGlauberFacN[iz] = csdn/csup;
230 
231  if(verboseLevel > 0) {
232  G4cout << "Z= " << iz << " A= " << theA[iz]
233  << " GFactorP= " << theGlauberFacP[iz]
234  << " GFactorN= " << theGlauberFacN[iz] << G4endl;
235  }
236  }
237 
238  theCoulombFacP[0] = theCoulombFacP[1] =
239  theCoulombFacN[0] = theCoulombFacN[1] = 1.0;
240  dp.SetDefinition(theProton);
241  dp.SetKineticEnergy(fLowEnergy);
242  for(G4int iz=2; iz<93; ++iz) {
243  theCoulombFacP[iz] = fNucleon->GetElasticCrossSection(&dp, iz)
244  /CoulombFactor(fLowEnergy, iz);
245  }
246  dp.SetDefinition(G4Neutron::Neutron());
247  for(G4int iz=2; iz<93; ++iz) {
248  theCoulombFacN[iz] = fNucleon->GetElasticCrossSection(&dp, iz)
249  /CoulombFactor(fLowEnergy, iz);
250 
251  if(verboseLevel > 0) {
252  G4cout << "Z= " << iz << " A= " << theA[iz]
253  << " CFactorP= " << theCoulombFacP[iz]
254  << " CFactorN= " << theCoulombFacN[iz] << G4endl;
255  }
256  }
257  }
258 }
259 
260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261 
262 G4double G4BGGNucleonElasticXS::CoulombFactor(G4double kinEnergy, G4int Z)
263 {
264  G4double res= 1.0;
265  if(isProton) {
266  res = G4NuclearRadii::CoulombFactor(Z, theA[Z], theProton, kinEnergy);
267  }
268  return res;
269 }
270 
271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
272 
273 void G4BGGNucleonElasticXS::CrossSectionDescription(std::ostream& outFile) const
274 {
275  outFile << "The Barashenkov-Glauber-Gribov cross section handles elastic\n"
276  << "scattering of protons and neutrons from nuclei using the\n"
277  << "Barashenkov parameterization below 91 GeV and the Glauber-Gribov\n"
278  << "parameterization above 91 GeV. n";
279 }
280 
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......