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