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